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

matlab实现牛顿迭代法求解非线性方程组已知非线性方程组如下3*x1-cos(x2*x3)-1/2=0x1^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保存到工作路径中:functionf=fun(x);%定义非线性方程组如下%变量x1x2x3%函数f1f2f3symsx1x2x3f1=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=[f1f2f3];————————————————————————————————建立函数dfun用来求方程组的雅克比矩阵将dfun.m保存到工作路径中:functiondf=dfun(x);%用来求解方程组的雅克比矩阵储存在dfun中f=fun(x);df=[diff(f,'x1');diff(f,'x2');diff(f,'x3')];df=conj(df');————————————————————————————————编程牛顿法求解非线性方程组将newton.m保存到工作路径中:functionx=newton(x0,eps,N);con=0;%其中x0为迭代初值eps为精度要求N为最大迭代步数con用来记录结果是否收敛fori=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;forj=1:length(x0);il(i,j)=x(j);endifnorm(x-x0)epscon=1;break;endx0=x;end%以下是将迭代过程写入txt文档文件名为iteration.txtfid=fopen('iteration.txt','w');fprintf(fid,'iteration');forj=1:length(x0)fprintf(fid,'x%d',j);endforj=1:ifprintf(fid,'\n%6d',j);fork=1:length(x0)fprintf(fid,'%10.6f',il(j,k));endendifcon==1fprintf(fid,'\n计算结果收敛！');endifcon==0fprintf(fid,'\n迭代步数过多可能不收敛！');endfclose(fid);————————————————————————————————运行程序在matlab中输入以下内容newton([0.10.1-0.1],0.00001,20)————————————————————————————————输出结果——————————————————————————————————————————在iteration中查看迭代过程iterationx1x2x3.mulStablePoint用不动点迭代法求非线性方程组的一个根function[r,n]=mulStablePoint(F,x0,eps)%非线性方程组：f%初始解：a%解的精度：eps%求得的一组解：r%迭代步数：nifnargin==2eps=1.0e-6;endx0=transpose(x0);n=1;tol=1;whiletolepsr=subs(F,findsym(F),x0);%迭代公式tol=norm(r-x0);%注意矩阵的误差求法，norm为矩阵的欧几里德范数n=n+1;x0=r;if(n100000)%迭代步数控制disp('迭代步数太多，可能不收敛！');return;endendx0=[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.mfunction[r,n,data]=budong(x0,tol)ifnargin=-1tol=1e-3：endx1=budongfun(x0)；n=1;while(norm(x1-x0))tol)&(n500)x0=x1;x1=budong_fun(x0)；n=n+1：data(:,n)=x1；endr=x1：new_ton.mfunction[x,n,data]=new_ton(x0,tol)ifnargin=-1tol=1e-8；endx1=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;endx=x1;budong_fun.mfunctionf=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.mfunctionf=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-Inf5.6541初始值为0，牛顿法计算结果为：x=0.5000-0.0000-0.5236初始值为1，牛顿法计算结果为：x=0.50000.0000-0.5236初始值为2，牛顿法计算结果为：x=0.50000.0000-0.5236