Matemática numérica II

6 Equações Diferenciais Parciais 6.2 Equação do calor Referências

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

6.3 Equação da onda

Consideramos a equação da onda com condições iniciais dadas e condições de contorno de Dirichlet homogêneas

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u_{tt}-\alpha u_{xx}=0,\leavevmode\nobreak\ 0<t\leq t_{f},\leavevmode\nobreak% \ a<x<b,$		(6.81)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u(0,x)=f(x),\leavevmode\nobreak\ a\leq x\leq b,$		(6.82)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u_{t}(0,x)=g(x),\leavevmode\nobreak\ a\leq x\leq b,$		(6.83)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u(t,a)=u(t,b)=0,\leavevmode\nobreak\ 0\leq t\leq t_{f},$		(6.84)

onde $u=u(t,x)$ é a incógnita com $f$ , $g$ e $\alpha>0$ dadas.

Para a aplicação do Método das Diferenças Finitas (MDF), assumimos as discretizações: no tempo, $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}t^{(k)}=kh_{t}$ , $j=0,1,\dotsc,n_{t}$ , $h_{t}=t_{f}/n_{t}$ ; no espaço $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}x_{i}=a+(i-1)% h_{x}$ , $i=1,2,\dotsc,n_{x}+1$ , $h_{x}=(b-a)/n_{x}$ . Então, assumindo a notação $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u^{(k)}_{i}% \approx u\left(t^{(k)},x_{i}\right)$ usando a fórmula de diferenças finitas central $D^{2}_{0,h^{2}}$ , obtemos a seguinte forma discreta da equação Eq. (6.81)

	$\displaystyle\frac{u^{(k-1)}_{i}-2u^{(k)}_{i}+\color[rgb]{0,0,1}\definecolor[% named]{pgfstrokecolor}{rgb}{0,0,1}u^{(k+1)}_{i}}{h_{t}^{2}}$		(6.85)
	$\displaystyle\quad-\alpha\frac{u^{(k)}_{i}-2u^{(k)}_{i}+u^{(k)}_{i+1}}{h_{x}^{% 2}}=0,$		(6.86)

para $i=2,3,\dotsc,n_{x}$ , $j=1,2,\dotsc,n_{t}-1$ . Denotando $\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\lambda:=% \alpha h_{t}^{2}/h_{x}^{2}$ , rearranjando os termos e aplicando as condições de contorno, obtemos

	$\displaystyle u^{(k+1)}_{2}=2\left(1-\lambda\right)u^{(k)}_{2}+\lambda u^{(k)}% _{3}-u^{(k-1)}_{2},$		(6.87)
	$\displaystyle u^{(k+1)}_{i}=\lambda u^{(k)}_{i}+2(1-\lambda)u^{(k)}_{i}+% \lambda u^{(k)}_{i+1}-u^{(k-1)}_{i},$		(6.88)
	$\displaystyle u^{(k+1)}_{n_{x}}=\lambda u^{(k)}_{n_{x}-1}+2(1-\lambda)u^{(k)}_% {n_{x}}-u^{(k-1)}_{n_{x}},$		(6.89)

para $i=2,3,\dotsc,n_{x}$ , $j=1,2,\dotsc,n_{t}-1$ . Ou, equivalentemente, na forma matricial

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\tilde{\bm{u}% }^{(k+1)}=A\tilde{\bm{u}}^{(k)}-\tilde{\bm{u}}^{(k-1)},

(6.90)

para $k=1,2,\dotsc,n_{t}-1$ , onde $\tilde{\bm{u}}^{(k)}=\left(u^{(k)}_{i}\right)_{i=2}^{n_{x}}$ e $A=[a_{i,j}]_{i,j=1}^{n_{x}-1,n_{x}-1}$ é a matriz tridiagonal de elementos

a_{i,j}=\left\{\begin{array}[]{ll}a_{i,i-1}=\lambda&,1<i\leq n_{x}-1,\\ a_{i,i}=2(1-\lambda)&,1\leq i\leq n_{x}-1,\\ a_{i,i+1}=\lambda&,1\leq i<n_{x}-1.\\ \end{array}\right.

(6.91)

Para a inicialização, a Eq. (6.90) requer que conhecemos $\tilde{\bm{u}}^{(0)}$ e $\tilde{\bm{u}}^{(1)}$ . A primeira, vem diretamente da condição inicial Eq. (6.82), i.e.

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\tilde{\bm{u}% }^{(0)}=f\left(\tilde{\bm{x}}\right),

(6.92)

onde $\tilde{\bm{x}}=\left(x_{i}\right)_{i=2}^{n_{x}}$ . Agora, aplicando a fórmula de diferenças finitas progressiva $D_{+,h}$ , temos da condição inicial Eq. (6.83)

\frac{\tilde{\bm{u}}^{(1)}-\tilde{\bm{u}}^{(0)}}{h_{t}}=g\left(\tilde{\bm{x}}\right)

(6.93)

ou, equivalentemente,

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\tilde{\bm{u}% }^{(1)}=\tilde{\bm{u}}^{(0)}+h_{t}g\left(\tilde{\bm{x}}\right).

(6.94)

De tudo isso, temos que a solução numérica da equação da onda pode ser computada com a seguinte iteração

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \tilde{\bm{u}}^{(0)}=f\left(\tilde{\bm{x}}\right),$		(6.95)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \tilde{\bm{u}}^{(1)}=\tilde{\bm{u}}^{(0)}+h_{t}g\left(\tilde{\bm{x}}\right),$		(6.96)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \tilde{\bm{u}}^{(k+1)}=A\tilde{\bm{u}}^{(k)}-\tilde{\bm{u}}^{(k-1)},$		(6.97)

para $k=1,2,\dotsc,n_{t}-1$ , com $\bm{u}^{(k)}=\left(0,\tilde{\bm{u}},0\right)$ .

Observação 6.3.1.(Estabilidade e erro de truncamento)

Pode-se mostrar a seguinte condição de estabilidade [3, p. 487]

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

(6.98)

Com isso e para $f$ e $g$ suficientemente suaves, o esquema numérica (6.95)-(6.97) tem erro de truncamento $O(h_{t}^{2}+h_{x}^{2})$ .

Exemplo 6.3.1.

Consideramos o seguinte problema

	$\displaystyle u_{tt}-u_{xx}=0,\leavevmode\nobreak\ 0<t\leq 2,\leavevmode% \nobreak\ 0<x<1,$		(6.99)
	$\displaystyle u(0,x)=0,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.100)
	$\displaystyle u_{t}(0,x)=\pi\operatorname{sen}(\pi x),\leavevmode\nobreak\ 0% \leq x\leq 1,$		(6.101)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 2.$		(6.102)

Sua solução exata é $u(t,x)=\operatorname{sen}(\pi t)\operatorname{sen}(\pi x)$ . A Figura 6.6 contém gráficos de comparação entra as soluções numérica e exata. Para a solução numérica, tomamos $n_{t}=40$ ( $h_{t}=0.05$ ) e $n_{x}=10$ ( $h_{x}=0.1$ ).

Refer to caption — Figura 6.6: Gráficos comparativos das soluções numérica e exata do problema de onda do Exemplo 6.3.1.

⬇

1import numpy as np

2from numpy import pi, sin, cos

4# params

5nt = 40

6ht = 2./nt

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

9nx = 10

10hx = 1./nx

11xx = np.linspace(0., 1., nx+1)

13# c.i.s

14def f(x):

15 return np.zeros_like(x)

17def g(x):

18 return pi*sin(pi*x)

20# axiliares

21lbda = ht**2/hx**2

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

24A[0,0] = 2*(1. - lbda)

25A[0,1] = lbda

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

27 A[i,i-1] = lbda

28 A[i,i] = 2*(1 - lbda)

29 A[i,i+1] = lbda

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

31A[nx-2,nx-2] = 2*(1 - lbda)

33# laço no tempo

34## c.i.s

35u0 = f(xx)

37u1 = u0.copy()

38u1[1:-1] = u0[1:-1] + ht*g(xx[1:-1])

40u = u1.copy()

41for k in range(1,nt):

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

45 u[1:-1] = A@u1[1:-1] - u0[1:-1]

47 u0 = u1.copy()

48 u1 = u.copy()

6.3.1 Exercício

E. 6.3.1.

Considere o problema

	$\displaystyle u_{tt}-u_{xx}=0,\leavevmode\nobreak\ 0<t\leq 1.5,\leavevmode% \nobreak\ 0<x<1,$		(6.103)
	$\displaystyle u(0,x)=0,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.104)
	$\displaystyle u_{t}(0,x)=\pi\operatorname{sen}(\pi x),\leavevmode\nobreak\ 0% \leq x\leq 1,$		(6.105)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.5.$		(6.106)

Sua solução exata é $u(t,x)=\operatorname{sen}(\pi t)\operatorname{sen}(\pi x)$ . Faça testes numéricos para determinar os passos $h_{t}$ e $h_{x}$ para os quais o esquema numérico (6.95)-(6.97) compute o valor de $u(1.5,0.5)$ com $5$ dígitos significativos corretos.

$h_{t}=2.5\mathrm{e}\!-3\!$ , $h_{x}=1.\mathrm{e}\!-2\!$

E. 6.3.2.

Considere o problema

	$\displaystyle u_{tt}-u_{xx}=e^{-t}(2+x-x^{2}),\leavevmode\nobreak\ 0<t\leq 1,% \leavevmode\nobreak\ 0<x<1,$		(6.107)
	$\displaystyle u(0,x)=x-x^{2},\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.108)
	$\displaystyle u_{t}(0,x)=x^{2}-x,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.109)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.110)

Sua solução exata é $u(t,x)=e^{-t}(x-x^{2})$ . Implemente um esquema numérico semelhante ao (6.95)-(6.97) para computar soluções numéricas desse problema.

E. 6.3.3.

Considere o problema

	$\displaystyle u_{tt}-u_{xx}=e^{-t}(2+x-x^{2}),\leavevmode\nobreak\ 0<t\leq 1,% \leavevmode\nobreak\ 0<x<1,$		(6.111)
	$\displaystyle u(0,x)=x-x^{2},\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.112)
	$\displaystyle u_{t}(0,x)=x^{2}-x,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.113)
	$\displaystyle u_{x}(t,0)=e^{-t},\leavevmode\nobreak\ 0\leq t\leq 1,$		(6.114)
	$\displaystyle u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 1.$		(6.115)

Sua solução exata é $u(t,x)=e^{-t}(x-x^{2})$ . Implemente um esquema numérico semelhante ao (6.95)-(6.97) para computar soluções numéricas desse problema.

E. 6.3.4.

Considere o problema

	$\displaystyle u_{tt}-u_{xx}=0,\leavevmode\nobreak\ 0<t\leq 2,\leavevmode% \nobreak\ 0<x<1,$		(6.116)
	$\displaystyle u(0,x)=0,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.117)
	$\displaystyle u_{t}(0,x)=\pi\operatorname{sen}(\pi x),\leavevmode\nobreak\ 0% \leq x\leq 1,$		(6.118)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 2.$		(6.119)

Sua solução exata é $u(t,x)=\operatorname{sen}(\pi t)\operatorname{sen}(\pi x)$ . Baseado em (6.95)-(6.97), desenvolva um novo esquema numérico substituindo o passo (6.96) por um esquema numérico de mais alta ordem.

Dica: use, por exemplo, um método de R-K-2.

E. 6.3.5.

Considere o problema

	$\displaystyle u_{tt}-\alpha u_{xx}=0,\leavevmode\nobreak\ 0<t\leq 1,% \leavevmode\nobreak\ 0<x<1,$		(6.120)
	$\displaystyle u(0,x)=x(1-x),\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.121)
	$\displaystyle u_{t}(0,x)=0,\leavevmode\nobreak\ 0\leq x\leq 1,$		(6.122)
	$\displaystyle u(t,0)=u(t,1)=0,\leavevmode\nobreak\ 0\leq t\leq 2.$		(6.123)

Use o esquema numérico (6.95)-(6.97) para fazer testes numéricos para $\alpha=1.,0.5,0.1,0.01$ . É necessário ajustar os parâmetros $h_{t}$ e $h_{x}$ ao variar o parâmetro $\alpha$ ? Justifique sua resposta.

Dica: consulte a Observação 6.3.1.

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.3 Equação da onda

Observação 6.3.1.(Estabilidade e erro de truncamento)

Exemplo 6.3.1.

6.3.1 Exercício

E. 6.3.1.

E. 6.3.2.

E. 6.3.3.

E. 6.3.4.

E. 6.3.5.

Envie seu comentário