matlab非线性方程求解解读.docx
《matlab非线性方程求解解读.docx》由会员分享，可在线阅读，更多相关《matlab非线性方程求解解读.docx（18页珍藏版）》请在冰豆网上搜索。
matlab非线性方程求解解读
非线性方程的解法
1引言
数学物理中的许多问题归结为解函数方程的问题，即，
f（x）=O（1.1）
这里，f（x）可以是代数多项式，也可以是超越函数。
若有数x为方程f（x）=0的根，或称函
数f（x）的零点。
设函数f（x）在［a,b］内连续，且f（a）f（b）<0。
根据连续函数的性质知道，方程f（x）=O在
区间［a,b］内至少有一个实根；我们又知道，方程f（x）=O的根，除了极少简单方程的根可以用解析式表达外，一般方程的根很难用一个式子表达。
即使能表示成解析式的，往往也很复杂，不便计算。
所以，具体求根时，一般先寻求根的某一个初始近似值，然后再将初始近似值逐步加工成满足精度要求为止。
如何寻求根的初始值呢？
简单述之，为了明确起见，不妨设f（x）在区间［a,b］内有一个实
的单根，且f（a）<0,f（b）0。
我们从左端出点冷二a出发，按某一预定的步长h一步一步地向右跨，每跨一步进行一次根的搜索”即检查每一步的起点Xk和xui（即，Xk+h）的函数值是否同号。
若有：
f（Xk）*f（Xkh）EO（1.2）
那么所求的根必在（Xk，xkh）内，这时可取Xk或Xkh作为根的初始近似值。
这种方法通常称为定步长搜索法”另外，还是图解法、近似方程法和解析法。
2迭代法
2.1迭代法的一般概念
迭代法是数值计算中一类典型方法，不仅用于方程求根，而且用于方程组求解，矩阵求
特征值等方面。
迭代法的基本思想是一种逐次逼近的方法。
首先取一个精糙的近似值，然后
用同一个递推公式，反复校正这个初值，直到满足预先给定的精度要求为止
对于迭代法，一般需要讨论的基本问题是：
迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。
这里，主要看看解方程迭代式的构造。
对方程（1.1），在区间［a,b］内，可改写成为：
x=%x）（2.1）
取X。
•[a,b]，用递推公式：
Xk1二:
（xk）,k=0,1,2,（2.2）
可得到序列：
X0NX2,Xk,二｛xjki（2.3）
当k》:
:
时，序列｛兀｝萬有极限~，且（x）在~附近连续，则在式（2.2）两边极限，得,
~x
-
*
X
*X
-
k
m:
x
即，~为方程（2.1）的根。
由于方式（1.1）和方程（2.1）等价，所以,
式（2.2）称为迭代式，也称为迭代公式；（x）可称为迭代函数。
称求得的序列｛xk｝k出为迭代序
列。
2.2程序和实例
下面是基于MATLAB的迭代法程序，用迭代格式png（xn），求解方程x=g（x），其中
初始值为po。
**************************************************************************
function[p,k,err,P]=fixpt（f1021,p0,tol,max1）
%f1021是给定的迭代函数。
%p0是给定的初始值。
%tol是给定的误差界。
%maxi是所允许的最大迭代次数。
%k是所进行的迭代次数加1。
%p是不动点的近似值。
%err是误差。
%P={p1,p2,…,pn}
P⑴=p0；
fork=2:
max1
P（k）=feval（'f1021',P（k-1））;
k,err=abs（P（k）-P（k-1））
p=P（k）;
if（errvtol）,
break;
end
ifk==max1
disp（'maximumnumberofiterationsexceeded'）;end
end
P=P;
****************************************************************************
例2.1用上述程序求方程sinx-x—0的一个近似解，给定初始值x^0.5，误差界为
10^。
解：
先用m文件先定义一个名为f1021.m的函数文件。
functiony=f1021（x）
y=sin（x）/x;
建立一个主程序prog1021.m
clcclearfixpt（
all
'f1021',0.5,10A（-5）,20）
然后在MATLAB命令窗口运行上述主程序，即:
>>prog1021
计算结果如下。
k=3
err=0.1052
k=4
err=0.0292
k=5
err=0.0078
k=6
err=0.0021
k=7
err=5.7408e-004
k=8
err=1.5525e-004
k=9
err=4.1975e-005
k=10
err=1.1350e-005
k=11
err=3.0688e-006
P=Columns1through6
0.5000
0.9589
0.8537
0.8829
0.8751
0.8772
Columns7through11
0.8766
0.8768
0.8767
0.8767
0.8767
ans=0.8767
3二分法
3.1二分法原理
二分法是方程求解最直观、最简单的方法。
二分法以连续函数的介值定理为基础的。
由
介值定理知道，若函数f（x）区间［a,b］上连续，且f（a）*f（b）：
：
：
0，即f（a）和f（b）负号相反,
则f（x）在［a,b］内一定有实根。
二分法的基本思想是：
用对分区间的方法根据分点处函数f（x）
的符号逐步将有限区间缩小，使在足够小的区间内，方程有且仅有一根。
下面简述其基本步骤。
首先记a^a.bo=b。
用中点x。
二也直将区间［a°,bo］等分成2个小区间：
⑻凶］和
2
［xo,bo］。
然后分析可能存在的三种情况：
如果f（a）*f（X。
）=0，则x。
是零点，也就是方程的根。
如果f（a）*f（xo）：
：
0，则区间［ao,x。
］内存在零点。
如果f（xo）*f（b）:
:
0，则区间［xo,bo］内存在零点。
对有根的新区间施行同样的操作，于是得到一系列有空的区间:
［a°,bo］二g’b』二a©］二二［ak，bk］
（3.1）
其中每1个区间的长度都是前一区间长度的一半，最后1个区间的长度为:
b「akb「a
2k
如果取最后1个区间［ak,bk］的中点:
（3.2）
bk—ak
x“2
（3.3）
作为f（x）=0根的近似值，则有误差估计式:
...bk理...b-a
对于所给精度；，若取k使得
b—a
则有，
（3.4）
（3.5）
（3.6）
3.2程序与实例
用二分法求解方程f（x）=0在有根区间［a,b］内的一个根，其中f（x）在［a,b］只有一个根的
情形
***************************************************************************
function[c,err,yc]=bisect（f1031,a,b,delta）
%f1031是所要求解的函数。
%a和b分别为有根区间的左右限。
%delta是允许的误差界。
%c为所求近似解。
%yc为函数f在c上的值。
%err是c的误差估计。
ya=feval（'f1031',a）;
yb=feval（'f',b）;
ifyb==0,c=b；returnendifya*yb>0,
disp（'（a,b）不是有根区间'）;return
end
max1=1+round（（log（b-a）-log（delta））/log
（2））;
fork=1:
max1
c=（a+b）/2;
yc=feval（'f',c）;ifyc==0a=c;b=c;
returnelseifyb*yc>0b=c;
yb=yc;elsea=c;ya=yc;
end
if（b-a）end
k;
c=（a+b）/2;
err=abs（b-a）;
yc=feval（'f',c）;
*****************************************************************************
例3.1用上述程序求方程x‘—x一1=0在区间［1.0,1.5］内的一个近似解，要求准确到小数
点后第2位。
解：
先用m文件先定义一个名为f1031.m的函数文件。
functiony=f1031（x）y=xA3-x-1;
建立一个主程序prog1031.m
clc
clearall
bisect（'f1031',1.0,1.5,0.0005）
然后在MATLAB命令窗口运行上述主程序，即:
>>prog1031
计算结果如下。
ans=
1.3247
4牛顿法
4.1牛顿法原理
从前面迭代法，我们知道，迭代函数（x）构造的好坏，不仅影响收敛速度，而且迭代格
式有可能发散。
怎样选择一个迭代函数才能保证迭代序列一定收敛呢？
构代迭代函数的一条重要途径是用近似方程来代替原方程去求根。
因此，如果能将非线性方程（1.1）用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。
牛顿（Newton）法就是一种将非线性方程线化的一种方法。
设Xk是方程（1.1）的一个近似根，把如果f（x）在Xk处作一阶Taylor展开，即：
f（x）：
f（xk）f'（Xk）（x-Xk）（4.1）
于是我们得到如下近似方程：
f（xjf'（Xk）（x-Xk）=0（4.2）
设f'（Xk）=0，则方程（10.2.1）的解为:
取~作为原方程（1.1）的新近似根X-，即令:
f'（Xk）
上式称为牛顿迭代格式。
用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法
牛顿法具有明显的几何意义。
方程：
（1045）
y=f（兀）f'（Xk）（x-Xk）
是曲线y=f（X）上点（Xk，f（Xk））处的切线方程。
迭代格式（4.4）就是用切线式（4.5）的零点来代替曲线（1.1）的零点。
正因为如此，牛顿法也称为切线法。
牛顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。
一般来说，牛顿法对初值Xo的要求较高，初值足够靠近X*时才能保证收敛。
若要保证初值在较大范围内收敛，则需对f（X）加一些条件。
如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛顿法作一些改时，即可以采用下面的迭代格式：
Xk1=Xk，k=0,1,2，（4.6）
f（Xk）
式中，0：
：
：
'：
：
：
1，称为下山因子。
因此，用这种方法求方程的根，也称为牛顿下山法。
牛顿法对单根收敛速度快，但每迭代一次，除需计算f（xj之外，还要计算f'（Xk）的值。
如果f（X）比较复杂，计算f'（Xk）的工作量就可能比较大。
为了避免计算导数值，我们可用差商来代替导数。
通常用如下几种方法：
（1）割线法。
如果用
f（Xk）-f（Xk4）
Xk_Xk4
代替f'（Xk），贝U得到割线法的迭代格式为:
（4.7）
⑵拟牛顿法。
如果用
f（xk）-■f（xk—■f（Xk」））f（Xk）
代替f'（Xk），贝U得到拟牛顿法的迭代格式为:
xk1
f2（Xk）
f（Xk）-f做-f（Xk』）
（4.8）
⑶Steffenson法。
如果用
f（Xkf（Xk））-f（Xk）
f（Xk）
代替f'（Xk），贝U得到拟牛顿法的迭代格式为:
4.2程序与实例
1.牛顿法的程序
f（x）=0。
给定初值po，用牛顿法格式Pk-Pw」12山，k=1,2/，求解非线性方程f'（P」
**************************************************************************
function[p1,err,k,y]=newton（f1041,df1041,p0,delta,max1）
%f1041是非线性函数。
%df1041是f1041的微商。
%p0是初始值。
%delta是给定允许误差。
%max1是迭代的最大次数。
%p1是牛顿法求得的方程的近似解。
%err是p0的误差估计。
%k是迭代次数。
%y=f（p1）
p0,feval（'f1041',p0）
fork=1:
max1
pl=p0-feval（'f1041',p0）/feval（'df1041',pO）;
err=abs（pl-pO）;
pO=p1;
p1,err,k,y=feval（'f1041',p1）
if（errbreak,
end
p1,err,k,y=feval（'f1041',p1）
end
****************************************************************************
例4.1用上述程序求方程x3-3x•2=0的一个近似解，给定初值为1.2,误差界为10$
解：
先用m文件先定义二个名为f1041.m和df1041.m的函数文件。
functiony=f1041（x）y=xA3-3*x+2;
functiony=df1041（x）
y=3*xA2-3;
建立一个主程序prog1041.m
clear
newton（'f1041','df1041',1.2,10八（-6）,18）
然后在MATLAB命令窗口运行上述主程序，即：
>>prog1041
计算结果如下。
p0=1.2000
ans=0.1280
p1=1.1030
err=0.0970
k=1
y=0.0329
p1=1.1030
err=0.0970
k=1
y=
0.0329
p1=
1.0524
err=
0.0507
k=
2
y=
0.0084
p1=
1.0524
err=
0.0507
k=
2
y=
0.0084
p1=
1.0264
err=
0.0260
k=
3
y=
0.0021
p1=
1.0264
err=
0.0260
k=
3
y=
0.0021
p1=
1.0133
err=
0.0131
k=
4
y=
5.2963e-004
p1=
1.0133
err=
0.0131
k=
4
y=
5.2963e-004
p1=
1.0066
err=
0.0066
k=
5
y=
1.3270e-004
p1=
1.0066
err=
0.0066
k=
5
y=
1.3270e-004
p1=
1.0033
err=
0.0033
k=
6
y=
3.3211e-005
p1=
1.0033
err=
0.0033
k=
6
y=
3.3211e-005
p1=
1.0017
err=
0.0017
k=
7
y=
8.3074e-006
p1=
1.0017
err=
0.0017
k=
7
y=
8.3074e-006
p1=
1.0008
err=
8.3157e-004
k=
8
y=
2.0774e-006
p1=
1.0008
err=
8.3157e-004
k=
8
y=
2.0774e-006
p1=
1.0004
err=
4.1596e-004
k=
9
y=
5.1943e-007
p1=
1.0004
err=
4.1596e-004
k=
9
y=
5.1943e-007
p1=
1.0002
err=
2.0802e-004
k=
10
y=
1.2987e-007
p1=
1.0002
err=
2.0802e-004
k=
10
y=
1.2987e-007
p1=
1.0001
err=
1.0402e-004
k=
11
y=
3.2468e-008
p1=
1.0001
err=
1.0402e-004
k=
11
y=
3.2468e-008
p1=
1.0001
err=
5.2014e-005
k=
12
y=
8.1170e-009
p1=
1.0001
err=
5.2014e-005
k=
12
y=
8.1170e-009
p1=
1.0000
err=
2.6008e-005
k=
13
y=
2.0293e-009
p1=
1.0000
err=
2.6008e-005
k=
13
y=
2.0293e-009
p1=
1.0000
err=
1.3004e-005
k=
14
y=
5.0732e-010
p1=
1.0000
err=
1.3004e-005
k=
14
y=
5.0732e-010
p1=
1.0000
err=
6.5020e-006
k=
15
y=
1.2683e-010
p1=
1.0000
err=
6.5020e-006
k=
15
y=
1.2683e-010
p1=
1.0000
err=
3.2510e-006
k=
16
y=
3.1708e-011
p1=
1.0000
err=3.2510e-006
k=
16
y=
3.1708e-011
p1=
1.0000
err=
1.6255e-006
k=
17
y=
7.9270e-012
p1=
1.0000
err=
1.6255e-006
k=
17
y=
7.9270e-012
p1=
1.0000
err=
8.1277e-007
k=
18
y=
1.9817e-012
ans
1.0000
这说明，经过18次迭代得到满足精度要求的值。
2.割线法的MATLAB实现
****************************************************************************
function[p1,err,k,y]=secant（f1042,p0,p1,delta,max1）%f1042是给定的非线性函数。
%p0,p1为初始值。
%delta为给定误差界。
%max1是迭代次数的上限。
%p1为所求得的方程的近似解。
%err为p1-p0的绝对值。
%k为所需的迭代次数。
%y=f（pi）
p0,p1,feval（'f1042',p0）,feval（'f1042',p1）,k=0;
fork=1:
max1
p2=pl-feval（'f1042',p1）*（p1-p0）/（feval（'f1042',p1）-feval（'f1042',p0））;err=abs（p2-p1）;
p0=p1;
pl=p2;
p1,err,k,y=feval（'f1042',p1）；
if（errbreak,
end
end
******************************************************************************
例4.2用上述程序求方程x‘-x•2=0的一个近似解，给定初值为-1.5,-1.52,误差界为
10B
解：
先用m文件先定义一个名为f1042.m的函数文件。
functiony=f1042（x）
y=xA3-x+2;
建立一个主程序prog1042.m
clc
clear
secant（'f1042',-1.5,-1.52,10八（-6）,11）
然后在MATLAB命令窗口运行上述主程序，即：
>>prog1042
计算结果如下。
p0=
-1.5000
p1=
-1.5200
ans=
0.1250
ans=
0.0082
p1=
-1.5214
err=
0.0014
k=
1
p1=
-1.5214
err=2.2961e-005
k=2p1=-1.5214err=2.4318e-008
k=3
-1.5214。
ans=-1.5214
这就表明，经过3次迭代得到满足精度要求的的近似解
升级会员 