matlab实现牛顿迭代法求解非线性方程组教学内容.doc

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

精品文档
matlab实现牛顿迭代法求解非线性方程组 
已知非线性方程组如下
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)
————————————————————————————————
输出结果 

——————————————————————————————————————————   
在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=[0 0 0];
  [r,n,data]=budong(x0);
   disp('不动点计算结果为')
   x1=[1 1 1];
   x2=[2 2 2];
   [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)
   if nargin=-1
   tol=1e-3：
   end
   x1=budong fun(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)
   if nargin=-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
   function f=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
   function f=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
   初始值为0，牛顿法计算结果为：
   x=
   0.5000   -0.0000   -0.5236
   初始值为1，牛顿法计算结果为：
   x=
   0.5000    0.0000   -0.5236
   初始值为2，牛顿法计算结果为：
   x=
   0.5000    0.0000   -0.5236
收集于网络，如有侵权请联系管理员删除