MATLAB应用第7章第4版 求解非线性方程.docx

MATLAB应用第7章第4版 求解非线性方程.docx

《MATLAB应用第7章第4版 求解非线性方程.docx》由会员分享，可在线阅读，更多相关《MATLAB应用第7章第4版 求解非线性方程.docx（24页珍藏版）》请在冰豆网上搜索。

MATLAB应用第7章第4版求解非线性方程

第7章求解非线性方程

7.1多项式运算在MATLAB中的实现

一、多项式的表达

n次多项式表达为：

，是n+1项之和

在MATLAB中，n次多项式可以用n次多项式系数构成的长度为n+1的行向量表示

[a0,a1,……an-1,an]

二、多项式的加减运算

设有两个多项式

和

。

它们的加减运算实际上就是它们的对应系数的加减运算。

当它们的次数相同时，可以直接对多项式的系数向量进行加减运算。

当它们的次数不同时，应该把次数低的多项式无高次项部分用0系数表示。

例2计算

a=[1,-2,5,3];b=[0,0,6,-1];c=a+b

例3设

，

，求f（x）+g（x）

f=[3,-5,2,-7,5,6];g=[3,5,-3];g1=[0,0,0,g];%为了和f的次数找齐

f+g1,f-g1

三、多项式的乘法运算

conv（p1,p2）

例4在上例中，求f（x）*g（x）

f=[3,-5,2,-7,5,6];g=[3,5,-3];

conv（f,g）

四、多项式的除法运算

[Q,r]=deconv（p1,p2）

表示p1除以p2，给出商式Q（x），余式r（x）。

Q,和r仍为多项式系数向量

例4在上例中，求f（x）/g（x）

f=[3,-5,2,-7,5,6];g=[3,5,-3];

[Q,r]=deconv（f,g）

五、多项式的导函数

p=polyder（P）：

求多项式P的导函数

p=polyder（P,Q）：

求P·Q的导函数

[p,q]=polyder（P,Q）：

求P/Q的导函数，导函数的分子存入p，分母存入q。

参数P,Q是多项式的向量表示，p,q也是多项式的向量表示。

例4求有理分式

的导函数

P=[3,5,0,-8,1,-5];%有理分式分子

Q=[10,5,0,0,6,0,0,7,-1,0,-100];%有理分式分母

[p,q]=polyder（P,Q）

六、多项式求根

多项式求根就是求满足多项式p（x）＝0的x值。

N次多项式应该有n个根。

这些根可能是实根，也可能是若干对共轭复根。

其调用格式是

x=roots（P）

其中P为多项式的系数向量，求得的根赋给向量x，即x

（1）,x

（2）,…,x（n）分别代表多项式的n个根。

该命令每次只能求一个一元多项式的根，该指令不能用于求方程组的解，必须把多项式方程变成Pn（x）=0的形式；

例4求方程

的解。

首先将方程变成Pn（x）=0的形式：

roots（[1-10-1]）

例5求多项式x4+8x3-10的根。

A=[1,8,0,0,-10];

x=roots（A）

若已知多项式的全部根，则可以用poly函数建立起该多项式，其调用格式为：

P=poly（x）

若x为具有n个元素的向量，则poly（x）建立以x为其根的多项式，且将该多项式的系数赋给向量P。

例6已知f（x）=3x5+4x3-5x2-7.2x+5

（1）计算f（x）=0的全部根。

（2）由方程f（x）=0的根构造一个多项式g（x），并与f（x）进行对比。

P=[3,0,4,-5,-7.2,5];

X=roots（P）%求方程f（x）=0的根

G=poly（X）%求多项式g（x）

将这个结果乘以3，就与f（x）一致

7.2求解非线性方程f（x）=0

方程求根的一般形式是求下列方程的根：

f（x）=0（l）

实际上，就是寻找使函数f（x）等于零的变量x，所以求方程（l）的根，也叫求函数f（x）的零点。

如果变量x是列阵，则方程（l）就代表方程组。

当方程（l）中的函数f（x）是有限个指数、对数、三角、反三角或幂函数的组合时，则方程（l）被称为超越方程，例如e-x-sin（πx/2）＋lnx=0就是超越方程。

当方程（l）中的函数f（x）是多项式时，即f（x）＝Pn（x）=anxn+an-1xn＋…＋alx+a0，则方程（l）就成为下面的多项式方程，也称代数方程：

Pn（x）=anxn+an-1xn＋…＋alx+a0=0

（2）

Pn（x）的最高次数n等于2、3时，用代数方法可以求出方程

（2）的解析解，但是，当n≥5时，伽罗瓦（Galois）定理已经证明它是没有代数求根方法的。

至于超越方程，通常很难求出其解析解。

所以，方程（l）的求解经常使用作图法或数值法，而计算机的发展和普及又为这些方法提供了广阔的发展前景，使之成为科学和工程中最实用的方法之一。

本章首先介绍求解f（x）=0的MATLAB符号法指令，然后介绍求方程数值解的基本原理，最后再介绍求解f（x）=0的MATLAB数值法指令。

一、符号方程求解

在MATLAB中，求解用符号表达式表示的方程可由函数solve实现，其调用格式为：

solve（s）：

求解符号表达式s的方程，求解变量为默认变量。

当方程右端为0时，方程可以不标出等号和0，仅标出方程的左端。

solve（s,v）：

求解符号表达式s的方程，求解变量为v。

solve（s1,s2,…,sn,v1,v2,…,vn）：

求解符号表达式s1,s2,…,sn组成的方程组，求解变量分别v1,v2,…,vn。

例1．解下列方程。

1．

x=solve（'1/（x+2）+4*x/（x^2-4）=1+2/（x-2）','x'）

2．

f=sym（'x-（x^3-4*x-7）^（1/3）=1'）

x=solve（f）

3．

x=solve（'2*sin（3*x-pi/4）=1'）

4．

x=solve（'x+x*exp（x）-10','x'）%仅标出方程的左端

二、求方程f（x）=0数值解的基本方法

并非所有的方程f（x）=0都能求出精确解或解析解，不存在这种解的方程就需要用数值解法求出近似解，有几种常见的数值解法基本原理：

二分法。

1求实根的二分法原理

设方程f（x）=0中的函数f（x）为实函数，且满足：

①函数f（x）在[a,b]上单调、连续；

②方程f（x）=0在（a,b）内只有一个实根x*。

则求方程f（x）=0的根，就是在（a,b）内找出使f（x）为零的点x*：

f（x*）=0，即求函数f（x）的零点。

因为f（x）单调连续，由连续函数的性质可知，若任意两点aj，bj[a,b]，而且满足条件f（aj）f（bj）<0，则闭区间[aj,bj]上必然存在方程的根x*，即x*[aj,bj]。

据此原理提出求实根的二分法如下图所示，

图1方程求根二分法原理示意图

先用中点

将区间[a,b]平分为两个子区间（a,b1）和（b1,b），方程的根必然在子区间两端点上函数值之积小于零的那一半中，即不在（a,b1）内，就在（b1,b）内，除非f（b1）=0，于是寻根的范围缩小了一半。

图1中的根x*在区间中点左侧，即x*（a,bl）。

再将新的含根区间（a,b1）分成两半，重复上述步骤确定出更新的含根子区间。

如此重复n次，设含根区间缩小为（an,bn），则方程的根x*（an,bn）,这一系列含根的子区间满足：

（a,b）D（al,bl）（a2,b2）…（a0，b0）…

由于含根区间范围每次减半，子区间的宽度为

（n=1,2,….），显然当n→时，（bn一an）→0，即子区间收敛于一点x*，这个点就是方程的根。

若n为有限整数，取最后一个子区间的中点

作为方程根的近似值，它满足f（xn）≈0，于是有：

这就是近似值xn的绝对误差限。

假定预先要求的误差为，由

便可以求出满足误差要求的最小等分次数n。

下面是二分法的程序

function[c,err,yc]=bisect（f,a,b,delta）

%Input-fisthefunctioninputasastring‘f’

%-aandbaretheleftandrightendpoints

%.-deltaisthetolerance

%Output-cisthezero

%-yc=f（c）

%-erristheerrorestimateforc

ya=feval（f,a）;

yb=feval（f,b）;

ifya*yb>0,break,end%表示无解，结束

maxl=l+round（（log（b-a）-log（delta））/log

（2））;%从误差表达式得到最小等分次数n

fork=1:

max1

c=（a+b）/2;%取区间中点

yc=feval（f,c）;

ifyc==0

a=c;

b=c;%这时解已经找到

elseifyb*yc>0

b=c;%区间减半

yb=yc;

else

a=c;

ya=yc;

end

ifb-a

end

c=（a+b）/2;

err=abs（b-a）;

yc=feval（f,c）

2迭代法

迭代法是计算数学中的一种重要方法，用途很广，求解线性方程组和矩阵特征值时也要用到它。

这里结合非线性方程的迭代法求解，介绍一下它的基本原理。

迭代法基本原理

迭代法的基本原理就是构造一个迭代公式，反复用它得出一个逐次逼近方程根的数列，数列中每个元素都是方程根的近似值，只是精度不同。

迭代法求解方程

f（x）=0

（1）

时，先把方程等价地变换成形式

f（x）=x－g（x）=0，

（2）

移项得出：

x=g（x）（3）

若函数g（x）连续，则称（3）为迭代函数。

用它构造出迭代公式：

xk+1=g（xk）,k=0,l,2，…（4）

从初始值x0出发，便可得出迭代序列：

{xk｝=x0,x1,x2,….xk,…..（5）

如果迭代序列（5）收敛，且收敛于x*，则由式（4）有：

可见x*便是方程（l）的根。

迭代法几何意义：

如下图所示，解方程f（x）=0可以等价地变换成求解x=g（x），

图4-2方程求根迭代法原理示意图

在几何上，就等价求曲线y＝x和y＝g（x）交点P*的坐标x*。

求迭代序列（5），就等于从图中x0点出发，由函数y＝g（x0）得出y＝P0，代入函数y＝x中得出Q1，再把Q1的x坐标x1代入方程y=g（x）得出P1，如此继续下去，便可在曲线y＝g（x）上得到一系列的点P0，P1，…，Pk，…，这些点的x坐标便是迭代数列xl,x2，…，xk，…，它趋向于方程（l）的根x*，数列的元素就是方程根的近似值。

数列的收敛就等价于曲线y＝x和y＝g（x）能够相交于一点。

迭代公式收敛定理

要想用迭代法求出方程根的近似值，迭代序列（4-5）必须收敛。

下面的定理给出了迭代法的收敛条件，同时也给出了迭代公式的误差。

收敛定理：

方程x=g（x）在（a,b）内有根x*，如果：

①当x[a,b]时，g（x）[a,b];

②g（x）可导，且存在正数q<1，使得对于任意x[a,b]都有|g’（x）|q<1，则有以下结论。

①方程x=g（x）在（a,b）内有唯一的根x*。

②迭代公式xk+1=g（xk）对（a,b）内任意初始近似根x0均收敛于x*。

③近似根xk的误差估计公式为：

（4-6）

3切线法

切线法就是从函数曲线上的一点出发，不断用曲线的切线代替曲线，求得收敛于根的数列。

切线法原理：

解非线性方程f（x）=0的切线法也称牛顿法，它是把方程线性化的一种近似方法，用函数f（x）的切线代替曲线产生一个收敛于方程根的迭代序列，从而得到方程的近似根。

把函数f（x）在某一初始值x0点附近展开成泰勒级数：

（4-7）

取其线性部分，近似地代替函数f（x）可得方程的近似式：

设

，解该近似方程可得：

把函数f（x）在xl点附近展开成泰勒级数，取其线性部分替代函数f（x），设

，得：

如此继续做下去，就可以得到牛顿迭代公式：

，k=0,l,2，…（8）

由式（8）得出的迭代序列xl,x2，…，xk…，在一定的条件下收敛于方程的根x*。

2．几何意义

图4一3方程求根切线法原理示意图

选取初值x0后，过

点作曲线

的切线，其方程为

。

设切线与X釉的交点为x1，则

，再过

作切线，与x轴的交点为

，如此不断作切线，求与x轴的交点，便可得出的一系列的交点x1，x2，…，xk，…，它们逐渐逼近方程的根x*。

3．切线法的收敛性

理论可以证明，在有根区间[a,b]上，如果

、

连续且不变号，则只要选取的初始近似根x0满足

f（x。

），切线法必定收敛。

它的收敛速度经推导可得出：

（9）

是个常数，式（9）表明用牛顿迭代公式在某次算得的误差，与上次误差的平方成正比，可见牛顿迭代公式的收敛速度很快。

4.2.4割线法（弦截法）

应用切线法的牛顿迭代公式时，每次都得计算导数

，若将该导数用差商代替，就成为割线法（有时称快速弦截法）的迭代公式：

，k=0,l,2，…（4一10）

割线法的几何意义也很明显。

如图所示，

图4一4方程求根割线法原理示意图

过点（x0，f（x0））和（x1，f（x1））作函数y＝f（x）曲线的割线，交X轴于点x2，再过点（x1，f（x1））和（x2，f（x2））作曲线的割线，交X轴于点x3，……一直做下去，则得到一系列割线与X轴的交点，这些交点序列将趋于方程的根x*。

非线性方程的数值解法还有许多，这里仅介绍了几种基本方法的原理。

二分法简单方便，但收敛速度慢；

迭代法虽然收敛速度稍微快点，但需要判断能否收敛；

只要初值选取得当，切线法具有恒收敛且收敛速度快的优点，但需要求出函数的导数；

弦截法不需要求导数，特别是前面介绍的快速弦截法，收敛速度很快，但是需要知道两个近似的初始根值才能作出弦，要求的初始条件较多。

这些方法各有千秋，需根据具体情况选用。

三、方程f（x）=0数值解的MATLAB实现

MATLAB中求方程数值解的办法很多，有的是专用指令，有的是根据方程性质而借用其他专用指令求得的。

4.3.2求函数零点指令fzero

求解方程f（x）=0的实数根也就是求函数f（x）的零点。

MATLAB中设有求函数f（x）零点的指令fzero，可用它来求方程的实数根。

该指令的使用格式为：

fzero（fun,x0,options）

①输入参数fun为函数f（x）的字符表达式、内联函数名或M函数文件名。

②输入参数x0为函数某个零点的大概位置（不要取零）或存在的区间［xi，xj］，要求函数f（x）在x0点左右变号，即f（xi）f（xj）<0。

③输入参数options可有多种选择，若用optimset（'disp','iter'）代替options时，将输出寻找零点的中间数据。

④该指令无论对多项式函数还是超越函数都可以使用，但是每次只能求出函数的一个零点，因此在使用前需摸清函数零点数目和存在的大体范围。

为此，一般先用绘图指令plot,fplot或ezplot画出函数f（x）的曲线，从图上估计出函数零点的位置。

例4一5求方程x2+4sin（x）=25的实数根（-2π＜x<2π）。

解：

一、fun为函数f（x）的字符表达式

（l）首先要确定方程实数根存在的大致范围。

为此，先将方程变成标准形式f（x）=x2+4sin（x）-25=0。

作f（x）的曲线图：

x=-2*pi:

0.1:

2*pi;

f=x.^2+4*sin（x）-25;

plot（x,f）;gridon;

从曲线上可以看出，函数的零点大约在x1≈-4和x2≈5附近。

（2）直接使用指令fzero求出方程在x1≈-4时的根。

x1=fzero（'x^2+4*sin（x）-25',-4）

若键入：

fzero（'x^2+4*sin（x）-25',-4,optimset（'disp','iter'）），将显示迭代过程。

中间数据表明，求根过程中不断缩小探测范围，最后得出-4附近满足精度的近似根。

（3）求x2≈5的根：

x2=fzero（'x^2+4*sin（x）-25',5）

二、fun为函数f（x）的M函数文件名

将方程x2+4sin（x）=25编成M函数文件（实用中在函数较为复杂、而又多次重复调用时，才这样做），用fzero求解。

（1）在M文件编辑调试窗中键入：

functionyy=li4_5（x）

yy=x^2+4*sin（x）-25;

以li4_5为文件名存盘，退出编辑调试窗，回到指令窗。

（2）确定根的大体位置;

（3）在指令窗中键入下述指令可求出-4附近的根：

x1=fzero（'li4_5',-4）

键入下述指令可求出5附近的根：

x2=fzero（'li4_5',5）

三、fun为函数f（x）的内联函数名

内联函数是MATLAB提供的一个对象（Object）。

它的性状表现和函数文件一样，但内联函数的创建比较容易。

inline（'CE'）

'CE'是字符串，CE为不包含赋值符号“=”的表达式。

上式把串表达式转化为输入宗量自动生成的内联函数。

上述调用格式将自动地对CE进行辫识，把CE中由字母／数字组成的连续字符认做变量，除“预定义变量名（如i,j,pi）”和“常用函数名（如sin）”以外的由字母／数字组成的连续字符将被认做变量。

但注意：

若连续字符后紧接“左圆括号”，那么将不被当作输入宗量。

如x

（1），就不会认做输入宗量处理。

inline（'CE',arg1,arg2,…）

上述调用格式把串表达式转化为arg1,arg2等指定输入宗量的内联函数；这种调用格式是创建内联函数的最稳妥、可靠途径。

输入宗量字符可表达得更自如。

将函数f（x）写成内联函数的形式：

f=inline（'x^2+4*sin（x）-25'）

这时内联函数名为f

分别求x1≈-4和x2≈5时的根：

x1=fzero（f,-4）,x2=fzero（f,5）

例4-6求f（x）=x-10x+2=0在x0=0.5附近的根。

从f（x）的曲线看（x=-2.5:

0.01:

0.5;fx=x-10.^x+2;plot（x,fx）），

曲线的零点有两个，一个在x=-2附近，另一个在x=0.5附近

（1）建立函数文件funx.m。

（function[输出变量列表]=函数名（输入变量列表））

functionfx=funx（x）

fx=x-10.^x+2;

（2）调用fzero函数求根。

z=fzero（'funx',0.5）

例2求

在

和

作为迭代初值时的零点。

从f（x）的曲线看（x=-7:

0.01:

2;f=x-1./x+5;plot（x,f）），

曲线的零点有两个，一个在x=-5附近，另一个在x=1附近

（1）建立函数文件fz.m。

functionf=fz（x）

f=x-1./x+5;

（2）调用fz.函数求根。

fzero（'fz',-5）%以-5作为迭代初值

fzero（'fz',1）%以1作为迭代初值

7.3求解非线性方程组数值解的迭代法

一、符号方程组求解

在MATLAB中，求解用符号表达式表示的方程组仍然可由函数solve实现，其调用格式与解用符号表达式表示的方程一样。

例解下列方程组。

1．

[xy]=solve（'1/x^3+1/y^3=28','1/x+1/y=4','x,y'）

2．

[uv]=solve（'u^3+v^3=98','u+v=2','u,v'）

3．

[xy]=solve（'x+y=98','x^（1/3）+y^（1/3）=2','x,y'）

回车后出现下面的提示

Warning:

Explicitsolutioncouldnotbefound.

>InD:

\MATLAB6p5\toolbox\symbolic\solve.matline136

如果做代换：

，方程3就变成方程2，就可解

这个问题说明，符号求解并不是万能的。

如果用MATLAB得出无解或未找到所期望的解时，应该用其它方法试探求解。

4．

[xy]=solve（'x^2+y^2=5','2*x^2-3*x*y-2*y^2'）%变量由默认规则确定

二、求解非线性方程组的基本方法

对于非线性方程组（以二元方程组为例，其他可以类推）

（11）

的数值解求法，跟一元非线性方程的切线法（牛顿法）雷同，也是把非线性函数线性化，近似替代原方程得出数值解，所以也叫作牛顿迭代法。

假设方程组（11）的初始估计值为（x0，y0），可以把方程组（11）中的两个函数f1（x，y）和f2（x，y）在（x0，y0）处用二元泰勒级数展开，只取线性部分，移项得出：

（12）

若系数矩阵行列式

，则方程组（12）的解为：

，

方程组（11）中的两个函数f1（x，y）和f2（x，y）在（x1，y1）处，再用二元泰勒级数展开，只取线性部分，……如此继续替代下去，直到方程组的根达到所要求的精度，就完成了方程组的求解。

求解方程组（11）还有许多其他办法，如“最速下降法”，它是利用方程组（11）构成所谓模函数

，通过求模函数极小值的方法得到方程组的数值解，诸如此类在此不再一一列举。

三、求方程组数值解的指令

fsolve是用最小二乘法求解非线性方程组F（X）＝0的指令，变量X可以是向量或矩阵，方程组可以由代数方程或者超越方程构成。

它的使用格式为：

fsolve（'fun',X0,OPTIONS）

①参数fun是编辑并存盘的M函数文件的名称，可以用＠代替单引号对它进行标识。

M函数文件主要内容是方程F（X）=0中的函数F（X），即方程左边的函数。

②参数X0是向量或矩阵，为探索方程组解的起始点。

求解将从X0出发，逐渐趋向，最终得到满足精度要求、最接近X0的近似根X*:

F（X*）≈0。

由于X0是向量或矩阵，无法用画图方法进行估计，实际问题中常常是根