Matemática Numérica II

6 Equações Diferenciais Parciais 6.1 Equação de Poisson 6.3 Equação da Onda

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

6.2 Equação do Calor

Consideramos a equação do calor com condição inicial dada e condições de contorno de Dirichlet homogêneas


	$\displaystyle u_{t}-\alpha u_{xx}=f(t,x),\leavevmode\nobreak\ 0<t\leq t_{f},% \leavevmode\nobreak\ a<x<b,$		(6.22a)
	$\displaystyle u(0,x)=u_{0}(x),\leavevmode\nobreak\ a<x<b,$		(6.22b)
	$\displaystyle u(t,a)=u(t,b)=0,\leavevmode\nobreak\ 0<t\leq t_{f},$		(6.22c)

onde $u=u(t,x)$ é a incógnita.

O problema (6.22) é um problema de valor inicial com condições de contorno. Uma das estratégias numéricas de solução é o chamado Método das Linhas, o qual trata separadamente as discretizações espacial e temporal. Aqui, vamos começar pela discretização espacial e, então, trataremos a discretização temporal.

1. Discretização Espacial.

Na discretização espacial, aplicamos o Método de Diferenças Finitas (MDF). Começamos considerando uma malha uniforme de nodos $x_{i}=a+(i-1)h_{x}$ , $i=1,2,\dotsc,n_{x}+1$ , com tamanho de malha $h_{x}=(b-a)/n_{x}$ , sendo $n_{x}$ o número de subintervalos. Denotando $u_{i}(t)\approx u(t,x_{i})$ e empregando a fórmula de diferenças finitas centrais $D^{2}_{0,h^{2}}$ , temos que a Eq. (6.22a) fica aproximada por

\frac{du_{i}}{dt}=\frac{\alpha}{h_{x}^{2}}u_{i-1}-\frac{2\alpha}{h_{x}^{2}}u_{% i}+\frac{\alpha}{h_{x}^{2}}u_{i+1}+f(t,x_{i}),

(6.23)

para $i=2,3,\dotsc,n_{x}$ . Agora, das condições de contorno (6.22c), temos $u_{1}=0$ e $u_{n}=0$ . Com isso, obtemos o seguinte sistema de equações diferenciais ordinárias


	$\displaystyle\frac{du_{2}}{dt}=-\frac{2\alpha}{h_{x}^{2}}u_{2}+\frac{\alpha}{h% _{x}^{2}}u_{3}+f(t,x_{2}),$		(6.24a)
	$\displaystyle\frac{du_{i}}{dt}=\frac{\alpha}{h_{x}^{2}}u_{i-1}-\frac{2\alpha}{% h_{x}^{2}}u_{i}+\frac{\alpha}{h_{x}^{2}}u_{i+1}+f(t,x_{i}),$		(6.24b)
	$\displaystyle\frac{du_{n}}{dt}=\frac{\alpha}{h_{x}^{2}}u_{n-2}-\frac{2\alpha}{% h_{x}^{2}}u_{n-1}+f(t,x_{n-1}),$		(6.24c)

onde $i=3,4,\dotsc,n-1$ e com condições iniciais dadas por (6.22b), i.e.

u_{i}(0)=u_{0}(x_{i}),

(6.25)

para $i=2,3,\dotsc,n$ . Este sistema pode ser escrito na seguinte forma matricial

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\frac{d\tilde% {\boldsymbol{u}}}{dt}=A\tilde{\boldsymbol{u}}+\tilde{f},

(6.26)

onde $\tilde{\boldsymbol{u}}(t)=(u_{2}(t),u_{3}(t),\dotsc,u_{n}(t))$ , $\tilde{f}(t)=(f(t,x_{2}),f(t,x_{3}),\dotsc,f(t,x_{n}))$ e $A$ é uma matriz $(n-1)\times(n-1)$ da forma

A=\begin{bmatrix}-\frac{2\alpha}{h^{2}}&\frac{\alpha}{h^{2}}&0&0&0&\cdots&0&0% \\ \frac{\alpha}{h^{2}}&-\frac{2\alpha}{h^{2}}&\frac{\alpha}{h^{2}}&0&0&\cdots&0&% 0\\ 0&\frac{\alpha}{h^{2}}&-\frac{2\alpha}{h^{2}}&\frac{\alpha}{h^{2}}&0&\cdots&0&% 0\\ 0&0&\ddots&\ddots&\ddots&\cdots&0&0\\ 0&0&&0&0&\cdots&\frac{\alpha}{h^{2}}&-\frac{2\alpha}{h^{2}}\end{bmatrix}.

(6.27)

2. Discretização Temporal.

Para a discretização temporal vamos usar o esquema- $\theta$ . Consideramos os tempos discretos $t^{(k)}=kh_{t}$ , com passo no tempo $h_{t}=t_{f}/n_{t}$ , para $k=0,1,2,\dotsc,n_{t}$ . Denotando $u^{(k)}_{i}\approx u\left(t^{(k)},x_{i}\right)$ , o esquema consiste nas iterações


	$\displaystyle\tilde{\boldsymbol{u}}^{(0)}=\tilde{\boldsymbol{u}}_{0}$		(6.28a)
	$\displaystyle\tilde{\boldsymbol{u}}^{(k+1)}=\tilde{\boldsymbol{u}}^{(k)}+(1-% \theta)h_{t}\left(A\tilde{\boldsymbol{u}}^{(k)}+\tilde{\boldsymbol{f}}^{(k)}\right)$
	$\displaystyle\qquad\qquad\qquad+\theta h_{t}\left(A\tilde{\boldsymbol{u}}^{(k+% 1)}+\tilde{\boldsymbol{f}}^{(k+1)}\right),$		(6.28b)

para $k=0,1,\dotsc,n_{t}-1$ e para um escolhido $0\leq\theta\leq 1$ . No caso, $f$ não depende de $u$ e a Eq. (6.28b) é equivalente ao sistema linear

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\left(I-% \theta h_{t}A\right)\tilde{\boldsymbol{u}}^{(k+1)}=\left[I+(1-\theta)h_{t}A% \right]\tilde{\boldsymbol{u}}^{(k)}+h_{t}\tilde{\boldsymbol{f}}_{\theta}^{(k)},

(6.29)

com $\tilde{\boldsymbol{f}}_{\theta}^{(k)}=(1-\theta)\tilde{\boldsymbol{f}}^{(k)}+% \theta\tilde{\boldsymbol{f}}^{(k+1)}$ .

Observação 6.2.1.

(Estabilidade e Erro de Truncamento.) Para $\theta=0$ (Método de Euler Explícito) o esquema numérico condicionalmente estável [2, Cap. 12, Seç. 2] para

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\alpha\frac{h% _{t}}{h^{2}}\leq\frac{1}{2}.

(6.30)

Para $\theta=1$ (Método de Euler Implícito) o esquema é incondicionalmente estável. Em ambos estes casos, o erro de truncamento é $O(h_{t}+h_{x}^{2})$ . Escolhendo-se $\theta=\frac{1}{2}$ (Método de Crank-Nicolson), o esquema numérico é incondicionalmente estável e com erro de truncamento $O(h_{t}^{2}+h_{x}^{2})$ .

Exemplo 6.2.1.

Consideramos o seguinte problema de calor


	$\displaystyle u_{t}-u_{xx}=(\pi^{2}-1)e^{-t}\operatorname{sen}(\pi x),% \leavevmode\nobreak\ 0<t\leq 1,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.31a)
	$\displaystyle u(0,x)=\operatorname{sen}(\pi x),\leavevmode\nobreak\ 0\leq x% \leq 1,$		(6.31b)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.31c)

Este problema tem solução exata $u(t,x)=e^{-t}\operatorname{sen}(\pi x)$ . A Figura 6.4 mostra o gráfico de superfície $u=u(t,x)$ da solução numérica. Na Figura 6.5, temos a comparação entre a solução numérica e a solução exata (isolinhas).

Refer to caption — Figura 6.4: Solução aproximada do problema de calor do Exemplo 6.2.1.

Código 21: calor.py

⬇

1import numpy as np

2import numpy.linalg as npla

4# params

5alpha = 1.

6theta = 0.5

8# malha temporal

9nt = 10

10ht = 1./nt

11tt = np.linspace(0., 1., nt+1)

13# malha espacial

14nx = 10

15h = 1./n

16xx = np.linspace(0., 1., n+1)

18# rhs

19def f(t, x):

20 return (np.pi**2-1)*np.exp(-t)*np.sin(np.pi*x)

22# axiliares

23lbda = alpha/h**2

25# matriz difusão

26A = np.zeros(((nx-1), (nx-1)))

27A[0,0] = -2*lbda

28A[0,1] = lbda

29for i in range(1,nx-2):

30 A[i,i-1] = lbda

31 A[i,i] = -2*lbda

32 A[i,i+1] = lbda

33A[nx-2,nx-3] = lbda

34A[nx-2,nx-2] = -2*lbda

36# matrizes auxiliares

37Jth = np.identity(A.shape[0]) - theta*ht*A

38J1th = np.identity(A.shape[0]) + (1-theta)*ht*A

40# c.i.

41u0 = np.sin(np.pi * xx)

43# laço no tempo

44u = u0.copy()

45for k in range(nt):

46 print(f'{k+1}: t = {tt[k+1]:f}')

47 fth = (1-theta)*f(tt[k],xx[1:-1]) + theta*f(tt[k+1],xx[1:-1])

48 u[1:-1] = npla.solve(Jth, J1th@u0[1:-1]+ht*fth)

49 u0 = u.copy()

6.2.1 Exercícios

E. 6.2.1.

Considere o problema


	$\displaystyle u_{t}-u_{xx}=0,\leavevmode\nobreak\ 0<t\leq 1,\leavevmode% \nobreak\ -\pi<x<\pi,$		(6.32a)
	$\displaystyle u(0,x)=\operatorname{sen}(x),\leavevmode\nobreak\ -\pi\leq x\leq\pi,$		(6.32b)
	$\displaystyle u(t,-\pi)=u(t,\pi)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.32c)

Sua solução exata é $u(t,x)=e^{-t}\operatorname{sen}(x)$ . Implemente o MDF com esquema- $\theta$ em uma malha uniforme de tamanho espacial $h_{x}$ e passo no tempo $h_{t}$ para obter uma solução numérica $\boldsymbol{u}_{h_{x},h_{t}}$ . Então, verifique a taxa de convergência do erro $\varepsilon_{h_{x},h_{t}}:=\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{2}$ para os diferentes esquemas:

a)

Euler Explícito: $\theta=0$ .
b)

Euler Implícito: $\theta=1$ .
c)

Crank-Nicolson: $\theta=\frac{1}{2}$ .

E. 6.2.2.

Considere o problema


	$\displaystyle u_{t}-\alpha u_{xx}=0,\leavevmode\nobreak\ 0<t\leq 1,\leavevmode% \nobreak\ -\pi<x<\pi,$		(6.33a)
	$\displaystyle u(0,x)=\operatorname{sen}(\alpha x),\leavevmode\nobreak\ -\pi% \leq x\leq\pi,$		(6.33b)
	$\displaystyle u(t,-\pi)=u(t,\pi)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.33c)

Sua solução exata é $u(t,x)=e^{-\alpha^{2}t}\operatorname{sen}(\alpha x)$ . Implemente o MDF com esquema- $\theta$ em uma malha uniforme. Faça testes numéricos para analisar a validade da condição de estabilidade (6.30) para os seguintes esquemas:

a)

Euler Explícito: $\theta=0$ .
b)

Euler Implícito: $\theta=1$ .
c)

Crank-Nicolson: $\theta=\frac{1}{2}$ .

E. 6.2.3.

Considere o problema


	$\displaystyle u_{t}-u_{xx}=\left(\frac{\pi^{2}}{4}-1\right)e^{-t}\cos\left(% \frac{\pi}{2}x\right),\leavevmode\nobreak\ 0<t\leq 1,\leavevmode\nobreak\ 0<x<1,$		(6.34a)
	$\displaystyle u(0,x)=\cos\left(\frac{\pi}{2}x\right),\leavevmode\nobreak\ 0% \leq x\leq 1,$		(6.34b)
	$\displaystyle u(t,0)=e^{-t},\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.34c)
	$\displaystyle u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.34d)

Sua solução exata é $u(t,x)=e^{-t}\cos(\pi x/2)$ . Implemente o MDF com esquema- $\theta$ em uma malha uniforme de tamanho espacial $h_{x}$ e passo no tempo $h_{t}$ para obter uma solução numérica $\boldsymbol{u}_{h_{x},h_{t}}$ . Então, verifique a taxa de convergência do erro $\varepsilon_{h_{x},h_{t}}:=\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{2}$ para os diferentes esquemas:

a)

Euler Explícito: $\theta=0$ .
b)

Euler Implícito: $\theta=1$ .
c)

Crank-Nicolson: $\theta=\frac{1}{2}$ .

E. 6.2.4.

Considere o problema


	$\displaystyle u_{t}-u_{xx}=\left(\frac{\pi^{2}}{4}-1\right)e^{-t}\cos\left(% \frac{\pi}{2}x\right),\leavevmode\nobreak\ 0<t\leq 1,\leavevmode\nobreak\ 0<x<1,$		(6.35a)
	$\displaystyle u(0,x)=\cos\left(\frac{\pi}{2}x\right),\leavevmode\nobreak\ 0% \leq x\leq 1,$		(6.35b)
	$\displaystyle u_{x}(t,0)=0,\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.35c)
	$\displaystyle u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.35d)

Sua solução exata é $u(t,x)=e^{-t}\cos(\pi x/2)$ . Implemente o MDF com o Método de Crank-Nicolson em uma malha uniforme para obter uma solução numérica $\boldsymbol{u}_{h_{x},h_{t}}$ . Então, verifique a taxa de convergência do erro $\varepsilon_{h_{x},h_{t}}:=\|\boldsymbol{u}_{h}-\boldsymbol{u}\|_{2}$ para os seguintes diferentes esquemas:

a)

empregando a diferença finita $D_{+,h_{x}}$ na condição de contorno de Neumann.
b)

empregando a diferença finita $D_{+,h_{x}^{2}}$ na condição de contorno de Neumann.

E. 6.2.5.

Considere o seguinte problema de calor


	$\displaystyle u_{t}-u_{xx}=(\pi^{2}-1)e^{-t}\operatorname{sen}(\pi x),% \leavevmode\nobreak\ 0<t\leq 1,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.36a)
	$\displaystyle u(0,x)=\operatorname{sen}(\pi x),\leavevmode\nobreak\ 0\leq x% \leq 1,$		(6.36b)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.36c)

Sua solução exata $u(t,x)=e^{-t}\operatorname{sen}(\pi x)$ . Faça implementações numéricas do Método das Linhas com MDF na discretização espacial e empregando os seguintes métodos de Runge-Kutta para resolver o sistema de EDOs associado:

a)

Método do Ponto Médio.
b)

Método de R-K-4.

E. 6.2.6.

(Equação de Burgers.) Considere o problema


	$\displaystyle u_{t}+uu_{x}=\alpha u_{xx},\leavevmode\nobreak\ 0<t\leq 1,% \leavevmode\nobreak\ 0<x<1,$		(6.37a)
	$\displaystyle u(0,x)=2\alpha\pi\frac{\operatorname{sen}(\pi x)}{2+\cos(\pi x)}% ,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.37b)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.37c)

Sua solução analítica é [9]

u(t,x)=2\alpha\pi\frac{e^{-\alpha\pi^{2}t}\operatorname{sen}(\pi x)}{2+e^{-% \alpha\pi^{2}t}\cos(\pi x)}.

(6.38)

Faça uma implementação numérica com MDF e com esquema- $\theta$ para resolver este problema. Teste os esquemas para $\alpha=1.,0.1,0.01,0.001$ .

Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Política de uso de dados

Política de uso de dados

Matemática Numérica II

6.2 Equação do Calor

Observação 6.2.1.

Exemplo 6.2.1.

6.2.1 Exercícios

E. 6.2.1.

E. 6.2.2.

E. 6.2.3.

E. 6.2.4.

E. 6.2.5.

E. 6.2.6.

Envie seu comentário