Implementación de algoritmos/Matemáticas/Método de Newton
El algoritmo del Método de Newton consiste en una relación de recurrencia.
Sea f: [a, b] -> R función derivable definida en el intervalo real [a, b]. Empezamos con un valor inicial x0 y definimos para cada número natural n
Donde f ' denota la derivada de f.
Ejemplo
editarConsideremos el problema de encontrar un número positivo x tal que cos(x) = x3. Podríamos tratar de encontrar el cero de f(x) = cos(x) - x3.
Sabemos que f '(x) = -sin(x) - 3x2. Ya que cos(x) ≤ 1 para todo x y x3 > 1 para x>1, deducimos que nuestro cero está entre 0 y 1. Comenzaremos probando con el valor inicial x0 = 0,5
Los dígitos correctos están subrayados. En particular, x6 es correcto para el número de decimales pedidos. Podemos ver que el número de dígitos correctos después de la coma se incrementa desde 2 (para x3) a 5 y 10, ilustando la convergencia cuadrática.
En pseudocódigo, esto es:
function newtonIterationFunction(x) {
return x - (cos(x) - x^3) / (-sin(x) - 3*x^2)
}
var x := 0,5
for i from 0 to 99 {
print "Iteraciones: " + i
print "Valor aproximado: " + x
xold := x
x := newtonIterationFunction(x)
if x = xold {
print "Solución encontrada!"
break
}
}
Código en Python
editarPrograma escrito en Python correspondiente al ejemplo f(x) = cos(x) - x3.
import math
def newtonIterationFunction(x):
return x - (math.cos(x) - x**3) / (-math.sin(x) - 3*(x**2))
x = 0.5
for i in range(100):
print "Iteraciones: ",str(i),"Valor aproximado: ", str(x)
xold = x
x = newtonIterationFunction(x)
if x == xold:
print "Solución encontrada!"
break
Para ejecutar escribe en una terminal:
$ python nombre-del-programa.py
Código en C
editarPrograma 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.
Código en Matlab
editarPrograma 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