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-INDIAfunction 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); endend 

% 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-INDIAfunction [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(); }}