fortran 2

Prévia do material em texto

Operações Matemáticas Básicas Diferenciação Numérica
Diferenciação Numérica
Como avaliar f ′(a)?
Discretização de f : Conhecemos f (xn) = fn num conjunto de pontos xn igualmente
espaçados:
xn = nh (n = 0,±1,±2, · · · ) ; fn = f (xn)
Série de Taylor de f em a: Expansão de f em torno de a
f deve ser infinitamente diferenciável em a.
f (x) = f (a) +
f ′(a)
1!
(x − a) + f
′′(a)
2!
(x − a)2 + f
′′′(a)
3!
(x − a)3 + · · ·
= Σ∞n=0
f (n)(a)
n!
(x − a)n
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 4/84 2015– 4 / 84
Operações Matemáticas Básicas Diferenciação Numérica
Exemplos
Por simplicidade, a = 0 (basta transladar para a arbitrário)
f±1 = f (x = ±h) = f0 ± hf ′0 + h2 f
′′
0
2!
± h3 f
′′′
0
3!
+O(h4)
f±2 = f (x = ±2h) = f0 ± 2hf ′0 + 2h2f ′′0 ± 43h
3
f
′′′
0 +O(h4)
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 5/84 2015– 5 / 84
Operações Matemáticas Básicas Diferenciação Numérica
Aproximações para f ′
Assuma que f e suas derivadas têm mesma ordem de magnitude (hipótese razoável)
Dois pontos (Euller) – (0,±h)
f
′
0 =
f1 − f0
h
+
1
h
O(h2f ′′0 ) (forward)
f
′
0 =
f0 − f−1
h
+
1
h
O(h2f ′′0 ) (backward)
Três pontos (−h,0,h)
f
′
0 =
f+1 − f−1
2h
− h
3
6h
f
′′′
0 +
1
h
O(h5f (5)0 )
| {z }
(nulo para P(x2))
Cinco pontos (−2h,−h,0,h,2h)
f
′
0 =
1
12h
[f−2 − 8(f−1 − f1)− f2] + h
5
30h
f
(5) +
1
h
O(h7)
| {z }
(nulo para P(x4))
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 6/84 2015– 6 / 84
Operações Matemáticas Básicas Diferenciação Numérica
Nota sobre precisão
Euller é exata para P(h)
Três pontos é exata para P(h2)
Cinco pontos é exata para P(h4)
Exemplo: Quanto menor h, mais exato é o cálculo?
Verifique como o erro na aproximação de f ′ varia com h. Sugestão:
Para vários valores de h, compute com os métodos de 2,3 e 5 pontos f ′(1),
onde f = sin(x) (f ′ = cos(x))
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 7/84 2015– 7 / 84
Operações Matemáticas Básicas Diferenciação Numérica
Aproximações para f ′′
Três pontos (−h,· · · ,h)
f
′′ ≈ 1
h2
[f−1 − 2f0 + f1]
Cinco pontos (−2h,· · · ,2h)
f
′′ ≈ 1
12h2
[−f−2 + 16 (f−1 + f1)− 30f0 − f2]
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 8/84 2015– 8 / 84
Operações Matemáticas Básicas Quadratura Numérica
Sumário
3 Operações Matemáticas Básicas
Diferenciação Numérica
Quadratura Numérica
Encontrar, numericamente, f (xi ) = 0
Optimization
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 9/84 2015– 9 / 84
Operações Matemáticas Básicas Quadratura Numérica
Quadratura Numérica
Como avaliar
R
f (x) dx?
Conhecer f no intervalo de avaliação
Exemplo: Para calcular
R b
a
f (x) dx , devemos conhecer f em [a,b].
É conveniente particionar o intervalo num número par de subintervalos, igualmente
espaçados
Ideia básica de avaliação de quadraturas
Aproximar f por uma função integrável, de forma exata, no subintervalo [−h, h]
Estratégia: Com uma fórmula no intervalo [−h, h], decompomos
R b
a
f (x) dx em:
Z b
a
f (x) dx =
Z a+2h
a
f (x) dx +
Z a+4h
a+2h
f (x) dx + · · ·+
Z b
b−2h
f (x) dx
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 10/84 2015– 10 / 84
Operações Matemáticas Básicas Quadratura Numérica
Regra Trapezoidal
Aproxima f por sua série de Taylor O(1) – uma função linear
������������
������������
������������
������������
������������
������������
f−h
f0 f+h
+h−h 0
Figura 1: A função f (linha contínua) é aproximada por retas (linha tracejada) no intervalo [−h, 0] e [0, h]
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 11/84 2015– 11 / 84
Operações Matemáticas Básicas Quadratura Numérica
Intervalo [−h, 0]:
Função: f (x) = f0 +
(f−h−f0)
h−0 x +O(h2)
Z h
0
f (x) dx =
Z h
0
f0 dx +
(f−h − f0)
h − 0
Z h
0
x dx +
Z h
0
O(h2) dx
= f0h +
(f−h − f0)
h
h2
2
+O(h3)
Intervalo [0, h]
Função: f (x) = f0 +
(fh−f0)
h−0 x +O(h2)
Z h
0
f (x) dx =
Z h
0
f0 dx +
(fh − f0)
h − 0
Z h
0
x dx +
Z h
0
O(h2) dx
= f0h +
(fh − f0)
h
h2
2
+O(h3)
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 12/84 2015– 12 / 84
Operações Matemáticas Básicas Quadratura Numérica
Intervalo [−h, h]
Z h
−h
f (x) dx = f0h +
(f−h − f0)
h
h2
2
+O(h3)
| {z }
[−h,0]
+ f0h +
(fh − f0)
h
h2
2
+O(h3)
| {z }
[0,h]
=
h
2
(f−h + 2f0 + fh) +O(h3)
Nota
Erro na aproximação de f é O(h2)
Erro na aproximação de
R
f é O(h3)
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 13/84 2015– 13 / 84
Operações Matemáticas Básicas Quadratura Numérica
Regra de Simpson
Aproxima f por sua série de Taylor O(2) – uma função quadrática
Intervalo [−h, h]. Aproxima um polinômio P(x2) considerando os três pontos do
intervalo:
f (x) = f0 + f
′
0
|{z}
fh−f−h
2h
x +
1
2!
f
′′
|{z}
fh−2f0+f−h
h2
x
2 +
1
3!
f
′′′
0 x
3 +
1
4!
f
(iv)
x
4
Z h
−h
f (x) dx = f0
Z h
−h
dx + f ′0
Z h
−h
x dx
| {z }
=0 (???)
+
1
2!
f
′′
0
Z h
−h
x
2
dx+
1
3!
f
′′′
0
Z h
−h
x
3
dx
| {z }
=0 (???)
+
1
4!
f
(iv)
0
Z h
−h
x
4
dx
=
h
3
(fh + 4f0 + f−h) +O(h5)
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 14/84 2015– 14 / 84
Operações Matemáticas Básicas Quadratura Numérica
Intervalo [−h, h]:
Z h
−h
f (x) dx ==
h
3
(fh + 4f0 + f−h) +O(h5)
Nota
Erro na aproximação de f é O(h2)
Erro na aproximação de
R
f é O(h5)
compare com a Regra do Trapézio
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 15/84 2015– 15 / 84
Operações Matemáticas Básicas Quadratura Numérica
Exemplo
Compute
R 1
0
ex dx :
Analiticamente
Regra do Trapézio
Regra de Simpson
Nota
Quanto maior a ordem do P(x) para interpolar f, melhor a integração?
Não. P(x) de ordem alta oscilam muito rapidamente → a interpolação fica ruim
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 16/84 2015– 16 / 84
Operações Matemáticas Básicas Encontrar, numericamente, f (xi ) = 0
Sumário
3 Operações Matemáticas Básicas
Diferenciação Numérica
Quadratura Numérica
Encontrar, numericamente, f (xi ) = 0
Optimization
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 17/84 2015– 17 / 84
Operações Matemáticas Básicas Encontrar, numericamente, f (xi ) = 0
Método da Busca por f (xi ) = 0
Quando f trocar de sinal . . .
<0 <0
>0 >0
Figura 2: Esquema do método de busca
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 18/84 2015– 18 / 84
Operações Matemáticas Básicas Encontrar, numericamente, f (xi ) = 0
Método de Busca
Sempre converge para alguma raíz se a função for contínua no intervalo
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 19/84 2015– 19 / 84
Operações Matemáticas Básicas Encontrar, numericamente, f (xi ) = 0
Método de Newton–Raphson (precisamos de f ′)
Assume que f é localmente linear próximo a cada xj
i
f(x )i
0
i
1f´(x )
f´(x )i
0
x2i i
1x x0i
f´(x )2i
ix
3
f´(x )3i
ix= ix
4
 
Figura 3: Esquema do Método de Newton-Raphson
g(x) = f (x ji ) + f
′(x ji )(x − x ji ); f ′(x ji ) =
0− f (x ji )
x
j+1
i − x ji
⇒ x j+1i = x ji −
f (x ji )
f ′(x ji )
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 20/84 2015– 20 / 84
Operações Matemáticas Básicas Encontrar, numericamente, f (xi ) = 0
Método da Secante (aproximamos f ′)
Assume que f é localmente linear próximo a cadaxj
i
g(x) = f (xi ) + f
′(xi)(x − xi ); f ′(xi) =
0− f (xi )
xi+1 − xi
⇒ xi+1 = xi −
f (xi)
f ′(xi)
�
f
′(xi ) =
f (xi)− f (xi−1)
xi − xi−1
�
xi+1 = xi − f (xi )
xi − xi−1
f (xi )− f (xi−1)
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 21/84 2015– 21 / 84
Operações Matemáticas Básicas Optimization
Sumário
3 Operações Matemáticas Básicas
Diferenciação Numérica
Quadratura Numérica
Encontrar, numericamente, f (xi ) = 0
Optimization
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 22/84 2015– 22 / 84
Operações Matemáticas Básicas Optimization
Figura 4: William H. Press e col. Numerical Recipes in FORTRAN (1992)
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 23/84 2015– 23 / 84
Operações Matemáticas Básicas Optimization
Método de Relaxação
Relaxation Method
x0 = f (x0) Se for possível avaliar f (xi 6= x0), obtemos xi+1 6= xi , pois xi não é solução.
Expandindo f (xi) em sua série de Taylor em torno de x0:
xi+1 = f (xi ) = f (x
0)
| {z }
=x0
+(xi − x0) f ′(x0) +O[(xi − x0)2]
xi+1 − x0 = (xi − x0) f ′(x0)
Portanto a solução converge se e somente se |f ′(x0)| < 1.
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 24/84 2015– 24 / 84
Operações Matemáticas Básicas Optimization
Método de Relaxação
E se |f ′(x0)| > 1 ? Pego meu chapeu e vou embora? Para a função inversa de f , i.e.:
g = f −1, temos que:
x = f (x)→ g(x) = g(f (x)) = x
xi+1 = g(xi ) = g(x
0) + (xi − x0) g ′(x0) +O[(xi − x0)2]
xi+1 − x0 = (xi − x0) g ′(x0)
Quem é g ′?
f (g(x) = y) = x
d
dy
f
d
dx
y =
d
dx
x = 1
f
′
g
′ = 1→ g ′ = 1
f ′
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 25/84 2015– 25 / 84
Operações Matemáticas Básicas Optimization
Método de Relaxação
Quando o problema converge?
Se |f ′| < 1, converge
Se |f ′| > 1, diverge, mas . . .
Se f tem inversa:
g
′ =
1
f ′
|g ′| < 1, converge: reformule o problema em termos de g = f −1
Se f não tem inversa . . .
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 26/84 2015– 26 / 84
Operações Matemáticas Básicas Optimization
Exemplo de implementação
call rx(x0,func4,tol,itmax)
subroutine rx(x0,func,tol,itmax)
it = 0
xo=x0 ! chute inicial
do
xn = func(xo)
if (abs(xn-xo)<tol) then
print *, ’convergiu’
exit
elseif (it > itmax) then
print *, ’não convergiu’
exit
endif
xo=xn ! atualiza valor de x
it = it +1
enddo
endsubroutine
Figura 5: Solução da eq. x = f (x) pelo método de relaxação.
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 27/84 2015– 27 / 84
Operações Matemáticas Básicas Optimization
Exemplos
1) x = 2− exp−x
2) x = exp (1− x2) ou x =
√
1− ln x (inversa)
3) x = x2 + sin (2x) ou x = 1/2 sin1(x − x2)
Raphael Longuinhos (DFI/UFLA) Curso de Métodos Computacionais 28/84 2015– 28 / 84