Considere a equação diferencial parcial para a condução de calor em uma barra metálica em regime transiente:
\frac{\partial T}{\partial t}=\frac{\partial^{2} T}{\partial x^{2}}
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:
\left.\frac{\partial T}{\partial x}\right|_{x=1}=0
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).

Passo 1

O enunciado fornece:
\frac{\partial T}{\partial t}=\frac{\partial^{2}T}{\partial x^{2}}
em 0\le x\le1 e t\ge0. Condições de contorno e inicial:
T(0,t)=1,
\left.\frac{\partial T}{\partial x}\right|_{x=1}=0,
T(x,0)=0.
O que é pedido:
  • Obter numericamente os perfis de temperatura T(x,t) usando diferenças finitas centradas no espaço e RK4 no tempo.
  • Usar pelo menos n=51 nós espaciais e passo temporal h=10^{-4}.
  • Plotar num único gráfico as curvas em t=0{,}0,0{,}025,0{,}05,\ldots,0{,}25 (11 perfis).
Anotei todas as informações importantes acima e quais curvas precisamos salvar.

Passo 2

Discretização espacial (definição dos nós e do espaçamento):
  • Vou tomar n=51 nós uniformes no intervalo [0,1]. O espaçamento é
\Delta x=\frac{1-0}{n-1}=\frac{1}{50}=0{,}02.
  • Enumero os nós por j=0,1,\ldots,50. O nó j tem coordenada x_j=j\Delta x.
  • A temperatura no nó j e tempo t será escrita como T_j(t). A condição de contorno T(0,t)=1 corresponde a T_0(t)=1 para todo t.
  • A condição de Neumann em x=1 é \partial_x T(1,t)=0. Vou impor isso pela aproximação por diferença:
\left.\frac{\partial T}{\partial x}\right|_{x=1}\approx\frac{T_{50}-T_{49}}{\Delta x}=0\quad\Rightarrow\quad T_{50}=T_{49}.
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
\frac{\partial^2 T}{\partial x^2}\bigg|_{x_j}\approx\frac{T_{j+1}-2T_j+T_{j-1}}{\Delta x^2}.
  • Então as equações diferenciais ordinárias para os nós internos ficam
\frac{dT_j}{dt}=\frac{T_{j+1}-2T_j+T_{j-1}}{\Delta x^2},\qquad j=1,\dots,49.
  • No caso j=49 usamos a condição T_{50}=T_{49}. Substituindo, a equação para j=49 fica
\frac{dT_{49}}{dt}=\frac{T_{49}-2T_{49}+T_{48}}{\Delta x^2}=\frac{T_{48}-T_{49}}{\Delta x^2}.
  • Para j=1 a equação usa T_0=1:
\frac{dT_1}{dt}=\frac{T_{2}-2T_{1}+T_{0}}{\Delta x^2}=\frac{T_{2}-2T_{1}+1}{\Delta x^2}.
  • Assim temos um sistema de 49 ODEs para o vetor
\mathbf{u}(t)=(T_1(t),T_2(t),\dots,T_{49}(t))^\top.
  • Escrevo o lado direito como uma função \mathbf{f}(\mathbf{u}) definida componente a componente conforme as fórmulas acima.

Passo 4

Explicação do método de integração temporal RK4 aplicado ao sistema:
  • RK4 é um integrador explícito de quarto ordem para ODEs. Para um sistema \dfrac{d\mathbf{u}}{dt}=\mathbf{f}(\mathbf{u},t) com passo h os estágios são:
\begin{aligned} \mathbf{k}_1&=\mathbf{f}(\mathbf{u}_n,t_n),\\\ \mathbf{k}_2&=\mathbf{f}\!\left(\mathbf{u}_n+\tfrac{h}{2}\mathbf{k}_1,\,t_n+\tfrac{h}{2}\right),\\\ \mathbf{k}_3&=\mathbf{f}\!\left(\mathbf{u}_n+\tfrac{h}{2}\mathbf{k}_2,\,t_n+\tfrac{h}{2}\right),\\\ \mathbf{k}_4&=\mathbf{f}\!\left(\mathbf{u}_n+h\mathbf{k}_3,\,t_n+h\right), \end{aligned}
e a atualização é
\mathbf{u}_{n+1}=\mathbf{u}_n+\frac{h}{6}\left(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4\right).
  • 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 é
N_{\text{steps}}=\frac{t_{\max}}{h}=\frac{0{,}25}{10^{-4}}=2500.
  • 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:
T_0=1,\quad T_j=\text{(da solução)}\\ (j=1,\dots,49),\quad T_{50}=T_{49}.çã

Passo 6

Condição de estabilidade e escolha do passo h:
  • 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):
  • Estrutura geral:
    1. Montar vetor x_j=j\Delta x para j=0,\ldots,50.
    2. Inicializar vetor incógnita \mathbf{u}=(T_1,\dots,T_{49})=\mathbf{0}.
    3. Definir função rhs(u) que retorna vetor de tamanho 49 com as derivadas:
      • Para índice interno i correspondente a nó j=i+1 (porque i=0 corresponde a j=1):
        • se j=1: dT1/dt = (T2 - 2 T1 + 1)/dx^2.
        • se 2\le j\le48: dTj/dt = (T_{j+1} - 2 T_j + T_{j-1})/dx^2.
        • se j=49: dT49/dt = (T_{48} - T_{49})/dx^2 (usando T50=T49).
    4. 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)
    5. 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.

Exercícios relacionados