Matemática Numérica II

6 Equações Diferenciais Parciais 6 Equações Diferenciais Parciais 6.2 Equação do Calor

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

6.1 Equação de Poisson

Consideramos a equação de Poisson¹¹1Siméon Denis Poisson, 1781 - 1840, matemático francês. Fonte: Wikipédia: Siméon Denis Poisson. (ou equação de Laplace²²2Pierre-Simon Laplace, 1749 - 1827, matemático francês. Fonte: Wikipédia: Pierre-Simon Laplace. heterogênea) no domínio retangular $D=(a,b)\times(c,d)$ com condições de contorno de Dirichlet homogêneas

	$\displaystyle\Delta u=f(x,y),\leavevmode\nobreak\ (x,y)\in D,$		(6.1)
	$\displaystyle u(x,y)=0,\leavevmode\nobreak\ \partial D,$		(6.2)

onde $u=u(x,y)$ é a incógnita, $\Delta u:=u_{xx}+u_{yy}$ e $\partial D$ é a fronteira do domínio $D$ .

A aplicação do Método de Diferenças Finitas para resolver este problema consiste dos mesmos passos usados para resolver problemas de valores de contorno (consulte Seção 5.1), a saber: 1. discretização do domínio, 2. discretização das equações, 3. resolução do problema discreto.

1. Discretização do Domínio (Malha).

Tratando-se do domínio retangular $\overline{D}=[a,b]\times[c,d]$ , podemos construir uma malha do produto cartesiano de partições uniformes dos intervalos $[a,b]$ e $[c,d]$ . Ou seja, tomamos

	$\displaystyle x_{i}:=a+(i-1)h_{x},$		(6.3)
	$\displaystyle y_{j}:=c+(j-1)h_{y},$		(6.4)

com $i=1,2,\dotsc,n_{x}+1$ , $j=1,2,\dotsc,n_{y}+1$ , sendo $n_{x}$ e $n_{y}$ o número de subintervalos escolhidos para as partições, respectivamente, e os passos $h_{x}=(b-a)/n_{x}$ e $h_{y}=(d-c)/n_{y}$ . O tamanho da malha é definido por $h:=\max\{h_{x},h_{y}\}$ .

Refer to caption — Figura 6.1: Malha bidimensional.

O produto cartesiano das partições em $x$ e $y$ nos fornece uma partição do domínio $\overline{D}$ da forma

P(\overline{D})=\{(x_{1},y_{1}),(x_{1},y_{2}),\dotsc,(x_{i},y_{j}),\dotsc,(x_{% n_{x}},y_{n_{y}})\},

(6.5)

cujos nodos $(x_{i},y_{j})$ podem ser enumerados (indexados) por $k=i+(j-1)(n_{x}+1)$ . Por simplicidade, no decorrer do texto, assumiremos $n_{x}=n_{y}=:n$ e, por conseguinte, $h_{x}=h_{y}=h$ e temos a enumeração

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}k=i+(j-1)(n+1).

(6.6)

Consulte a Figura 6.1.

2. Discretização das Equações.

Usando a fórmula de diferenças finitas central de ordem $h^{2}$ para a segunda derivada, temos

	$\displaystyle u_{xx}(x,y)=\frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^{2}}+O(h^{2}),$		(6.7)
	$\displaystyle u_{yy}(x,y)=\frac{u(x,y+h)-2u(x,y)+u(x,y-h)}{h^{2}}+O(h^{2}).$		(6.8)

Daí, denotando $u_{ij}\approx u(x_{i},y_{j})$ temos

	$\displaystyle u_{xx}(x_{i},y_{j})=\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}+O% (h^{2}),$		(6.9)
	$\displaystyle u_{yy}(x_{i},y_{j})=\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^{2}}+O% (h^{2}).$		(6.10)

Então, da (6.1) temos

	$\displaystyle\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^{2}}$
	$\displaystyle\text{}\qquad+\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^{2}}+O(h^{2})% =f(x_{i},y_{j}).$		(6.11)

Agora, com base na enumeração (6.6) denotamos $u_{k}:=u_{i+(j-1)(n+1)}$ , desprezando o erro de truncamento e rearranjando os termos, obtemos

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\frac{1}{h^{2% }}u_{k-n}+\frac{1}{h^{2}}u_{k-1}-\frac{4}{h^{2}}u_{k}+\frac{1}{h^{2}}u_{k+1}+% \frac{1}{h^{2}}u_{k+n}=f_{k},

(6.12)

para $k=i+(j+1)(n+1)$ com $i,j=2,3,\dotsc,n$ (nodos internos). Isto é, esta última expressão nos fornece um sistema de $(n-1)^{2}$ equações para $(n+1)^{2}$ incógnitas $\boldsymbol{u}=(u_{k})_{k=1}^{(n+1)^{2}}$ . Para fechar o sistema, usamos as condições de contorno (6.2)

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}u_{k}=0

(6.13)

para $k=i+(j+1)(n+1)$ com $i=1,n+1$ e $j=1,2,\dotsc,n+1$ , ou $i=2,3,\dotsc,n$ e $j=1,n+1$ .

Com isso, o problema discreto obtido da aplicação do MDF consiste no sistema linear de $(n+1)^{2}\times(n+1)^{2}$ (6.12)-(6.13).

3. Resolução do Problema Discreto.

O problema discreto (6.12)-(6.13) pode ser escrito na forma matricial

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}A\boldsymbol{% u}=\boldsymbol{b},

(6.14)

onde o vetor da incógnitas é $\boldsymbol{u}=(u_{k})_{k=1}^{(n+1)^{2}}$ . A matriz dos coeficientes $A=[a_{l,m}]_{l,m=1}^{(n+1)^{2},(n+1)^{2}}$ e o vetor dos termos contantes $\boldsymbol{b}=(b_{k})_{k=1}^{(n+1)^{2}}$ têm elementos não nulos

	$\displaystyle i=1,n+1,\leavevmode\nobreak\ j=1,2,\dotsc,n+1:$
	$\displaystyle\text{}\qquad a_{k,k}=1,$		(6.15)
	$\displaystyle\text{}\qquad b_{k}=0,$		(6.16)

	$\displaystyle i=1,2,\dotsc,n+1,\leavevmode\nobreak\ j=1,n+1:$
	$\displaystyle\text{}\qquad a_{k,k}=1,$		(6.18)
	$\displaystyle\text{}\qquad b_{k}=0,$		(6.19)

$\displaystyle i,j=2,3,\dotsc,n$	$\displaystyle:$
$\displaystyle\text{}\qquad a_{k,k-n}=\frac{1}{h^{2}},$		(6.21)
$\displaystyle\text{}\qquad a_{k,k-1}=\frac{1}{h^{2}},$		(6.22)
$\displaystyle\text{}\qquad a+{k,k}=-\frac{4}{h^{2}},$		(6.23)
$\displaystyle\text{}\qquad a_{k,k+1}=\frac{1}{h^{2}},$		(6.24)
$\displaystyle\text{}\qquad a_{k,k+n}=\frac{1}{h^{2}},$		(6.25)
$\displaystyle\text{}\qquad b_{k}=f(x_{i},y_{j}).$		(6.26)

Assim sendo, basta empregarmos um método apropriado para resolver o sistema linear (6.14) para obter a solução aproximada de $u$ nos nodos $(x_{i},y_{j})$ .

Exemplo 6.1.1.

Consideramos o seguinte problema

	$\displaystyle\Delta u=-2\pi^{2}\operatorname{sen}(\pi x)\operatorname{sen}(\pi y% ),\leavevmode\nobreak\ (x,y)\in(0,1)^{2},$		(6.27)
	$\displaystyle u=0,\leavevmode\nobreak\ (x,y)\in\partial D.$		(6.28)

A solução exata é $u(x,y)=\operatorname{sen}(\pi x)\operatorname{sen}(\pi y)$ .

A Figura 6.2 mostra o gráfico de superfície da solução aproximada obtida pelo MDF com $h=10^{-1}$ . A Figura 6.3 mostra a comparação entre os gráficos de contorno das soluções numérica e exata (linhas brancas).

Código 20: mdf_poisson.py

⬇

1import numpy as np

2import numpy.linalg as npla

4# malha

5n = 10

6h = 1./n

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

8yy = np.linspace(0., 1., n+1)

10# rhs

11def f(x,y):

12 return -2*np.pi**2*np.sin(np.pi*x)*np.sin(np.pi*y)

14# problema discreto

15A = np.zeros(((n+1)**2, (n+1)**2))

16b = np.empty((n+1)**2)

18# c.c.

19for j in range(n+1):

20 # i = 0

21 k = j*(n+1)

22 A[k,k] = 1.

23 b[k] = 0.

24 # i = n

25 k = n + j*(n+1)

26 A[k,k] = 1.

27 b[k] = 0.

29for i in range(1,n):

30 # j = 0

31 k = i

32 A[k,k] = 1.

33 b[k] = 0.

34 # j = n

35 k = i + n*(n+1)

36 A[k,k] = 1.

37 b[k] = 0.

39# pts internos

40for i in range(1,n):

41 for j in range(1,n):

42 k = i + j*(n+1)

43 A[k,k-n-1] = 1./h**2

44 A[k,k-1] = 1./h**2

45 A[k,k] = -4./h**2

46 A[k,k+1] = 1./h**2

47 A[k,k+n+1] = 1./h**2

48 b[k] = f(xx[i],yy[j])

50# resol p.d.

51u = npla.solve(A, b)

6.1.1 Exercícios

E. 6.1.1.

Use o MDF para encontrar uma solução aproximada do seguinte problema de Poisson

	$\displaystyle\Delta u=-2\pi^{2}\operatorname{sen}(\pi x)\operatorname{sen}(\pi y% ),\leavevmode\nobreak\ (x,y)\in D=(-1,1)^{2},$		(6.29)
	$\displaystyle u=0,\leavevmode\nobreak\ (x,y)\in\partial D.$		(6.30)

A solução exata é $u(x,y)=\operatorname{sen}(\pi x)\operatorname{sen}(\pi y)$ . Faça uma comparação gráfica entre as soluções numérica e exata no caso de $h=10^{-1}$ (malha uniforme). Compare o erro $\varepsilon_{h}:=\|\tilde{\boldsymbol{u}}_{h}-\boldsymbol{u}\|_{2}$ para $n=10,20,40,80,160$ (número de subintervalos na malha uniforme). A taxa de convergência é a esperada? Justifique sua resposta.

E. 6.1.2.

Use o MDF para encontrar uma solução aproximada do seguinte problema de Laplace

	$\displaystyle\Delta u=0,\leavevmode\nobreak\ (x,y)\in(0,1)^{2},$		(6.31)
	$\displaystyle u(0,y)=u(1,y)=y^{2}-y,\leavevmode\nobreak\ 0\leq y\leq 1,$		(6.32)
	$\displaystyle u(x,0)=u(x,1)=x-x^{2},\leavevmode\nobreak\ 0\leq x\leq 1.$		(6.33)

A solução exata é $u(x,y)=x-x^{2}+y^{2}-y$ . Faça uma comparação gráfica entre as soluções numérica e exata no caso de $h=10^{-1}$ . Compare o erro $\varepsilon_{h}:=\|\tilde{\boldsymbol{u}}_{h}-\boldsymbol{u}\|_{2}$ para $n=10,20,40,80,160$ (número de subintervalos na malha uniforme).

E. 6.1.3.

Considere o problema

	$\displaystyle\Delta u=-2\pi^{2}\operatorname{sen}(\pi x)\operatorname{sen}(\pi y% ),\leavevmode\nobreak\ (x,y)\in(0,1)^{2},$		(6.34)
	$\displaystyle u(0,y)=0,\leavevmode\nobreak\ 0\leq y\leq 1,$		(6.35)
	$\displaystyle u_{x}(1,y)=0,\leavevmode\nobreak\ 0\leq y\leq 1,$		(6.36)
	$\displaystyle u(x,0)=u(x,1)=0,\leavevmode\nobreak\ 0\leq x\leq 1.$		(6.37)

A solução exata é $u(x,y)=\operatorname{sen}(\pi x)\operatorname{sen}(\pi y)$ . Com uma malha uniforme, obtenha uma solução aproximada com o MDF empregando, na fronteira com condições de Neumann³³3Carl Gottfried Neumann, 1832 - 1925, matemático alemão. Fonte: Wikipédia: Carl Neumann.:

a)

$D_{-,h}$ fórmulas diferença regressiva de ordem $h$ .
b)

$D_{-,h^{2}}$ diferença regressiva de ordem $h^{2}$ .

Compare a taxa de convergência do erro $\varepsilon_{h}:=\|\tilde{\boldsymbol{u}}_{h}-\boldsymbol{u}\|_{2}$ entre essas duas formulações.

E. 6.1.4.

Considere o problema

	$\displaystyle\Delta u=-2\pi^{2}\operatorname{sen}(\pi x)\operatorname{sen}(\pi y% ),\leavevmode\nobreak\ (x,y)\in(0,1)^{2},$		(6.38)
	$\displaystyle u(0,y)=u(1,y)=0,\leavevmode\nobreak\ 0\leq y\leq 1,$		(6.39)
	$\displaystyle u_{y}(x,0)=u_{y}(x,1)=0,\leavevmode\nobreak\ 0\leq x\leq 1.$		(6.40)

A solução exata é $u(x,y)=\operatorname{sen}(\pi x)\operatorname{sen}(\pi y)$ . Com uma malha uniforme, obtenha uma solução aproximada com o MDF empregando, nas fronteiras com condições de Neumann:

a)

fórmulas de diferenças finitas de $O(h)$ .
b)

fórmulas de diferenças finitas de $O(h^{2})$ .

Compare a taxa de convergência do erro $\varepsilon_{h}:=\|\tilde{\boldsymbol{u}}_{h}-\boldsymbol{u}\|_{2}$ entre essas duas formulações.

E. 6.1.5.