matlab实现牛顿迭代法求解非线性方程组.doc
《matlab实现牛顿迭代法求解非线性方程组.doc》由会员分享，可在线阅读，更多相关《matlab实现牛顿迭代法求解非线性方程组.doc（5页珍藏版）》请在冰点文库上搜索。
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）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=[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
初始值为0，牛顿法计算结果为：
x=
0.5000 -0.0000 -0.5236
初始值为1，牛顿法计算结果为：
x=
0.5000 0.0000 -0.5236
初始值为2，牛顿法计算结果为：
x=
0.5000 0.0000 -0.5236
升级会员 