MATLAB应用求解非线性方程

1、第7章求解非线性方程7.1 多项式运算在MATLAB中的实现一、多项式的表达n次多项式表达为：p(x) =a0xn a/n-1an-1x - an，是n+1项之和在MATLAB中，n次多项式可以用n次多项式系数构成的长度为n+1的行 向量表示a0, a1,an-1,a n二、多项式的加减运算设 有两个 多项式 p1(x)二a0xn a/"-1an-1x - an和p2(x)二boXm b1Xm-1 bm-1X bm。它们的加减运算实际上就是它们的对应系数的加减运算。当它们的次数相同时，可以直接对多项式的系数向量进行加减运 算。当它们的次数不同时，应该把次数低的多项式无高次项部分用0系

2、数表示。例 2 计算 x3 -2x2 5x 3 6x-1a=1, -2, 5, 3; b=0, 0, 6,-1; c=a+b例 3 设 f x = 3x5 -5x° 2x3 -7x2 5x 6 , g x = 3x2 5x -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;con v(f, g)四、多项式的除法运算Q, r=d

3、econv(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=dec onv(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求有理分式f x二的导函数3x5 + 5x4 8x? + x 一510x10 5x9 6x6 7x3

4、 x? 一100%有理分式分子%有理分式分母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个根。该命令每次只能求一个一元多项式的根，该指令不能用于求方程组的解，必须把多项式方程变成Pn (x) = 0的形式；例4求方程x3 = x2

5、1的解首先将方程变成Pn (x) = 0的形式：X3 - X2 - 1 = 0roots(1 -1 0 -1)例5求多项式x4+8x3-10的根。A=1,8,0,0,-10; x=roots(A)poly函数建立起该多项式，其调用格式若已知多项式的全部根，贝U可以用为：P=poly(x)若x为具有n个元素的向量，则poly(x)建立以x为其根的多项式，且将该 多项式的系数赋给向量P。例 6 已知 f(x)=3x5+4x3-5x2-7.2x+5(1) 计算f(x)=0的全部根。由方程f(x)=0的根构造一个多项式g(x)，并与f(x)进行对比。P=3,0,4,-5,-7.2,5;X=roots(

6、P)G=poly(X)%求方程f(x)=0的根%求多项式g(x)将这个结果乘以3,就与f(x)一致7.2求解非线性方程f ( x ) = 0方程求根的一般形式是求下列方程的根：f ( x ) = 0(l)实际上，就是寻找使函数f ( x)等于零的变量x，所以求方程(I)的根, 也叫求函数f ( x)的零点。如果变量x是列阵，则方程(I )就代表方程组。当方程(I )中的函数f (x)是有限个指数、对数、三角、反三角或幕函数 的组合时，则方程(I)被称为超越方程，例如 e-x - sin (n / 2 ) + Inx = 0就是 超越方程。当方程(I)中的函数f (X)是多项式时，即f ( x)

7、= Pn (x) = anxn + an-1X° + ax + ao，则方程(I)就成为下面的多项式方程，也称代数方程：Pn (x) = anxn + an-ixn + + ax + ao = 0( 2 )Pn (x)的最高次数n等于2、3时，用代数方法可以求出方程(2)的解析 解，但是，当n > 5时，伽罗瓦(GaIois)定理已经证明它是没有代数求根方法 的。至于超越方程，通常很难求出其解析解。所以，方程(I )的求解经常使用 作图法或数值法，而计算机的发展和普及又为这些方法提供了广阔的发展前景， 使之成为科学和工程中最实用的方法之一。本章首先介绍求解f ( x ) = 0

8、的MATLAB符号法指令，然后介绍求方程数 值解的基本原理，最后再介绍求解f ( x ) = 0的MATLAB数值法指令。一、符号方程求解在MATLAB中，求解用符号表达式表示的代数方程可由函数soIve实现,其调用格式为：solve(s):求解符号表达式s的代数方程，求解变量为默认变量。 当方程右端为0时，方程可以不标出等号和0,仅标出方程的左端。soIve(s,v)：求解符号表达式s的代数方程，求解变量为v。soIve(s1,s2, ,sn,v1,v2,:求解符号表达式 s1,s2,sfi成的代数方程 组，求解变量分别v1,v2,vn例1.解下列方程。丄玉胡丄2 1x 2 x2 -4 x

9、-2x= soIve('1/(x+2)+4*x/(xA2-4)=1+2/(x-2)', 'x')2. x -即 x2sin 3x 一 | = 1<4丿x= soIve('2*si n(3*x-pi/4)=1') x xex -10=0x= solve('x+x*exp(x)-10', 'x')% 仅标出方程的左端二、求方程f( x ) = 0数值解的基本方法并非所有的方程f ( x ) = 0都能求出精确解或解析解，不存在这种解的方程 就需要用数值解法求出近似解，有几种常见的数值解法基本原理：二分法。1求实根

10、的二分法原理设方程f (x) =0中的函数f ( x)为实函数，且满足： 4x - 7 =1f=sym('x-(xA3-4*x-7)A(1/3)=1') x= soIve(f) 函数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)f (bj) < 0，则闭区间aj , bj上必然存 在方程的

11、根x*，即x* aj , bj。据此原理提出求实根的二分法如下图所示，22根必然在子区间两端点上函数值之积小于零的那一半中，即不在(a ,bi)内，就在(bi ,b )内，除非f(bi) = 0，于是寻根的范围缩小了一半。图1中的根x*在区 间中点左侧，即x*(a , b)。再将新的含根区间(a , bi)分成两半，重复上述步 骤确定出更新的含根子区间。如此重复n次，设含根区间缩小为(an, bn)，贝昉程的根( an, bn),这一系列含根的子区间满足：(a , b ) D 二(ai ,bi )二(a2,b2 )二 二：(ao，bo)二：由于含根区间范围每次减半，子区间的宽度为0 - a.二

12、宁 (n = 1,2,.)，显然当n-：时，(bn一 an)-0,即子区间收敛于一点x*，这个点就是方程的根。若n为有限整数，取最后一个子区间的中点作为方程根的近似值，它2满足f ( Xn)旳，于是有:*1 b a b a& 一X近亍二尹这就是近似值Xn的绝对误差限。假定预先要求的误差为:，由；：:瞎便可以求2 n-n出满足误差要求的最小等分次数n。下面是二分法的程序fun cti on c,err,yc =bisect (f,a,b,delta)%ln put - f is the fun cti on in put as a stri ng '% - a and b are

13、 the left and right end points%. - delta is the tolera nce%Output - c is the zero% - yc=f(c)% - err is the error estimate for cya=feval (f,a);yb=feval (f,b);if ya*yb>0, break, end%表示无解，结束maxl=l+round( (log (b-a) -log (delta)/log (2); % 从误差表达式得到最小等分次数nfor k=1:max1c=(a+b)/2; %取区间中点yc=feval (f,c);if

14、 yc=0a=c;b=c; %这时解已经找到elseif yb*yc>0b=c; %区间减半yb=yc;else?a=c;ya=yc;endif b-a < delta, break, endendc=(a+b)/2;err=abs(b-a);yc=feval (f, c)2迭代法迭代法是计算数学中的一种重要方法，用途很广，求解线性方程组和矩阵特 征值时也要用到它。这里结合非线性方程的迭代法求解，介绍一下它的基本原理。迭代法基本原理迭代法的基本原理就是构造一个迭代公式，反复用它得出一个逐次逼近方程根的数列，数列中每个元素都是方程根的近似值，只是精度不同。迭代法求解方程f ( X )

15、 = 0(1)时，先把方程等价地变换成形式f ( X ) = x g(x) = 0 ,(2)移项得出：x = g(x)( 3 )若函数g (x)连续，则称(3)为迭代函数。用它构造出迭代公式：xk+i= g ( xk), k = 0,1,2 ,(4 )从初始值xo出发，便可得出迭代序列： x k = Xo, xi, X2,.Xk,.(5 )如果迭代序列(5 )收敛，且收敛于x*，则由式(4)有：pm g x -Xki 二 g x* -x* 二 f x* =0可见x*便是方程(I)的根。迭代法几何意义：如下图所示，解方程f ( X ) = 0可以等价地变换成求解 x = g ( x )，在几何上

16、，就等价求曲线y=x和y = g ( x )交点P*的坐标x*。求迭代序列(5)， 就等于从图中xo点出发，由函数y= g ( xo)得出y = Po，代入函数y=x中得出 Qi，再把Qi的x坐标xi代入方程y= g ( x )得出Pi，如此继续下去，便可在曲 线y=g ( x)上得到一系列的点Po，Pi，Pk，这些点的x坐标便是迭 代数列XI , x2 ，Xk，它趋向于方程(I )的根X*，数列的元素就 是方程根的近似值。数列的收敛就等价于曲线 y= x和y=g ( x )能够相交于一 点。迭代公式收敛定理要想用迭代法求出方程根的近似值，迭代序列(4 - 5)必须收敛。下面的定 理给出了迭代

17、法的收敛条件，同时也给出了迭代公式的误差。收敛定理：方程x = g ( 乂)在(a , b )内有根x*，如果： 当 x a,b时，g( x)a,b; g ( x)可导，且存在正数q < i，使得对于任意x a,b都有|g'(x )|乞q < i,则有以下结论。 方程x = g ( x)在(a , b)内有唯一的根x*。 迭代公式xk+i = g ( xk)对(a , b)内任意初始近似根 xo均收敛于x*。 近似根Xk的误差估计公式为：kXn X* 兰一|xi Xo(4 - 6)i -q3切线法切线法就是从函数曲线上的一点出发，不断用曲线的切线代替曲线，求得 收敛于根的数

18、列。切线法原理：解非线性方程f( x ) = 0的切线法也称牛顿法，它是把方程线性化的一种近 似方法，用函数f (x)的切线代替曲线产生一个收敛于方程根的迭代序列，从而 得到方程的近似根。把函数f ( x)在某一初始值x。点附近展开成泰勒级数：2 f "fx f x = f x0 X - Xo f XoX - Xo - 丄(4 - 7)2!取其线性部分，近似地代替函数f (x)可得方程的近似式：f X : f XoX - Xo f Xo =0设f Xo - o，解该近似方程可得:f X。f Xo把函数f ( X )在Xl点附近展开成泰勒级数，取其线性部分替代函数f (x)，设 f X

19、i o，得：f(Xi )X2 =Xi -f X1如此继续做下去，就可以得到牛顿迭代公式:xk 1= xkf Xk由式(8)得出的迭代序列Xl , X2 ,(8 ),Xk，在一定的条件下收敛于f(Xi )方程的根X2 .几何意义(忌，f (Xo)/J 3Y图4 一 3方程求根切线法原理示意图选取初值Xo后，过Xo, f Xo 点作曲线y = f x的切线，其方程为y - fX。= x - XofXo。设切线与X釉的交点为xi，则捲=Xo- fXo，再过f(X。)Xf为作切线，与X轴的交点为X2 = X1 _ f X1 ，如此不断作切线，求与X f(Xi )轴的交点，便可得出的一系列的交点 X1，

20、X2，Xk，它们逐渐逼近方程的 根X*。3 .切线法的收敛性理论可以证明，在有根区间a, b上,如果f Xo -o、r Xo -o连续且不变 号，则只要选取的初始近似根xo满足f Xo f Xo o f ( X。)，切线法必定收敛。 它的收敛速度经推导可得出：f "(x* ) #2Xk i -X*Xk - x*( 9)2f x*f是个常数，式（9）表明用牛顿迭代公式在某次算得的误差，与上次误2f x*差的平方成正比，可见牛顿迭代公式的收敛速度很快。4.2.4割线法（弦截法）应用切线法的牛顿迭代公式时，每次都得计算导数 f Xk，若将该导数用差商代替，就成为割线法（有时称快速弦截法）的

21、迭代公式：Xk 1 = xk -xk xkd.fXk- f Xkf Xk ，（4 一 10 ）割线法的几何意义也很明显。如图 所示,过点（X0，f（ Xo）和（XI，f （Xi）作函数y = f（x ）曲线的割线，交X轴于点X2, 再过点（Xi，f （Xi）和（X2，f（X2）作曲线的割线，交X轴于点X3，一直做下去，则得到一系列割线与 X轴的交点，这些交点序列将趋于方程的根x*。非线性方程的数值解法还有许多，这里仅介绍了几种基本方法的原理。二分法简单方便，但收敛速度慢；迭代法虽然收敛速度稍微快点，但需要判断能否收敛；只要初值选取得当，切线法具有恒收敛且收敛速度快的优点， 但需要求出函 数的导

22、数；弦截法不需要求导数，特别是前面介绍的快速弦截法，收敛速度很快，但是 需要知道两个近似的初始根值才能作出弦，要求的初始条件较多。这些方法各有千秋，需根据具体情况选用。三、方程 f(x) = 0 数值解的 MATLAB 实现MATLAB 中求方程数值解的办法很多，有的是专用指令，有的是根据方程 性质而借用其他专用指令求得的。4 . 3 . 2 求函数零点指令 fzero求解方程 f ( x ) = 0 的实数根也就是求函数 f ( x )的零点。 MATLAB 中设有 求函数f (x)零点的指令fzero,可用它来求方程的实数根。该指令的使用格式为： fzero (fun, x0, optio

23、ns) 输入参数fun为函数f (x)的字符表达式、内联函数名或M函数文件名。 输入参数X0为函数某个零点的大概位置(不要取零)或存在的区间Xi , Xj,要求函数f(X)在X0点左右变号，即f (Xi)f (Xj) < 0。 输入参数options可有多种选择，若用 optimset ('disp', 'iter')代替options 时，将输出寻找零点的中间数据。 该指令无论对多项式函数还是超越函数都可以使用,但是每次只能求出 函数的一个零点,因此在使用前需摸清函数零点数目和存在的大体范围。为此, 一般先用绘图指令plot, fplot或ezplot画

24、出函数f (x)的曲线，从图上估计出函 数零点的位置。2例 4 一 5 求方程 x + 4sin(x) = 25 的实数根(-2 n<x < 2 n)。 解：一、fun为函数f (x。的字符表达式(l。首先要确定方程实数根存在的大致范围。为此，先将方程变成标准形式 f(x) = x2 + 4sin(x) - 25 = 0。作 f(x)的曲线图： x=-2*pi:0.1:2*pi; f=x.A2+4*si n(x)-25;plot(x,f);grid on;从曲线上可以看出，函数的零点大约在X1-4和X25附近。(2) 直接使用指令fzero求出方程在xi4时的根。x1= fzero

25、 ('xA2+4*sin(x)-25',-4)若键入：fzero ('xA2+4*sin(x)-25',-4, optimset('disp', 'iter')，将显示迭代过程。 中间数据表明，求根过程中不断缩小探测范围，最后得出 - 4附近满足精度 的近似根。(3) 求X2 "5的根：x2= fzero ('xA2+4*sin(x)-25',5)二、fun为函数f (x。的M函数文件名将方程 x2 + 4sin(x) = 25 编成 M 函数文件(实用中在函数较为复杂、而又 多次重复调用时，才这样做)，

26、用fzero求解。(1) 在 M 文件编辑调试窗中键入：function yy = li4_5(x)yy= xA2+4*sin(x)-25;以 li4_5 为文件名存盘，退出编辑调试窗，回到指令窗。(2) 确定根的大体位置 ;(3) 在指令窗中键入下述指令可求出 - 4 附近的根： x1= fzero (' li4_5',-4)键入下述指令可求出 5 附近的根：x2= fzero (' Ii4_5',5)三、fun为函数f (x)的内联函数名内联函数是MATLAB提供的一个对象(Object)。它的性状表现和函数文件 一样，但内联函数的创建比较容易。in li

27、ne('CE')'CE'是字符串，CE为不包含赋值符号“=”的表达式。上式把串表达式转化为输入宗量自动生成的内联函数。上述调用格式将自动地对CE进行辫识，把CE中由字母/数字组成的连续 字符认做变量，除“预定义变量名(如 i , j , pi ) ”和“常用函数名(如sin )” 以外的由字母/数字组成的连续字符将被认做变量。但注意：若连续字符后紧接“左圆括号”，那么将不被当作输入宗量。如x ( 1 )，就不会认做输入宗量处理。inlin e('CE',arg1,arg2,)上述调用格式把串表达式转化为arg1,arg2等指定输入宗量的内联函数；

28、这 种调用格式是创建内联函数的最稳妥、可靠途径。输入宗量字符可表达得更自如。将函数f (x)写成内联函数的形式：f = in li ne ('xA2+4*si n(x)-25')这时内联函数名为f分别求xi q- 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.Ax+2; plot(x,fx)，曲线的零点有两个，一个在x=-2附近，另一个在x=0.5附近(1) 建立函数文件funx.m。(functi

29、on 输出变量列表=函数名(输入变量列表) fun cti on fx=fun x(x)fx=x-10.Ax+2; 调用fzero函数求根。z=fzero('f un x',0.5)“ 1例2求f x二x-在x。=5和x。=1作为迭代初值时的零点。x从f(x)的曲线看(x=-7:0.01:2; f= x-1./x+5; plot(x,f)，曲线的零点有两个，一个在建立函数文件fz.m。 fun cti on f=fz(x) f= x-1./x+5;(2)调用fz.函数求根。 fzero('fz',-5) fzero('fz',1)x=-5附近，另

30、一个在x=1附近%以-5作为迭代初值%以1作为迭代初值7. 3求解非线性方程组数值解的迭代法、符号方程组求解在MATLAB中，求解用符号表达式表示的方程组仍然可由函数 其调用格式与解用符号表达式表示的方程一样。例解下列方程组。solve实现,1 1x y4、一 x yx y=solve('1/xA3+1/yA3=28', '1/x+1/y=4', 'x,y') I3 +v3 =982.丿u +v = 2u v=solve('uA3+vA3=98', 'u+v=2', 'u,v')；x + y = 9

31、83.対x +V7=2x y=solve('x+y=98', 'xA(1/3)+yA(1/3)=2', 'x,y')回车后出现下面的提示Warning: Explicit solution could not be found.> In D:MATLAB6p5toolboxsymbolicsolve.m at li ne 136如果做代换：x. u3,yv3，方程3就变成方程2，就可解这个问题说明，符号求解并不是万能的。如果用MATLAB得出无解或未找到所期望的解时，应该用其它方法试探求 解。4.丿x2 y2 =52 2 2x -3xy-2

32、y =Ox y=solveCxA2+yA2=5', '2*xA2-3*x*y-2*yA2'）% 变量由默认规则确定二、求解非线性方程组的基本方法对于非线性方程组（以二元方程组为例，其他可以类推）(11)；fi（x, y ）=O 2（x,y ）=O的数值解求法，跟一元非线性方程的切线法（牛顿法）雷同，也是把非线性函数 线性化，近似替代原方程得出数值解，所以也叫作牛顿迭代法。假设方程组（ii）的初始估计值为（xo，yo），可以把方程组（ii）中的两 个函数fi（x，y）和f2（x，y）在（xo，yo）处用二元泰勒级数展开，只取线性部分， 移项得出：也(Xo,y°L

33、、丄 ofi(xo, yo,、(12):(x X。)+(y y。)= fi(X。, yo ) exdyf2(Xo,yO)(X- Xo )+ 芯(Xo,yo)(y_yo)=_ f? (x° , y° )&xdy若系数矩阵行列式Jo汴xo,yo；:xfXo, yo:x戲i(x°, y°)弓的2(x°, y°)cy=0 ,贝U方程组（ 的解为:新i(x°, yo )1exJocf2(Xo, yo )fl Xo, yof2 Xo, yoXi - Xo 'iyi = yo ' Jofi Xo, yof2 Xo,

34、 yoy）和f2（x, y）在（xi, yi）处，再用二元泰勒方程组（11）中的两个函数fi（x，级数展开，只取线性部分，如此继续替代下去，直到方程组的根达到所要求的精度，就完成了方程组的求解。求解方程组（11）还有许多其他办法，如“最速下降法”，它是利用方程组（11）构成所谓模函数门X二l-f1 x, y F f2 x, y 2，通过求模函数极小值的方法 得到方程组的数值解，诸如此类在此不再一一列举。三、求方程组数值解的指令fsolve是用最小二乘法求解非线性方程组F （X ） = O的指令，变量X可 以是向量或矩阵，方程组可以由代数方程或者超越方程构成。它的使用格式为：fsolve ( &

35、#39;fun' , X0 , OPTIONS ) 参数fun是编辑并存盘的M函数文件的名称，可以用代替单引号对 它进行标识。M函数文件主要内容是方程 F ( X ) = 0中的函数F (X)，即方程 左边的函数。 参数X0是向量或矩阵，为探索方程组解的起始点。求解将从X0出发， 逐渐趋向，最终得到满足精度要求、最接近 X0的近似根X* : F ( X* ) 0。由 于X0是向量或矩阵，无法用画图方法进行估计，实际问题中常常是根据专业知识、物理意义等进行估计。 该指令输出一个与X0同维的向量或矩阵，为方程组的近似数值解。 参数OPTIONS为设置选项，用它可以设置过程显示与否、误差、算

36、法，具体内容可用help查阅。通常可以省略该项内容。f 22x + y +z + 7=10x例 求方程组 xy2 2z在X0= 1，y°= 1，Z0= 1附近的数值解。x2 + y2 +z2 = 3yI解(l)在文本编辑调试窗中编辑 M函数文件。首先将方程组变换成F ( X )"x2 _10x+y2 +z+7=0=0的形式，* xy2-2z=0x, y, z看成向量X的三个分量。2 2 2x + y 3y+z =0fun cti on ms=li4_7(X)ms(1) = X(1).A2 - 10* X(1) + X(2).A2 + X(3) +7;ms(2) = X(1)

37、.*X(2).A2 - 2* X(3);ms(3) = X(1).A2 + X(2)八2 - 3* X(2) + X(3).A2;这里，X (l )= x , X (2)= y , X (3)= z，输出参量 ms也有三个分量用“Ii4_7”为M函数文件名存盘。(2) 在指令窗中键入：fsolve (' li4_7',1 1 1)若键入：x = fsolve (li4_7,1 1 1, optimset('Display', 'iter')则得出求解过程。一该方程也可以用 MATLAB的符号指令solve求解，但结果非常冗长。3x =co&yz )+0.5例 求解方程组 <2x2 81(y+0.1f+