sábado, 14 de janeiro de 2017

Mais um pouco de Python: usando o método de Newton

Fonte aqui
O método de Newton para encontrar a raiz de uma função $F(x)$ pode ser expresso por:
$$ x_{k+1} = x_k - \frac{F(x_k)}{dF(x_k)/dx}$$
sendo $x_0$ um "chute" inicial para a raiz da função (ver aqui). Existem algumas condições para que ocorra a convergência (ver aqui). O código Python abaixo calcula a raiz da função da $f(x) = exp(-x) - (1/2)sen(x/2)$:

import numpy as np
import matplotlib.pyplot as plt

# funçao:
def ff(x):
    fx = np.exp(-x) - 0.5*np.sin(x/2);
    return fx

# derivada da função:
def fdv(x):
    h = 1e-6;
    dv = (ff(x+h) - ff(x-h))/(2*h);
    return dv

# Calculando o zero da função usando o
# Médodo de Newton com derivda numérica.
xk = 0.6;
erro = 1;
tol = 1e-12;
N = 0;
vx = [xk, 0, 0, 0, 0, 0, 0];
vf = [ff(xk), 0, 0, 0, 0, 0, 0];
while abs(erro) > tol:
    erro = ff(xk)/fdv(xk);
    xk = xk - erro;
    print(N+1, xk,erro);
    N = N + 1;
    vx[N] = xk;
    vf[N] = ff(xk);
    if (N > 10):
        erro = 0;

# Gráfico:
t = np.arange(0,20,0.01);
y = np.exp(-t) - 0.5*np.sin(t/2);
plt.plot(t,y,vx,vf,'o'); plt.grid();
plt.title('Função $ e^{-t} - (1/2)*sin(t/2) $. Localizando o 1o. zero.')
plt.show();


Resultados:
1 1.10917754392 -0.509177543917
2 1.23185126487 -0.122673720958
3 1.23768979912 -0.00583853424458
4 1.23770234859 -1.25494666513e-05
5 1.23770234864 -5.78292315977e-11
6 1.23770234864 -0.0

Gráfico:

Nenhum comentário:

Postar um comentário