Rkmultid

Prévia do material em texto

61 2.8 – O método de Runge-Kutta multidimensional
0
5
10
15
20
25
1 2 3 4 5
Γ
(x
)
x
Γ(x)— intgama1.py
Γ(x)— Gnuplot
Figura 2.12: Cálculo de Γ(x) com o método de Runge-Kutta: saída de intgam1 com-
parada com a função Γ(x) do programa de plotagem (Gnuplot)
2.11 Resolva, usando o método de Euler de ordem 2, e compare com a solução analítica:
d y
d x
+ y = x2 exp(−x), y(0) = 1.
2.12 Na equação
d y
d x
+
y
x
= sen
�
2πx
L
�
, y(0) = 0,
estude a sensibilidade do h, necessário para produzir ε= 0,001, ao valor L.
2.13 Utilizando um método implícito semi-analítico, resolva
d y
d x
+
y
x
=
ex
x
, y(x) = 0.
2.14 Resolva, utilizando Runge-Kutta:
d y
d x
+ y = sen(x), y(0) = 1.
2.8 –O método de Runge-Kutta multidimensional
Vamos, na sequência, generalizar o método de Runge-Kutta para que ele seja capaz
de resolver sistemas de equações diferenciais ordinárias do tipo
d y
d x
= f (x , y).
Note que y e f são vetores, enquanto que x permanece sendo um escalar.
Matemática Aplicada 62
Listagem 2.24: Listas em Python
1 #!/usr/bin/python
2 # -*- coding: iso-8859-1 -*-
3 # ----------------------------------------------------------
4 # tentalista: tenta manipular uma lista como se ela fosse
5 # um vetor matemático.
6 # ----------------------------------------------------------
7 from __future__ import print_function
8 uu = [1,2] # lista com dois elementos
9 vv = [3,4] # outra lista com dois elementos
10 ww = uu + vv # soma para ver no que dá
11 kk = 3*uu # multiplica pelo escalar 3 para ver no
12 # que dá
13 print('ww = ', ww) ;
14 print('kk = ',kk) ;
A base para a solução de sistemas de equações diferenciais ordinárias com o mé-
todo de Runge-Kutta é muito simples: basta reconhecer que as equações também
“funcionam” vetorialmente! De fato, podemos escrever
k1 = f (xn, yn),
k2 = h f (xn+
h
2
, yn+
1
2
k1),
k3 = h f (xn+
h
2
, yn+
1
2
k2),
k4 = h f (xn+ h, yn+ k3),
yn+1 = yn+
1
6
k1+
1
3
k2+
1
3
k3+
1
6
k4.
As questões fundamentais que nós teremos aqui do ponto de vista de implemen-
tação computacional do método de Runge-Kutta multidimensional são como:
• Representar um vetor computacionalmente.
• Multiplicar um vetor por um escalar.
• Somar dois vetores.
A primeira idéia que vem à cabeça é representar um vetor em Python por uma
lista; assim, por exemplo, o vetor (1,−1) seria representado pela lista [1,-1]. É
preciso lidar com o fato de que em geral nós escrevemos para as componentes de um
vetor em matemática:
u = (u1, u2, u3)
(os subscritos começam em 1), enquanto que os elementos de uma lista u em Python
são u[0], u[1], e u[2].
Vamos tentar então utilizar listas, no programa tentalista.py, mostrado na lis-
tagem 2.24.
O resultado de rodar tentalista.py é
ww = [1, 2, 3, 4]
kk = [1, 2, 1, 2, 1, 2]
63 2.8 – O método de Runge-Kutta multidimensional
Listagem 2.25: tentaarray.py: arrays com Numpy
1 #!/usr/bin/python
2 # -*- coding: iso-8859-1 -*-
3 # ----------------------------------------------------------
4 # tentaarray: usa arrays, em vez de listas
5 # ----------------------------------------------------------
6 from __future__ import print_function
7 from numpy import array
8 uu = array([1,2]) # lista com dois elementos
9 vv = array([3,4]) # outra lista com dois elementos
10 ww = uu + vv # soma para ver no que dá
11 kk = 3*uu # multiplica pelo escalar 3 para ver no
12 # que dá
13 print('ww = ', ww) ;
14 print('kk = ',kk) ;
Que não é o que esperávamos! O operador + aplicado a listas não soma com-
ponente a componente, mas sim concatena as duas listas; e o produto por 3 não
multiplica componente a componente, mas sim repete a lista 3 vezes.
Felizmente, existe uma maneira de se obter o comportamento esperado — mas
não com listas. Um módulo extra-Python (que tem que ser instalado separadamente),
denominado numpy, faz a mágica. A principal contribuição de numpy é que ele provê
um novo tipo chamado array, que se comporta como um vetor. Vamos vê-lo em ação
no próximo programa, tentaarray.py mostrado na listagem 2.25, cujo resultado é
ww = [4 6]
kk = [3 6]
Como podemos ver, agora as coisas “funcionam”.
Para usar numpy, é claro, você tem que baixar o pacote, e ler o manual: procure
em http://numpy.scipy.org.
Vamos agora resolver um caso para o qual possuímos solução analítica. Dado o
sistema
du1
d x
= u2,
du2
d x
= u1,
a sua solução é
u1(x) = k1e
−x + k2e
x ,
u2(x) =−k1e−x + k2ex .
O programa que resolve este sistema é o rungek4v.py mostrado na listagem 2.26;
o mais interessante, e um dos pontos fortes de Python, é que a rotina rk4 segue
inalterada: observe que ela é capaz de receber como entrada um array (y), assim
como uma função que agora devolve array’s (ff), e que ela própria agora devolve
um array.
Com um h= 0.1, os erros relativos médios de u1 e u2 são extremamente pequenos:
ε = 0.000004 em ambos os casos. Graficamente, temos a resultado mostrado na
figura 2.13. Note que cosh x ≈ sinh x para x ¦ 2,5.
http://numpy.scipy.org
Matemática Aplicada 64
Listagem 2.26: rungek4v.py: o método de Runge-Kutta multidimensional
1 #!/usr/bin/python
2 # -*- coding: iso-8859-1 -*-
3 # ----------------------------------------------------------
4 # rungek4v: resolve as equações diferenciais
5 #
6 # du1/dx = u2,
7 # du2/dx = u1,
8 #
9 # usando o método de Runge-Kutta de ordem 4
10 # ----------------------------------------------------------
11 from __future__ import print_function
12 from numpy import array
13 h = 0.1 # passo em x
14 x = [0.0] # x inicial
15 y = [array([1.0,0.0])] # y inicial
16 n = int(10/h) # número de passos
17 from math import sinh, cosh
18 def ff(x,y):
19 return array([y[1],y[0]])
20 def rk4(x,y,h,ff):
21 '''
22 rk4 implementa um passo do método de Runge-Kutta de ordem 4
23 '''
24 k1 = h*ff(x,y)
25 k2 = h*ff(x+h/2,y+k1/2)
26 k3 = h*ff(x+h/2,y+k2/2)
27 k4 = h*ff(x+h,y+k3)
28 yn = y + k1/6.0 + k2/3.0 + k3/3.0 + k4/6.0
29 return yn
30 for i in range(0,n): # loop da solução numérica
31 xn = (i+1)*h
32 yn = rk4(x[i],y[i],h,ff)
33 x.append(xn)
34 y.append(yn)
35 erro0 = 0.0 # calcula o erro relativo médio
36 erro1 = 0.0
37 for i in range(1,n+1):
38 yana0 = cosh(x[i])
39 yana1 = sinh(x[i])
40 erro0 += abs( (y[i][0] - yana0)/yana0 )
41 erro1 += abs( (y[i][1] - yana1)/yana1 )
42 erro0 /= n
43 erro1 /= n
44 print ( 'erro relativo médio = ', '%10.6f %10.6f' % (erro0,erro1) )
45 fou = open('rungek4v.out','wt')
46 for i in range(0,n+1): # imprime o arquivo de saída
47 fou.write( '%12.6f %12.6f %12.6f\n' % (x[i],y[i][0],y[i][1]) )
48 fou.close()
65 2.8 – O método de Runge-Kutta multidimensional
0.0
5.0
10.0
15.0
20.0
25.0
30.0
0 0.5 1 1.5 2 2.5 3 3.5 4
u
(x
)
x
delta x = 0.5
u1(x) (num)
u1(x) (ana)
u2(x) (num)
u2(x) (ana)
Figura 2.13: Solução numérica pelo Método de Runge-Kutta de um sistema de 2
equações diferenciais ordinárias
h(x , t)
x
z
Figura 2.14: Drenagem de um maciço poroso semi-infinito inicialmente totalmente
saturado.