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\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u_{t}-\alpha u_{xx}=f(t,x),\leavevmode\nobreak\ 0<t\leq t_{f},\leavevmode% \nobreak\ a<x<b,$		(6.43)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u(0,x)=u_{0}(x),\leavevmode\nobreak\ a<x<b,$		(6.44)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u(t,a)=u(t,b)=0,\leavevmode\nobreak\ 0<t\leq t_{f},$		(6.45)

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

O problema (6.43)-(6.45) é 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.43) 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.46)

para $i=2,3,\dotsc,n_{x}$ . Agora, das condições de contorno (6.45), 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.47)
	$\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.48)
	$\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.49)

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

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

(6.50)

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% {\bm{u}}}{dt}=A\tilde{\bm{u}}+\tilde{f},

(6.51)

onde $\tilde{\bm{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.52)

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\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \tilde{\bm{u}}^{(0)}=\tilde{\bm{u}}_{0}$		(6.53)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \tilde{\bm{u}}^{(k+1)}=\tilde{\bm{u}}^{(k)}+(1-\theta)h_{t}\left(A\tilde{\bm{u% }}^{(k)}+\tilde{\bm{f}}^{(k)}\right)$
	$\displaystyle\text{}\qquad+\theta h_{t}\left(A\tilde{\bm{u}}^{(k+1)}+\tilde{% \bm{f}}^{(k+1)}\right),$		(6.54)

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.54) é equivalente ao sistema linear

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

(6.55)

com $\tilde{\bm{f}}_{\theta}^{(k)}=(1-\theta)\tilde{\bm{f}}^{(k)}+\theta\tilde{\bm{% 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.56)

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.57)
	$\displaystyle u(0,x)=\operatorname{sen}(\pi x),\leavevmode\nobreak\ 0\leq x% \leq 1,$		(6.58)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.59)

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# auxiliares

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.60)
	$\displaystyle u(0,x)=\operatorname{sen}(x),\leavevmode\nobreak\ -\pi\leq x\leq\pi,$		(6.61)
	$\displaystyle u(t,-\pi)=u(t,\pi)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.62)

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 $\bm{u}_{h_{x},h_{t}}$ . Então, verifique a taxa de convergência do erro $\varepsilon_{h_{x},h_{t}}:=\|\bm{u}_{h}-\bm{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.63)
	$\displaystyle u(0,x)=\operatorname{sen}(\alpha x),\leavevmode\nobreak\ -\pi% \leq x\leq\pi,$		(6.64)
	$\displaystyle u(t,-\pi)=u(t,\pi)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.65)

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.56) 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.66)
	$\displaystyle u(0,x)=\cos\left(\frac{\pi}{2}x\right),\leavevmode\nobreak\ 0% \leq x\leq 1,$		(6.67)
	$\displaystyle u(t,0)=e^{-t},\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.68)
	$\displaystyle u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.69)

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 $\bm{u}_{h_{x},h_{t}}$ . Então, verifique a taxa de convergência do erro $\varepsilon_{h_{x},h_{t}}:=\|\bm{u}_{h}-\bm{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.70)
	$\displaystyle u(0,x)=\cos\left(\frac{\pi}{2}x\right),\leavevmode\nobreak\ 0% \leq x\leq 1,$		(6.71)
	$\displaystyle u_{x}(t,0)=0,\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.72)
	$\displaystyle u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.73)

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 $\bm{u}_{h_{x},h_{t}}$ . Então, verifique a taxa de convergência do erro $\varepsilon_{h_{x},h_{t}}:=\|\bm{u}_{h}-\bm{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.74)
	$\displaystyle u(0,x)=\operatorname{sen}(\pi x),\leavevmode\nobreak\ 0\leq x% \leq 1,$		(6.75)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.76)

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.77)
	$\displaystyle u(0,x)=2\alpha\pi\frac{\operatorname{sen}(\pi x)}{2+\cos(\pi x)}% ,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.78)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.79)

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.80)

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.(Estabilidade e erro de truncamento)

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.(Equação de Burgers)

Envie seu comentário