sábado, 26 de mayo de 2012

Método de Newton-Raphson

En análisis numérico, el método de Newton (conocido también como el método de Newton-Raphson o el método de Newton-Fourier) es un algoritmo eficiente para encontrar aproximaciones de losceros o raíces de una función real. También puede ser usado para encontrar el máximo o mínimo de una función, encontrando los ceros de su primera derivada.

Descripción del método

La función ƒ es mostrada en azul y la línea tangente en rojo. Vemos que xn+1 es una mejor aproximación que xn para la raíz x de la función f.
El método de Newton-Raphson es un método abierto, en el sentido de que su convergencia global no está garantizada. La única manera de alcanzar la convergencia es seleccionar un valor inicial lo suficientemente cercano a la raíz buscada. Así, se ha de comenzar la iteración con un valor razonablemente cercano al cero (denominado punto de arranque o valor supuesto). La relativa cercanía del punto inicial a la raíz depende mucho de la naturaleza de la propia función; si ésta presenta múltiples puntos de inflexión o pendientes grandes en el entorno de la raíz, entonces las probabilidades de que el algoritmo diverja aumentan, lo cual exige seleccionar un valor supuesto cercano a la raíz. Una vez que se ha hecho esto, el método linealiza la función por la recta tangente en ese valor supuesto. La abscisa en el origen de dicha recta será, según el método, una mejor aproximación de la raíz que el valor anterior. Se realizarán sucesivas iteraciones hasta que el método haya convergido lo suficiente. f'(x)= 0 Sea f : [a,b-> R función derivable definida en el intervalo real [ab]. Empezamos con un valor inicial x0 y definimos para cada número natural n
x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.
Donde f ' denota la derivada de f.
Nótese que el método descrito es de aplicación exclusiva para funciones de una sola variable con forma analítica o implícita cognoscible. Existen variantes del método aplicables a sistemas discretos que permiten estimar las raíces de la tendencia, así como algoritmos que extienden el método de Newton a sistemas multivariables, sistemas de ecuaciones, etc.

Código en C

Programa escrito en C correspondiente al ejemplo f(x) = cos(x) - x3.
#include <stdio.h>
#include <math.h>
 
double FuncionIteracion(double x);
 
int main() {
   int i;
   double x, xvieja;
 
   xvieja = 0.5; // valor inicial
 
   printf("Iteracion      valor\n");
   for(i=0; i<8; i++){
      x = FuncionIteracion(xvieja);
      xvieja = x;
      printf("%5.0d  %20.12f\n", i+1, x);
   }
 
   return 0;
}
 
double FuncionIteracion(double x)
{
   double valorFuncion, funcionOriginal, funcionDerivada;
 
   funcionOriginal = cos(x) - pow(x, 3);
   funcionDerivada = -sin(x) - 3*pow(x, 2);
   valorFuncion = x - funcionOriginal / funcionDerivada;
 
   return valorFuncion;
}
Para compilar en GNU/Linux con compilador de GNU, se escribe en una terminal:
$ gcc programa.c -lm -o programa
donde la opción -lm le indica al compilador que utilice la biblioteca de matemáticas.

[editar]Codigo en Matlab

Programa escrito en Matlab para hallar las raíces usando el método de NEWTON-RAPHSON
% Al escribir la función, usar x como variable.
x0=input('Ingrese el valor inicial: ');
tol=input('Ingrese el porcentaje de error: ');
f=input('Ingrese la función: ');
i=1;
fx(i)=x0;
syms x;
f1=subs(f,x,fx(i));
z=diff(f);
d=subs(z,x,fx(i));
ea(1)=100;
while abs(ea(i))>=tol;
    fx(i+1)=fx(i)-f1/d; f1=subs(f,x,fx(i+1)); d=subs(z,x,fx(i+1));
    ea(i+1)=abs((fx(i+1)-fx(i))/fx(i+1)*100);
    i=i+1;
end
fprintf('i     fx(i)         Error aprox (i) \n');
for j=1:i;
    fprintf('%2d \t %11.7f \t %7.3f \n',j-1,fx(j),ea(j));
end
El programa siguiente hace el cálculo para una superficie.
syms x1
syms x2
syms x3
V=['sin(x1)+2^x2+log(x3)-7';'3*x1+2*x2-x3^3+1      ';'x1+x2+x3-5            '];
%se calcula el jacobiano:
DV(1,:)=[diff(V(1,:),x1), diff(V(1,:),x2),diff(V(1,:),x3)];
DV(2,:)=[diff(V(2,:),x1), diff(V(2,:),x2),diff(V(2,:),x3)];
DV(3,:)=[diff(V(3,:),x1), diff(V(3,:),x2),diff(V(3,:),x3)];
%se da el valor de partida:
x1=0;
x2=4;
x3=2;
x_1o=[x1;x2;x3];
%se calcula H en ese punto
Vo(1,:)=eval(V(1,:));
Vo(2,:)=eval(V(2,:));
Vo(3,:)=eval(V(3,:));
%Se calcula el Hessiano en ese punto
DV1=eval(DV);
%se calcula la Inversa del Hessiano en ese punto
DV_1=DV1^-1;
%se calcula el siguiente valor de iteración
x_1=[x1;x2;x3]-DV_1*Vo;
%cantidad de iteraciones maxima:
n=50; 
%se define a = n, si se cumple condicion de error antes,  cambia
a=n; 
 
for i=1:n  
%error relativo entre aproximaciones sucecivas
    er=norm(x_1-x_1o)/norm(x_1);
    if er<.0001
        a=i;
       'break;' end
 
    x1=x_1(1);
    x2=x_1(2);
    x3=x_1(3);
    x_1o=[x1;x2;x3];
    Vo(1,:)=eval(V(1,:));
    Vo(2,:)=eval(V(2,:));
    Vo(3,:)=eval(V(3,:));
    DV1=eval(DV);
    DV_1=DV1^-1;
    x_1=[x1;x2;x3]-DV_1*Vo;  
 
end
    a
    x_1

Otros documento en PDF del mismo método:
http://www.epsem.upc.edu/~fpq/ale-hp/modulos/aplicaciones/newton.pdf
http://www.mty.itesm.mx/etie/deptos/m/ma00-130/lecturas/m130-18.pdf
http://www.mty.itesm.mx/etie/deptos/m/ma00-130/lecturas/m130-18.pdf

No hay comentarios:

Publicar un comentario