牛頓迭代法（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))。

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 

一個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(); } }