matlab實現牛頓迭代法求解非線性方程組

http://hi.baidu.com/aillieo/blog/item/0800e2a10ac9a59647106493.html

已知非線性方程組如下
3*x1-cos(x2*x3)-1/2=0

x1^2-81*(x2+0.1)^2+sin(x3)+1.06=0

exp(-x1*x2)+20*x3+(10*pi-3)/3=0

求解要求精度達到0.00001

————————————————————————————————

首先建立函數fun

儲存方程組編程如下將fun.m儲存到工作路徑中:

function f=fun(x);

%定義非線性方程組如下

%變數x1 x2 x3

%函數f1 f2 f3

syms x1 x2 x3

f1=3*x1-cos(x2*x3)-1/2;

f2=x1^2-81*(x2+0.1)^2+sin(x3)+1.06;

f3=exp(-x1*x2)+20*x3+(10*pi-3)/3;

f=[f1 f2 f3];

————————————————————————————————

建立函數dfun

用來求方程組的雅克比矩陣將dfun.m儲存到工作路徑中:

function df=dfun(x);

%用來求解方程組的雅克比矩陣儲存在dfun中

f=fun(x);

df=[diff(f,'x1');diff(f,'x2');diff(f,'x3')];

df=conj(df');

————————————————————————————————

編程牛頓法求解非線性方程組將newton.m儲存到工作路徑中:

function x=newton(x0,eps,N);

con=0;

%其中x0為迭代初值eps為精度要求N為最大迭代步數con用來記錄結果是否收斂

for i=1:N;

    f=subs(fun(x0),{'x1' 'x2' 'x3'},{x0(1) x0(2) x0(3)});

    df=subs(dfun(x0),{'x1' 'x2' 'x3'},{x0(1) x0(2) x0(3)});

    x=x0-f/df;

    for j=1:length(x0);

        il(i,j)=x(j);

    end

    if norm(x-x0)<eps

        con=1;

        break;

    end

    x0=x;

end

 

%以下是將迭代過程寫入txt文檔檔案名稱為iteration.txt

fid=fopen('iteration.txt','w');

fprintf(fid,'iteration');

for j=1:length(x0)

    fprintf(fid,'         x%d',j);

end

for j=1:i

    fprintf(fid,'\n%6d     ',j);

    for k=1:length(x0)

        fprintf(fid,' %10.6f',il(j,k));

    end

end

if con==1

    fprintf(fid,'\n計算結果收斂！');

end

if con==0

    fprintf(fid,'\n迭代步數過多可能不收斂！');

end

fclose(fid);

————————————————————————————————

運行程式

在matlab中輸入以下內容

newton([0.1 0.1 -0.1],0.00001,20)

————————————————————————————————

輸出結果

ans =

 

    0.5000    0.0000   -0.5236

 

——————————————————————————————————————————————————

在iteration中查看迭代過程

 

iteration         x1         x2         x3

     1        0.490718   0.031238  -0.519661

     2        0.509011   0.003498  -0.521634

     3        0.500928   0.000756  -0.523391

     4        0.500227   0.000076  -0.523550

     5        0.500019   0.000018  -0.523594

     6        0.500005   0.000002  -0.523598

     7        0.500000   0.000000  -0.523599

計算結果收斂！