Matemática numérica II

5 Problema de valor de contorno 5.3 Método de Volumes Finitos 6 Equações Diferenciais Parciais

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

5.4 Problemas não-lineares

Vamos estudar a resolução de problemas não-lineares de valores de contorno da forma

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u_{xx}=f\left(x,u,u_{x}\right),\leavevmode\nobreak\ a\leq x\leq b,$		(5.139)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u(a)=u_{a},$		(5.140)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u(b)=u_{b},$		(5.141)

onde $f=f(x,u,u_{x})$ é uma função não linear para $u$ ou $u_{x}$ .

Empregando o método de diferenças finitas (MDF), começamos assumindo uma malha uniforme de $n$ -subintervalos com nodos $x_{i}=a+(i-1)h$ , tamanho de malha $h=(b-a)/n$ , $i=1,2,\dotsc,n+1$ . Denotando $u_{i}\approx u(x_{i})$ e aplicando fórmulas de diferenças finitas centrais para $u_{xx}$ e $u_{x}$ , a Eq. (5.139) fornece

	$\displaystyle\frac{1}{h^{2}}u_{i-1}-\frac{2}{h^{2}}u_{i}\frac{1}{h^{2}}u_{i+1}=$
	$\displaystyle\text{}\qquad f\left(x_{i},u_{i},\frac{1}{2h}u_{i+1}-\frac{1}{2h}% u_{i-1}\right),$		(5.142)

para $i=2,3,\dotsc,n$ . As condições de contorno Eqs. (5.140)-(5.141), fornecem as equações de fechamento

	$\displaystyle u_{1}=u_{a},$		(5.143)
	$\displaystyle u_{n+1}=u_{b}.$		(5.144)

Com isso, temos que o problema discreto associado consiste em: encontrar $\bm{u}=(u_{i})_{i=1}^{n+1}$ solução do seguinte sistema de equações não-lineares

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u_{1}-u_{a}=0,$		(5.145)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% -\frac{1}{h^{2}}u_{i-1}+\frac{2}{h^{2}}u_{i}-\frac{1}{h^{2}}u_{i+1}$
	$\displaystyle\text{}\qquad\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor% }{rgb}{0,0,1}+f\left(x_{i},u_{i},\frac{1}{2h}u_{i+1}-\frac{1}{2h}u_{i-1}\right% )=0,$		(5.146)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% u_{n+1}-u_{b}=0.$		(5.147)

A resolução do problema discreto (LABEL:cap_pvc_sec_pnlin:eq:pd) pode ser feito com o Método de Newton⁵⁵5Isaac Newton, 1642 - 1727, matemático, físico, astrônomo, teólogo e autor inglês. Fonte: Wikipédia: Isaac Newton.. Para tando, observamos que o sistema tem a forma vetorial

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}F(\bm{u})=\bm% {0},

(5.148)

onde $F(\bm{u})=\left(f_{i}(\bm{u})\right)_{i=1}^{(n+1)}$ é a função vetorial de componentes

	$\displaystyle f_{1}(\bm{u})=u_{1}-u_{a},$		(5.149)
	$\displaystyle f_{i}(\bm{u})=-\frac{1}{h^{2}}u_{i-1}+\frac{2}{h^{2}}u_{i}-\frac% {1}{h^{2}}u_{i+1}$
	$\displaystyle\text{}\qquad+f\left(x_{i},u_{i},\frac{1}{2h}u_{i+1}-\frac{1}{2h}% u_{i-1}\right),$		(5.150)
	$\displaystyle f_{n+1}(\bm{u})=u_{n+1}-u_{b},$		(5.151)

com $i=2,3,\dotsc,n$ . A iteração do Método de Newton consiste em

	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \bm{u}^{(0)}=\text{aprox. inicial},$		(5.152)
	$\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% \bm{u}^{(k+1)}=\bm{u}^{(k)}+\bm{\delta}^{(k)},$		(5.153)

onde $\bm{\delta}^{(k)}$ é a atualização de Newton computada por

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}J_{F}\left(% \bm{u}^{(k)}\right)\bm{\delta}^{(k)}=-F\left(\bm{u}^{(k)}\right),

(5.154)

para $k=0,1,2,\ldots$ até que um critério de parada seja satisfeito. A matriz jacobiana⁶⁶6Carl Gustav Jakob Jacobi, 1804 - 1851, matemático alemão. Fonte: Wikipédia: Carl Gustav Jakob Jacobi. é denotada por $J_{F}\left(\bm{u}^{(k)}\right)=[\jmath_{i,j}]_{i,j=1}^{n+1,n+1}$ e tem elementos não nulos

\jmath_{1,1}=1,

(5.155)

	$\displaystyle\jmath_{i,i-1}=-\frac{1}{h^{2}}-\frac{1}{2h}f_{u_{x}}\left(x_{i},% u_{i},\frac{1}{2h}u_{i+1}-\frac{1}{2h}u_{i-1}\right),$		(5.156)
	$\displaystyle\jmath_{i,i}=\frac{2}{h^{2}}+f_{u}\left(x_{i},u_{i},\frac{1}{2h}u% _{i+1}-\frac{1}{2h}u_{i-1}\right),$		(5.157)
	$\displaystyle\jmath_{i,i+1}=-\frac{1}{h^{2}}+\frac{1}{2h}f_{u_{x}}\left(x_{i},% u_{i},\frac{1}{2h}u_{i+1}-\frac{1}{2h}u_{i-1}\right),$		(5.158)

para $i=2,3,\dotsc,n$ e

\jmath_{n+1,n+1}=1.

(5.159)

Exemplo 5.4.1.

Vamos considerar o seguinte PVC

	$\displaystyle uu_{x}-u_{xx}=\pi\operatorname{sen}(\pi x)\left[\pi+\cos(\pi x)% \right],\leavevmode\nobreak\ 0<x<1,$		(5.160)
	$\displaystyle u(0)=u(1)=0.$		(5.161)

Rearranjando os termos, podemos escrevê-lo na forma da Eq. (LABEL:cap_pvc_sec_pnlin:eq:pvc), com $u_{a}=u_{b}=0$ e

f(x,u,u_{x})=uu_{x}-\pi\operatorname{sen}(\pi x)\left[\pi+\cos(\pi x)\right].

(5.162)

Com isso, calculamos

	$\displaystyle f_{u}(x,u,u_{x})=u_{x},$		(5.163)
	$\displaystyle f_{u_{x}}(x,u,u_{x})=u.$		(5.164)

Então, a aplicação do MDF-Newton com $h=10^{-1}$ fornece o resultado da Figura 5.4. A solução exata é $u(x)=\sin(\pi x)$ .

Refer to caption — Figura 5.4: Resultado da aplicação do MDF-Newton para o PVC do Exemplo 5.4.1.

Código 19: mdf-newton.py

⬇

1import numpy as np

2import numpy.linalg as npla

3from numpy import pi, sin, cos

5# parâmetros

6n = 10

7h = 1./n

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

10# c.c. Dirichlet

11ua = 0.

12ub = 0.

14def f(x, u, ux):

15 return u*ux - pi*sin(pi*x)*(pi + cos(pi*x))

17def fu(x, u, ux):

18 return ux

20def fux(x, u, ux):

21 return u

23# rhs

24def F(u):

25 y = np.empty(n+1)

26 # f_1

27 y[0] = u[0] - ua

28 # f_i

29 for i in range(1,n):

30 ux = u[i+1]/(2*h) - u[i-1]/(2*h)

31 y[i] = -1./h**2*u[i-1] + 2./h**2*u[i] - 1./h**2*u[i+1] \

32 + f(xx[i], u[i], ux)

33 # f_{n+1}

34 y[n] = u[n] - ub

36 return y

38# jacobiana

39def J(u):

40 J = np.zeros((n+1,n+1))

41 J[0,0] = 1.

42 for i in range(1,n):

43 ux = 1./(2*h)*u[i+1] - 1./(2*h)*u[i-1]

44 J[i,i-1] = -1./h**2 - 1/(2*h)\

45 * fux(xx[i], u[i], ux)

46 J[i,i] = 2/h**2 + fu(xx[i], u[i], ux)

47 J[i,i+1] = -1./h**2 + 1/(2*h)\

48 * fux(xx[i], u[i], ux)

49 J[n,n] = 1.

51 return J

53# aprox inicial

54u = np.zeros(n+1)

56# iterações de Newton

57maxiter = 10

58for k in range(maxiter):

60 # passo de Newton

61 dlta = npla.solve(J(u), -F(u))

63 # atualização

64 u += dlta

66 ndlta = npla.norm(dlta)

67 print(f'{k+1}: norm = {ndlta:.2e}')

68 if (ndlta < 1e-10):

69 print('convergiu.')

70 break

5.4.1 Exercícios

E. 5.4.1.

Considere o PVC

	$\displaystyle u^{2}-u_{xx}=\cos^{2}(\pi x)+\pi^{2}\cos(\pi x),\leavevmode% \nobreak\ 0<x<1,$		(5.165)
	$\displaystyle u(0)=1,$		(5.166)
	$\displaystyle u(1)=-1.$		(5.167)

Este problema tem solução analítica $u(x)=\cos(\pi x)$ . Use o MDF-Newton para computar $u_{h}$ aproximações de $u$ para $h=10^{-1}$ , $10^{-2}$ , $10^{-3}$ , $10^{-4}$ . Então, verifique a convergência com base no erro $\varepsilon_{h}:=\|\tilde{\bm{u}}-\bm{u}\|_{2}$ . A convergência tem a taxa esperada? Justifique sua resposta.

$h$	$\varepsilon_{h}$
$10^{-1}$	$4.0\mathrm{e}\!-\!03$
$10^{-2}$	$1.3\mathrm{e}\!-\!04$
$10^{-3}$	$4.0\mathrm{e}\!-\!06$
$10^{-4}$	$1.3\mathrm{e}\!-\!07$

E. 5.4.2.

Considere o PVC

	$\displaystyle uu_{x}-u_{xx}=2+x(1-x)(1-2x),\leavevmode\nobreak\ 0<x<1,$		(5.168)
	$\displaystyle u_{x}(0)=1,$		(5.169)
	$\displaystyle u(1)=0.$		(5.170)

Este problema tem solução analítica $u(x)=x(1-x)$ . Use o MDF-Newton para computar $u_{h}$ aproximações de $u$ para $h=10^{-1}$ , $10^{-2}$ , $10^{-3}$ :

a)

aplicando as diferenças finitas $D_{0,h^{2}}u(x)$ para $0<x<1$ e $D_{+,h}u(x)$ para $x=0$ .
b)

aplicando as diferenças finitas $D_{0,h^{2}}u(x)$ para $0<x<1$ e $D_{+,h^{2}}u(x)$ para $x=0$ .

Qual dessas formulações tem a melhor taxa de convergência do erro em relação ao passo de malha $h$ ? Justifique e verifique sua resposta.

b) tem melhor taxa de convergência.

E. 5.4.3.

Desenvolva uma versão do método MEF-Newton (Método de Elementos Finitos com o Método de Newton) para computar a solução aproximada do PVC dado no Exemplo 5.4.1. Implemente-o e verifique a convergência do método para $h=10^{-1}$ , $10^{-2}$ e $10^{-3}$ .

E. 5.4.4.

Desenvolva uma versão do método MVF-Newton (Método de Volumes Finitos com o Método de Newton) para computar a solução aproximada do PVC dado no Exemplo 5.4.1. Implemente-o e verifique a convergência do método para $h=10^{-1}$ , $10^{-2}$ e $10^{-3}$ .

E. 5.4.5.

Desenvolva uma versão do método MEF-Newton (Método de Elementos Finitos com o Método de Newton) para computar a solução aproximada do PVC dado no E.5.4.2. Implemente-o e verifique a convergência do método para $h=10^{-1}$ , $10^{-2}$ e $10^{-3}$ .

E. 5.4.6.

Desenvolva uma versão do método MVF-Newton (Método de Volumes Finitos com o Método de Newton) para computar a solução aproximada do PVC dado no E.5.4.2. Implemente-o e verifique a convergência do método para $h=10^{-1}$ , $10^{-2}$ e $10^{-3}$ .

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