matlab实现牛顿迭代法求解非线性方程组

matlab实现牛顿迭代法求解非线性方程组matlab实现牛顿迭代法求解非线性方程组/NUMPAGES10matlab实现牛顿迭代法求解非线性方程组matlab实现牛顿迭代法求解非线性方程组matlab实现牛顿迭代法求解非线性方程组

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

x1^2-81*(x2+0.1)^2+sin(x3)+1.06=0exp(-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);endfor

j=1:ifprintf(fid,'\n%6d

',j);for

k=1:length(x0)fprintf(fid,'

%10.6f',il(j,k));end

endif

con==1fprintf(fid,'\n计算结果收敛！');

endif

con==0fprintf(fid,'\n迭代步数过多可能不收敛！');endfclose(fid);————————————————————————————————运行程序在matlab中输入以下内容

newton([0.1

0.1

-0.1],0.00001,20)————————————————————————————————输出结果

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

在iteration中查看迭代过程

iteration

x1

x2

x3

.mulStablePoint用不动点迭代法求非线性方程组的一个根

function

[r,n]=mulStablePoint(F,x0,eps)

%非线性方程组：f

%初始解：a

%解的精度：eps

%求得的一组解：r

%迭代步数：n

if

nargin==2

eps=1.0e-6;

end

x0

=

transpose(x0);

n=1;

tol=1;

while

tol>eps

r=

subs(F,findsym(F),x0);

%迭代公式

tol=norm(r-x0);

%注意矩阵的误差求法，norm为矩阵的欧几里德范数

n=n+1;

x0=r;

if(n>100000)

%迭代步数控制

disp('迭代步数太多，可能不收敛！');

return;

end

end

x0=[000];

[r,n,data]=budong(x0);

disp('不动点计算结果为')

x1=[111];

x2=[222];

[x,n,data]=new_ton(x0);

disp(’初始值为0，牛顿法计算结果为：’)

[x,n,data]=new_ton(x1);

disp('初始值为1，牛顿法计算结果为：')

[x,n,data]=new_ton(x2);

disp('初始值为2，牛顿法计算结果为：')

budong.m

function[r,n,data]=budong(x0,tol)

ifnargin=-1

tol=1e-3：

end

x1=budongfun(x0)；

n=1;

while(norm(x1-x0))tol)&(n500)

x0=x1;

x1=budong_fun(x0)；

n=n+1：

data(:,n)=x1；

end

r=x1：

new_ton.m

function[x,n,data]=new_ton(x0,tol)

ifnargin=-1

tol=1e-8；

end

x1=x0-budong_fun(x0)/df1(x0);

n=1;

while(norm(x1-x0))tol)

x0=x1;

x1=x0-budong_fun(x0)/df1(x0);

n=n+1;

data(:,n)=x1;

end

x=x1;

budong_fun.m

functionf=budong_fun(x)

f(1)=3*x(1)-cos(x(2)*x(3))-1/2;

f(2)=x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;

f(3)=exp(-x(1)*x(2))+20*x(3)+10*pi/3-1;

f=[f(1)*f(2)*f(3)]；

df1.m

functionf=df1(x)

f=[3sin(x(2)*x(3))*x(3)sin(x(2)*x(3))*x(2)

2*x(1)-162*(x(2)+0.1)cos(x(3))

exp(-x(1)*x(2))*(-x(2))exp(-x(1)*x(2))*(-x(1))20];

结果：

不动点计算结果为

r=

1.0e+012*

NaN

-Inf

5.6541