O metal está sujeito a uma temperatura constante na extremidade x=0, onde T(0, t)=1. Enquanto na outra extremidade, em x=1, o fluxo de calor é nulo. Portanto:
No instante inicial, a temperatura da barra é uniforme T(x, 0)=0. Encontre o perfil de temperatura do sólido utilizando o método das diferenças finitas centradas para o domínio espacial e o método de Runge-Kutta 4 (RK4) para a integração temporal. Você deve precisar de pelo menos n=51 nós de discretização e um passo de integração h=10^{-4}. Você deve plotar em um único gráfico as curvas de t=0.0,0.025,0.05, \ldots, 0.25 (portanto, 11 perfis de temperatura).
Observação: eu decidi tratar apenas os nós internos j=1,\dots,49 como incógnitas para integrar no tempo; T_0 é fixo e T_{50} é obtido por T_{50}=T_{49}.
Passo 3
Discretização da segunda derivada (diferença finita centrada) e sistema de ODEs:
Para um nó interior genérico j (com 1\le j\le49) a segunda derivada espacial é aproximada por
Aqui \mathbf{f} é o procedimento que, dado um vetor \mathbf{u} (valores em j=1\ldots49), computa as derivadas \dfrac{dT_j}{dt} usando as fórmulas do Passo 3. Sempre que for preciso T_0 uso 1 e para T_{50} uso T_{49} (impondo Neumann).
RK4 é explícito: cada avaliação de \mathbf{f} combina vizinhos imediatos, o que é eficiente para este problema.
Passo 5
Condições iniciais, passos de tempo e instantes de salvamento:
Condição inicial para o vetor incógnita:
T_j(0)=0\quad\text{para }j=1,\dots,49.
Passo temporal pedido: h=10^{-4}.
O tempo final que precisamos cobrir é t_{\max}=0{,}25. O número de passos é
Instantes para plotar são t=0{,}0,0{,}025,0{,}05,\dots,0{,}25. Cada intervalo de 0{,}025 corresponde a 0{,}025/10^{-4}=250 passos. Então devo salvar o perfil a cada 250 passos, inclusive no passo 0.
Ao salvar um perfil, reconstruo o vetor completo de temperaturas em todos os nós j=0,\dots,50 por:
Para esquemas explícitos aplicados à equação do calor, a razão adimensional
\lambda=\frac{h}{\Delta x^2}
é relevante para estabilidade. Aqui \Delta x^2=(0{,}02)^2=4\times10^{-4}. Logo
\lambda=\frac{10^{-4}}{4\times10^{-4}}=0{,}25.
Esse valor é seguro para integração explícita (RK4) neste problema — está abaixo de 0{,}5 que costuma aparecer como limite para esquemas simples; RK4 tem região de estabilidade compatível com este \lambda. Portanto a escolha h=10^{-4} é adequada.
Passo 7
Detalhes práticos de implementação (pseudocódigo em Python + dicas):
se j=49: dT49/dt = (T_{48} - T_{49})/dx^2 (usando T50=T49).
Loop em n=0…2499:
calcular k1=f(u)
calcular k2=f(u + h/2*k1)
calcular k3=f(u + h/2*k2)
calcular k4=f(u + h*k3)
atualizar u <- u + h/6*(k1 + 2k2 + 2k3 + k4)
a cada 250 passos salvar perfil completo: T_full=[1, u[0], u[1],…,u[48], u[48]] (observe T50=T49)
Depois plotar os 11 perfis salvos num único gráfico com legenda indicando o tempo.
Pseudocódigo (esquemático):
n = 51
dx = 1.0/(n-1) # 0.02
m = n-2 # 49 incógnitas (T1...T49)
h = 1e-4
Nsteps = int(0.25/h) # 2500
save_every = 250
# inicializa
u = zeros(m) # T1..T49 = 0
profiles = []
times = []
def rhs(u):
dudt = zeros(m)
# j=1 -> i=0
dudt[0] = (u[1] - 2*u[0] + 1)/(dx*dx)
# 2<=j<=48 -> i=1..46
for i in range(1, m-1):
dudt[i] = (u[i+1] - 2*u[i] + u[i-1])/(dx*dx)
# j=49 -> i=m-1
dudt[m-1] = (u[m-2] - u[m-1])/(dx*dx)
return dudt
# salvar t=0
T_full = [1.0] + list(u) + [u[-1]]
profiles.append(T_full)
times.append(0.0)
t = 0.0
for step in range(1, Nsteps+1):
k1 = rhs(u)
k2 = rhs(u + 0.5*h*k1)
k3 = rhs(u + 0.5*h*k2)
k4 = rhs(u + h*k3)
u = u + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)
t += h
if step % save_every == 0:
T_full = [1.0] + list(u) + [u[-1]]
profiles.append(T_full)
times.append(t)
# depois usar matplotlib para plotar x vs cada perfil em profiles,
# com legendas times.
Observações de implementação:
Use arrays NumPy para eficiência.
Ao montar rhs, cuidado com índices: u[0] corresponde a T_1, u[m-1] a T_{49}.
Ao plotar, crie vetor x = np.linspace(0,1,n) e plote cada profiles[k] contra x.
Configure labels, título e legenda com os tempos, e use cores distintas.
Passo 8
Como plotar e o que esperar visualmente:
Para plot no Python com matplotlib:
Após obter profiles (lista de 11 vetores cada um de tamanho 51) e times (lista de 11 tempos), faça:
x = np.linspace(0,1,51)
Para cada perfil p, tt em zip(profiles,times) faça plt.plot(x,p,label=f"t={tt:.3f}")
plt.legend(), plt.xlabel("x"), plt.ylabel("T"), plt.title("Perfis de temperatura"), plt.show().
O comportamento físico esperado:
Inicialmente T(x,0)=0 em toda a barra exceto x=0 que está fixo em 1.
Com o tempo o calor se difunde para a direita: os perfis mostram T decrescendo de 1 em x=0 até valores menores em direção a x=1.
A condição de Neumann em x=1 garante fluxo nulo no extremo direito. No regime estacionário a solução de \dfrac{d^2T}{dx^2}=0 com T(0)=1 e T'(1)=0 é T(x)\equiv1. Assim, com tempo a barra tende para temperatura uniforme 1.
As curvas salvas em t=0{,}0{,}025,\dots,0{,}25 mostrarão progressiva elevação de temperatura até se aproximarem de 1.
Passo 9
Verificações rápidas e dicas para a prova:
Cheque sempre:
índice e correspondência entre vetor numérico e nós físicos (off-by-one é a fonte comum de erro);
aplicação correta de T_0=1 e da condição T_{50}=T_{49};
o cálculo de \Delta x e uso de \Delta x^2 no denominador;
que os snapshots estão sendo salvos nos tempos corretos (a cada 250 passos).
Dica para a prova: escreva explicitamente as equações discretas para os três casos: nó adjacente a Dirichlet, nó interior genérico, nó adjacente a Neumann — isso mostra entendimento e evita erros.
Dica numérica: se você implementar e observar instabilidade (oscilações explosivas), reduza h até que \lambda=h/\Delta x^2 esteja bem dentro da região estável.
Passo 10
Encerramento motivacional e garantia de que com isso você terá o gráfico:
Seguindo o passo a passo acima (montagem de \mathbf{f}, RK4, salvamento a cada 250 passos) você obterá os 11 perfis pedidos e conseguirá plotá-los num único gráfico. Se quiser, eu posso te passar o código completo para copiar e executar (Python/NumPy/Matplotlib) ou adaptar para MATLAB — me diga qual prefere e eu te dou o código pronto.
Resposta
Os 11 perfis de temperatura T(x,t) em t=0{,}0{,}025,0{,}05,\ldots,0{,}25 são obtidos numericamente com n=51 e h=10^{-4} executando o algoritmo descrito (diferenças centradas espaciais, sistema de 49 ODEs para T_1,\dots,T_{49}, RK4 no tempo, impondo T_0=1 e T_{50}=T_{49}); ao executar o pseudocódigo fornecido você terá o gráfico único com todas as curvas.