MATLAB應用求解非線性方程

1、實用標準第 7 章 求解非線性方程7.1多項式運算在 MATLAB中的實現(xiàn)一、多項式的表達n 次多項式表達為：p(x)a0 x na1 x n -1an-1 xan ，是 n+1 項之和在 MATLAB 中，n 次多項式可以用 n 次多項式系數(shù)構成的長度為 n+1 的行向量表示a0, a1,an-1,an二、多項式的加減運算設有兩 個多項式 p1(x) ax nax n-1an-1xa和01np2(x)b0 x mb1x m-1b m-1xbm 。它們的加減運算實際上就是它們的對應系數(shù)的加減運算。 當它們的次數(shù)相同時， 可以直接對多項式的系數(shù)向量進行加減運算。當它們的次數(shù)不同時，應該把次數(shù)低的

2、多項式無高次項部分用0 系數(shù)表示。例 2 計算 x32x 25x36x1a=1, -2, 5, 3; b=0, 0, 6, -1; c=a+b例3設f x3x 55x 42x 37x 25x6 ，g x3x25x3 ， 求f(x)+g(x)f=3, -5, 2, -7, 5, 6; g=3, 5, -3; g1=0, 0, 0, g;%為了和 f 的次數(shù)找齊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=dec

3、onv(p1, p2)表示 p1 除以 p2 ，給出商式 Q(x) ，余式 r(x) 。Q,和 r 仍為多項式系數(shù)向量例 4 在上例中，求 f(x)/g(x) f=3, -5, 2, -7, 5, 6; g=3, 5, -3;Q, r=deconv(f, g)五、多項式的導函數(shù)p=polyder(P)：求多項式 P 的導函數(shù)p=polyder(P,Q) ：求 PQ 的導函數(shù)p,q=polyder(P,Q) ：求 P/Q 的導函數(shù)，導函數(shù)的分子存入p ，分母存入q 。參數(shù) P,Q 是多項式的向量表示，p,q 也是多項式的向量表示。例 4 求有理分式 f x3x55x 48x2x 55x 96 x

4、67x 3的導函數(shù)10 x10x 2 100P=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 個根。這些根可能是實根，也可能是若干對共軛復根。其調(diào)用格式是x=roots(P)其中 P 為多項式的系數(shù)向量， 求得的根賦給向量x，即 x(1),x(2),x(n) 分別代表多項式的 n 個根。該命令每次只能求一個一元多項式的根，該指令不能用于求方程組的解，必須把多項

5、式方程變成Pn (x) = 0 的形式；例 4 求方程 x 3x 21 的解。首先將方程變成 Pn (x) = 0 的形式： x3x 210roots(1 -1 0 -1)例 5 求多項式 x4+8x 3-10 的根。A=1,8,0,0,-10;x=roots(A)若已知多項式的全部根，則可以用poly 函數(shù)建立起該多項式，其調(diào)用格式為：P=poly(x)若 x 為具有 n 個元素的向量，則 poly(x) 建立以 x 為其根的多項式，且將該多項式的系數(shù)賦給向量 P。例 6 已知 f(x)=3x 5 +4x 3 -5x 2-7.2x+5(1) 計算 f(x)=0 的全部根。(2) 由方程 f(

6、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)實際上，就是尋找使函數(shù)f ( x）等于零的變量x，所以求方程（ l）的根，也叫求函數(shù)f ( x ）的零點。如果變量x 是列陣，則方程（ l ）就代表方程組。當方程（ l ）中的函數(shù)f (x）是有限個指數(shù)、對數(shù)、三角、反三角或冪函數(shù)的組合時，則方程（ l

7、）被稱為超越方程，例如e-x - sin （x / 2 )lnx = 0就是超越方程。當方程（l）中的函數(shù) f（ x）是多項式時， 即 f（ x） Pn（x）= a n xn + a n-1 xn alx + a 0 ，則方程（ l）就成為下面的多項式方程，也稱代數(shù)方程：nnn+ an-1xnl0= 0( 2 )P（ x） = ax a x + aPn（x）的最高次數(shù) n 等于 2 、3時，用代數(shù)方法可以求出方程（2 ）的解析解，但是，當 n 5時，伽羅瓦（ Galois ）定理已經(jīng)證明它是沒有代數(shù)求根方法的。至于超越方程，通常很難求出其解析解。所以，方程（l ）的求解經(jīng)常使用作圖法或數(shù)值法，

8、而計算機的發(fā)展和普及又為這些方法提供了廣闊的發(fā)展前景，使之成為科學和工程中最實用的方法之一。本章首先介紹求解f ( x ) = 0的 MATLAB符號法指令，然后介紹求方程數(shù)值解的基本原理，最后再介紹求解f ( x ) = 0的 MATLAB數(shù)值法指令。文檔實用標準一、符號方程求解在 MATLAB 中，求解用符號表達式表示的代數(shù)方程可由函數(shù)solve 實現(xiàn)，其調(diào)用格式為：solve(s) ：求解符號表達式s 的代數(shù)方程，求解變量為默認變量。當方程右端為 0 時，方程可以不標出等號和0，僅標出方程的左端。solve(s,v) ：求解符號表達式s 的代數(shù)方程，求解變量為v。solve(s1,s2,

9、sn,v1,v2,vn) ：求解符號表達式s1,s2,sn 組成的代數(shù)方程組，求解變量分別v1,v2,vn 。例 1 解下列方程。14x21x241x 2x 2x= solve(1/(x+2)+4*x/(x2-4)=1+2/(x-2), x)2 x3x 34x71f=sym(x-(x3-4*x-7)(1/3)=1)x= solve(f)3 2 sin 3x41x= solve(2*sin(3*x-pi/4)=1)4 x xex100x= solve(x+x*exp(x)-10, x)% 僅標出方程的左端二、求方程 f ( x ) = 0數(shù)值解的基本方法并非所有的方程f ( x ) = 0都能求

10、出精確解或解析解，不存在這種解的方程就需要用數(shù)值解法求出近似解，有幾種常見的數(shù)值解法基本原理：二分法。文檔實用標準1 求實根的二分法原理設方程f (x) =0 中的函數(shù)f ( x ）為實函數(shù)，且滿足： 函數(shù) f (x ）在 a , b 上單調(diào)、連續(xù)； 方程 f (x) = 0在（ a , b ）內(nèi)只有一個實根x* 。則求方程f (x) = 0的根，就是在（ a, b ）內(nèi)找出使 f (x ）為零的點 x* ：f (x*) =0 ，即求函數(shù)f ( x )的零點。因為f (x ）單調(diào)連續(xù)，由連續(xù)函數(shù)的性質(zhì)可知，若任意兩點 aj ， bj a , b，而且滿足條件f (aj) f (bj) 0, b

11、reak, end%表示無解，結束maxl=l+round( (log (b-a) -log (delta)/log (2); %從誤差表達式得到最小等分次數(shù) nfor k=1:max1c=(a+b)/2; %取區(qū)間中點yc=feval (f,c);if yc=0a=c;b=c; % 這時解已經(jīng)找到elseif yb*yc0b=c; % 區(qū)間減半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 迭代法迭代法是計算數(shù)學中的一種重要方法，用途很廣，求解線性

12、方程組和矩陣特征值時也要用到它。這里結合非線性方程的迭代法求解， 介紹一下它的基本原理。迭代法基本原理迭代法的基本原理就是構造一個迭代公式，反復用它得出一個逐次逼近方程根的數(shù)列，數(shù)列中每個元素都是方程根的近似值，只是精度不同。迭代法求解方程f ( x ) = 0（1）時，先把方程等價地變換成形式f ( x ) = x g(x) = 0 ，（2 ）移項得出：x = g(x)( 3 )若函數(shù) g (x ）連續(xù)，則稱（ 3）為迭代函數(shù)。用它構造出迭代公式：xk+1 = g ( x k ） , k = 0 , l , 2，( 4 )從初始值 x0 出發(fā)，便可得出迭代序列： x k = x 0, x 1

13、, x2 , .xk, .（ 5 ）如果迭代序列（ 5）收斂，且收斂于x* ，則由式（ 4 ）有：lim g x xk 1g x * x *f x *0k可見 x* 便是方程（ l）的根。迭代法幾何意義：如下圖所示，解方程f ( x ) = 0可以等價地變換成求解x = g ( x ) ，文檔實用標準圖 4 - 2 方程求根迭代法原理示意圖在幾何上，就等價求曲線 y x 和 y g ( x ）交點 P* 的坐標 x* 。求迭代序列（5) ，就等于從圖中 x0 點出發(fā)，由函數(shù) yg ( x 0）得出 y P0 ，代入函數(shù) y x 中得出Q1 ，再把 Q 1 的 x 坐標 x1 代入方程 y= g

14、 ( x）得出 P1，如此繼續(xù)下去，便可在曲線 y g ( x ) 上得到一系列的點P0 ，P1 ， Pk ，這些點的 x 坐標便是迭代數(shù)列xl , x2 ， x k，它趨向于方程（ l ) 的根 x* ，數(shù)列的元素就是方程根的近似值。數(shù)列的收斂就等價于曲線y x 和 y g ( x）能夠相交于一點。迭代公式收斂定理要想用迭代法求出方程根的近似值，迭代序列（4 - 5 ）必須收斂。下面的定理給出了迭代法的收斂條件，同時也給出了迭代公式的誤差。收斂定理：方程x = g ( x）在（ a , b）內(nèi)有根x* ，如果： 當 xa,b 時， g( x ）a,b; g ( x ）可導，且存在正數(shù)q 1

15、，使得對于任意 xa,b 都有 |g ( x )|q 1 ，則有以下結論。 方程 x = g ( x ）在（ a , b ）內(nèi)有唯一的根 x* 。 迭代公式xk+1 = g ( x k）對（ a , b ）內(nèi)任意初始近似根x0 均收斂于x* 。文檔實用標準 近似根xk 的誤差估計公式為：q kx1 x0(4-6)xn x *1q3 切線法切線法就是從函數(shù)曲線上的一點出發(fā)，不斷用曲線的切線代替曲線，求得收斂于根的數(shù)列。切線法原理：解非線性方程 f( x ) = 0的切線法也稱牛頓法，它是把方程線性化的一種近似方法，用函數(shù)f (x ）的切線代替曲線產(chǎn)生一個收斂于方程根的迭代序列，從而得到方程的近似

16、根。把函數(shù) f ( x ）在某一初始值x 。點附近展開成泰勒級數(shù)：f xf x0x x0 f x0x x02 fx0(4-7)2!文檔實用標準取其線性部分，近似地代替函數(shù) f(x ）可得方程的近似式：f xf x0x x0 f x00設 fx00 ，解該近似方程可得：x1fx0x0x0f把函數(shù) f ( x）在 xl 點附近展開成泰勒級數(shù)，取其線性部分替代函數(shù)f (x) ，設 fx10 ，得：x2fx1x1x1f如此繼續(xù)做下去，就可以得到牛頓迭代公式：f xk，k = 0 , l , 2，(8)xk 1 xkxkf由式（ 8）得出的迭代序列 xl , x2， xk，在一定的條件下收斂于方程的根

17、x* 。2 幾何意義圖 4 一 3 方程求根切線法原理示意圖文檔實用標準選取初值 x0 后，過x0 , f x0 點作曲線 yfx 的切線，其方程為y fx0xx0 fx0。設切 線與 X 釉的交點為 x1 ，則 x1x0f x0 ，再過f x0x1 , fx1作切線，與 x 軸的交點為 x2x1fx1， 如此不斷作切線，求與 xfx1軸的交點，便可得出的一系列的交點x1 ，x2 ， ，xk， ，它們逐漸逼近方程的根 x* 。3 切線法的收斂性理論可以證明，在有根區(qū)間 a, b 上，如果 fx00 、 fx00 連續(xù)且不變號，則只要選取的初始近似根 x0 滿足 fx0 fx00 f ( x。）

18、，切線法必定收斂。它的收斂速度經(jīng)推導可得出：xk 1x*f x *xk x * 2( 9)2 f x *f x * 是個常數(shù)，式(9) 表明用牛頓迭代公式在某次算得的誤差， 與上次誤2 fx *差的平方成正比，可見牛頓迭代公式的收斂速度很快。4 . 2 . 4割線法（弦截法）應用切線法的牛頓迭代公式時，每次都得計算導數(shù)fxk，若將該導數(shù)用差商代替，就成為割線法（有時稱快速弦截法）的迭代公式：xk 1xkxk 1f xk ，k = 0 , l , 2，(4一10)xkf xk 1f xk割線法的幾何意義也很明顯。如圖所示，文檔實用標準p0p1p2p3p4x4x3x2x1x 0圖 4 一 4 方程

19、求根割線法原理示意圖過點（ x0，f( x 0 ）和（x1 ，f (x 1)）作函數(shù) yf(x ）曲線的割線，交 X 軸于點 x2，再過點（ x1 ，f (x 1 )）和（ x2，f (x 2 )）作曲線的割線，交X 軸于點 x3 ， 一直做下去，則得到一系列割線與X 軸的交點，這些交點序列將趨于方程的根x* 。非線性方程的數(shù)值解法還有許多，這里僅介紹了幾種基本方法的原理。二分法簡單方便，但收斂速度慢；迭代法雖然收斂速度稍微快點，但需要判斷能否收斂；只要初值選取得當， 切線法具有恒收斂且收斂速度快的優(yōu)點，但需要求出函數(shù)的導數(shù)；弦截法不需要求導數(shù)， 特別是前面介紹的快速弦截法，收斂速度很快， 但

20、是需要知道兩個近似的初始根值才能作出弦，要求的初始條件較多。這些方法各有千秋，需根據(jù)具體情況選用。文檔實用標準三、方程 f(x) = 0數(shù)值解的 MATLAB實現(xiàn)MATLAB 中求方程數(shù)值解的辦法很多，有的是專用指令，有的是根據(jù)方程性質(zhì)而借用其他專用指令求得的。4 . 3 . 2求函數(shù)零點指令fzero求解方程 f ( x ) = 0的實數(shù)根也就是求函數(shù)f ( x ）的零點。 MATLAB 中設有求函數(shù) f (x ）零點的指令 fzero ，可用它來求方程的實數(shù)根。該指令的使用格式為：fzero (fun, x0 , options) 輸入?yún)?shù) fun 為函數(shù) f (x ）的字符表達式、內(nèi)聯(lián)函

21、數(shù)名或 M 函數(shù)文件名。 輸入?yún)?shù) x0 為函數(shù)某個零點的大概位置 （不要取零） 或存在的區(qū)間 x i，xj，要求函數(shù) f (x ）在 x0 點左右變號，即 f (x i )f (x j) 0 。 輸入?yún)?shù) options可有多種選擇，若用optimset (disp, iter）代替options時，將輸出尋找零點的中間數(shù)據(jù)。 該指令無論對多項式函數(shù)還是超越函數(shù)都可以使用，但是每次只能求出函數(shù)的一個零點，因此在使用前需摸清函數(shù)零點數(shù)目和存在的大體范圍。為此，一般先用繪圖指令plot, fplot或 ezplot畫出函數(shù) f (x ）的曲線，從圖上估計出函數(shù)零點的位置。例 4 一 5 求方程

22、x2 + 4sin(x) = 25的實數(shù)根（ -2 x In D:MATLAB6p5toolboxsymbolicsolve.m at line 136如果做代換： xu3 , yv3 ，方程 3 就變成方程 2，就可解這個問題說明，符號求解并不是萬能的。如果用 MATLAB 得出無解或未找到所期望的解時，應該用其它方法試探求解。4x 2y2523xy2 y202xx y=solve(x2+y2=5,2*x2-3*x*y-2*y2)% 變量由默認規(guī)則確定二、求解非線性方程組的基本方法對于非線性方程組（以二元方程組為例，其他可以類推）f1x, y0f2x, y( 11)0的數(shù)值解求法，跟一元非線

23、性方程的切線法（牛頓法）雷同，也是把非線性函數(shù)線性化，近似替代原方程得出數(shù)值解，所以也叫作牛頓迭代法。文檔實用標準假設方程組（ 11 ）的初始估計值為（ x0 ，y0），可以把方程組（ 11 ）中的兩個函數(shù) f1 (x，y) 和 f 2 (x，y)在（ x 0， y0 ）處用二元泰勒級數(shù)展開，只取線性部分，移項得出：f1x0, y0f1x0, y0y y0f1 x0 , y0xx x0y(12)f 2x0, y0f2x0 , y0y y0f2 x0 , y0xx x0yf1 x0 , y0f1 x0 , y0若系數(shù)矩陣行列式 J0xy0 ，則方程組（ 12 ）的解為：f2 x0 , y0f 2

24、 x0 , y0xyf1 x0, y0f1 x0 , y01 f1 x0 , y0f1 x0 , y0x1x01xy0xf2x0, y0， y1J0 f 2 x0 , y0f2x0 , y0f 2 x0 , y0J0yy方程組（ 11 ）中的兩個函數(shù)f1 (x ，y)和 f2 (x，y) 在(x1， y1 )處，再用二元泰勒級數(shù)展開，只取線性部分， 如此繼續(xù)替代下去，直到方程組的根達到所要求的精度，就完成了方程組的求解。求解方程組（ 11 ）還有許多其他辦法，如“最速下降法”，它是利用方程組（ 11 ）構成所謂模函數(shù)xf1 x, y 2f 2 x, y 2 ，通過求模函數(shù)極小值的方法得到方程組

25、的數(shù)值解，諸如此類在此不再一一列舉。三、求方程組數(shù)值解的指令fsolve 是用最小二乘法求解非線性方程組F (X） 0 的指令，變量 X 可以是向量或矩陣，方程組可以由代數(shù)方程或者超越方程構成。它的使用格式為：fsolve ( fun , X0 , OPTIONS ) 參數(shù) fun是編輯并存盤的M 函數(shù)文件的名稱，可以用代替單引號對它進行標識。M 函數(shù)文件主要內(nèi)容是方程F ( X ) = 0中的函數(shù)F (X) ，即方文檔實用標準程左邊的函數(shù)。 參數(shù) X0 是向量或矩陣， 為探索方程組解的起始點。 求解將從 X0 出發(fā)，逐漸趨向，最終得到滿足精度要求、最接近X0 的近似根 X*:F(X*） 0

26、。由于 X0 是向量或矩陣，無法用畫圖方法進行估計，實際問題中常常是根據(jù)專業(yè)知識、物理意義等進行估計。 該指令輸出一個與 X0 同維的向量或矩陣，為方程組的近似數(shù)值解。 參數(shù) OPTIONS 為設置選項，用它可以設置過程顯示與否、誤差、算法 ，具體內(nèi)容可用help 查閱。通常可以省略該項內(nèi)容。x2y2z710x例 求方程組xy22z在 x0 1 ，y 0 1 ，z0 1 附近的數(shù)值解。x2y 2z23y解（l）在文本編輯調(diào)試窗中編輯M 函數(shù)文件。首先將方程組變換成 F ( X )x210xy2z70= 0 的形式， xy22z0x, y, z 看成向量 X 的三個分量。x2y23 yz20fu

27、nction ms=li4_7(X)ms(1) = X(1).2 - 10* X(1) + X(2).2 + X(3) +7;ms(2) = X(1).*X(2).2 - 2* X(3);ms(3) = X(1).2 + X(2).2 - 3* X(2) + X(3).2;這里， X (l ） x , X (2 ） y , X (3 ） z，輸出參量ms 也有三個分量用“ li4_7 ”為 M 函數(shù)文件名存盤。(2）在指令窗中鍵入：fsolve ( li4_7,1 1 1)文檔實用標準若鍵入：x = fsolve (li4_7,1 1 1, optimset(Display, iter)則得出求解過程。該方程也可以用MATLAB的符號指令 solve 求解，但結果非常冗長。3xcos yz0.5例 求解方程組 2x281 y0.1 2sin z 1.06 0，在 x0 0.1 ，y 0 0.1 和e xy20z1013z0 - 0.1 附近的數(shù)值解。解：首先將方程組變換成 fj (x, y, z） f (X ） 0 ( j 1 , 2 , 3 ）的形式，設 X 為一個三維向量，令 X (l ） x , X (2 ） y , X (3 ） z，則三維向量 yy3 f (X ） fj (x, y, z ），然后編程計算。(l) 在文本編輯調(diào)試窗中編輯