Colabore! Saiba mais screen_rotation

Ao navegar por este site , você concorda com a Política de Uso de Dados.

| | | | |

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

6.2 Equação do Calor

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

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.2) é 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}% \pgfsys@color@rgb@stroke{0}{0}{1}\pgfsys@color@rgb@fill{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}% \pgfsys@color@rgb@stroke{0}{0}{1}\pgfsys@color@rgb@fill{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}% \pgfsys@color@rgb@stroke{0}{0}{1}\pgfsys@color@rgb@fill{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

As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a Política de Uso de Dados para mais informações. Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

| | | |