Newton-Raphson 法求解非线性方程组

牛顿迭代法（Newton's method）又称为牛顿-拉夫逊方法（Newton-Raphson method），它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。

设r是f(x) = 0的根，选取x0作为r初始近似值，过点（x0,f(x0)）做曲线y = f(x)的切线L，L的方程为y = f(x0) f'(x0)(x-x0)，求出L与x轴交点的横坐标 x1 = x0-f(x0)/f'(x0)，称x1为r的一次近似值。过点（x1,f(x1)）做曲线y = f(x)的切线，并求该切线与x轴的横坐标 x2 = x1-f(x1)/f'(x1)，称x2为r的二次近似值。重复以上过程，得r的近似值序列，其中x(n+1)=x(n)－f(x(n))/f'(x(n))，称为r的n+1次近似值，上式称为牛顿迭代公式。

解非线性方程f(x)=0的牛顿法是把非线性方程线性化的一种近似方法。把f(x)在x0点附近展开成泰勒级数 f(x) = f(x0)+(x－x0)f'(x0)+(x－x0)^2*f''(x0)/2! +… 取其线性部分，作为非线性方程f(x) = 0的近似方程，即泰勒展开的前两项，则有f(x0)+f'(x0)(x－x0)=f(x)=0 设f'(x0)≠0则其解为x1=x0－f(x0)/f'(x0) 这样，得到牛顿法的一个迭代序列：x(n+1)=x(n)－f(x(n))/f'(x(n))。

关于Newton's method的基本使用方法参见：http://blog.csdn.net/iamskying/archive/2009/08/22/4471071.aspx

Metlab源程序：

function homework

[P,iter,err]=newton('f','JF',[7.8e-001;4.9e-001; 3.7e-001],0.01,0.001,1000);

disp(P);

disp(iter);

disp(err);

function Y=f(x,y,z)

Y=[x^2+y^2+z^2-1;

2*x^2+y^2-4*z;

3*x^2-4*y+z^2];

function y=JF(x,y,z)

f1='x^2+y^2+z^2-1';

f2='2*x^2+y^2-4*z';

f3='3*x^2-4*y+z^2';

df1x=diff(sym(f1),'x');

df1y=diff(sym(f1),'y');

df1z=diff(sym(f1),'z');

df2x=diff(sym(f2),'x');

df2y=diff(sym(f2),'y');

df2z=diff(sym(f2),'z');

df3x=diff(sym(f3),'x');

df3y=diff(sym(f3),'y');

df3z=diff(sym(f3),'z');

j=[df1x,df1y,df1z;df2x,df2y,df2z;df3x,df3y,df3z];

y=(j);

function [P,iter,err]=newton(F,JF,P,tolp,tolfp,max)

%输入P为初始猜测值，输出P则为近似解

%JF为相应的Jacobian矩阵

%tolp为P的允许误差

%tolfp为f(P)的允许误差

%max：循环次数

Y=f(F,P(1),P(2),P(3));

for k=1:max

J=f(JF,P(1),P(2),P(3));

Q=P-inv(J)*Y;

Z=f(F,Q(1),Q(2),Q(3));

err=norm(Q-P);

P=Q;

Y=Z;

iter=k;

if (err<tolp)||(abs(Y)<tolfp||abs(Y)<0.0001)

break

end

end

来自Metlab官网的方法：

% This m-file takes care of synthetic division.
% By giving one polynomial and one root this function returns
% the polynomial formed with the other roots of the given polynomial excluding the given root.
% Keerthi Venkateswara Rao-M.Tech-NITC-INDIA
function coeff_second=syn_division(coeff_function,fun_root_new)
order_fun=size((coeff_function),2);
coeff_second=0;
for index=1:size((coeff_function),2)-1
if index==1
coeff_second(index)=coeff_function(index);
else
coeff_second(index)=coeff_function(index)+fun_root_new*coeff_second(index-1);
end
end

% This m-file calculates the derivative of the function, the limitation of
% this function is, it can calculate only the derivatives of power(x,n)....
% Keerthi Venkateswara Rao-M.tech-NITC-INDIA.
function coeff_derivative=derivate(coeff_function)
der_order=size((coeff_function),2)-1;
coeff_derivative=0;
for index=1:size((coeff_function),2)-1
coeff_derivative(index)=der_order*coeff_function(index);
der_order=der_order-1;
end

% this m-file calculates the real roots of the given polynomial using
% newton raphson technique.this m-file calls the functions in the two m-files named as syn_division and
% derivate.
% coeff_function is the 1xn matrix conatining the coeff of the polynomial.
% Keerthi Venkateswara Rao-M.Tech-NITC-INDIA
function [final_roots,functionvalue] = newton(coeff_function,initial_guess,error_tolerance,max_iterations)
iterations=0;
max_no_roots=size(coeff_function,2);
final_roots=0;
functionvalue=0;
for no_roots=1:max_no_roots-1
fun_root_new=initial_guess;
flag=1;
coeff_der_function=derivate(coeff_function);
order_fun=size(coeff_function,2);
order_der_fun=size(coeff_der_function,2);
while flag==1
fun_root_old=fun_root_new;
fx=0;
dfx=0;
nonzero=1;
while nonzero==1
powers=order_fun-1;
for index=1:order_fun
fx=fx+coeff_function(index)*fun_root_old^powers;
powers=powers-1;
end
powers=order_der_fun-1;
for index=1:order_der_fun
dfx=dfx+coeff_der_function(index)*fun_root_old^powers;
powers=powers-1;
end
if dfx==0
fun_root_old=fun_root_old+1;
else
nonzero=0;
end
end
iterations = iterations + 1;
fun_root_new = fun_root_old - fx/dfx;
if iterations >= max_iterations
flag=0;
elseif  abs(fun_root_new-fun_root_old)<=error_tolerance
flag=0;
final_roots(no_roots)=fun_root_new;
functionvalue(no_roots)=fx;
end
end
coeff_function=syn_division(coeff_function,fun_root_new);
end

参见：http://www.mathworks.com/matlabcentral/fileexchange/4313

一个C程序实例：

#include <math.h>
#include <iomanip.h>
#include <iostream.h>
#include <process.h>
void usrfun(double x[16],double alpha[16][16],double beta[16])
{
int np;
np = 15;
alpha[1][1] = -2.0 * x[1];
alpha[1][2] = -2.0 * x[2];
alpha[1][3] = -2.0 * x[3];
alpha[1][4] = 1.0;
alpha[2][1] = 2.0 * x[1];
alpha[2][2] = 2.0 * x[2];
alpha[2][3] = 2.0 * x[3];
alpha[2][4] = 2.0 * x[4];
alpha[3][1] = 1.0;
alpha[3][2] = -1.0;
alpha[3][3] = 0.0;
alpha[3][4] = 0.0;
alpha[4][1] = 0.0;
alpha[4][2] = 1.0;
alpha[4][3] = -1.0;
alpha[4][4] = 0.0;
beta[1]= x[1]*x[1]+x[2]*x[2]+x[3]*x[3]-x[4];
beta[2]=-x[1]*x[1]-x[2]*x[2]-x[3]*x[3]-x[4]*x[4]+1.0;
beta[3]=-x[1]+x[2];
beta[4]=-x[2]+x[3];
}
void main()
{
//program d10r13
//driver for routine mnewt
int i,j,kk,k,ntrial,np,n;
double tolx,tolf,xx;
ntrial = 5;
tolx = 0.000001;
n = 4;
tolf = 0.000001;
np = 15;
double x[16],alpha[16][16],beta[16];
for (kk = -1; kk<= 1; kk=kk+2)
{
for (k =1; k<=3; k++)
{
xx = 0.2 * k * kk;
cout<< "Starting vector number "<<k<<endl;
for (i = 1; i<=4; i++)
{
x[i] = xx + 0.2 * i;
cout<<"x("<<i<<")= "<<x[i]<<endl;
}
for (j = 1;j <=ntrial; j++)
{
mnewt(1,x,n,tolx,tolf);
usrfun(x,alpha,beta);
cout<<"  i        x(i)            f"<<endl;
for (i = 1; i<=n; i++)
{
cout<< setw(3)<< i;
cout<< setw(14)<< x[i];
cout<< setw(16)<< -beta[i]<<endl;
}
}
}
}
}

void mnewt(int ntrial,double x[16],int n,double tolx,double tolf)
{
double alpha[16][16],beta[16];
int k,i,indx[16];
double errf,d,errx;
for (k = 1; k<=ntrial; k++)
{
usrfun(x,alpha,beta);
errf = 0.0;
for (i = 1; i<=n; i++)
errf = errf + fabs(beta[i]);
if (errf <= tolf)  _c_exit();
ludcmp(alpha,n,indx,d);
lubksb(alpha,n,indx,beta);
errx = 0.0;
for (i = 1; i<=n; i++)
{
errx = errx + fabs(beta[i]);
x[i] = x[i] + beta[i];
}
if (errx <= tolx) _c_exit();
}
}