matlab教學(xué)第一章方程(組)模型

1、第一章第一章 方程（組）模型方程（組）模型2022-3-161本章學(xué)習(xí)目的：本章學(xué)習(xí)目的：l本章我們就是要學(xué)習(xí)求解線性方程組、非線性方程（組）的方法，以及利用數(shù)學(xué)軟件和計算機對方程和方程組進行求解。l考慮求方程f(x)=0的解，我們通常采用這樣的幾種方法：因式分解法、圖形放大法、數(shù)值迭代逼近法。2022-3-162本章學(xué)習(xí)目的：本章學(xué)習(xí)目的：l復(fù)習(xí)求解方程的基本原理和方法，掌握解方程的圖形放大法和迭代算法；l熟練掌握用MATLAB軟件的函數(shù)來求解方程和方程組；l通過范例展現(xiàn)求解實際問題的初步建模過程和MATLAB程序設(shè)計。2022-3-1631.1 圖形放大法：l由于計算機的廣泛應(yīng)用，可以非常

2、方便地作出函數(shù)函數(shù)f(x)的圖形的圖形（曲線），找出曲線與x軸的交點的橫坐標(biāo)值，就可求出f(x)0的近似根。l這些值盡管不精確，但是直觀，方程有多少個根、在什么范圍，一目了然。并且可以借助于計算機使用圖形局部放大圖形局部放大功能，將根定位定位得更加準(zhǔn)確準(zhǔn)確。2022-3-164用圖形放大法求解方程用圖形放大法求解方程f(x)0的步驟：的步驟：l建立坐標(biāo)系，作曲線f(x)；l觀察f(x)與x軸的交點；l將其中的一個交點進行局部放大；l該交點的橫坐標(biāo)值就是方程的一個根；l對所有的交點進行相同的處理，就得到方程的所有解。2022-3-165例例1.1 求方程求方程 所有的根及大致分所有的根及大致分布

3、范圍，欲尋求其中的一個實根，并且達到一布范圍，欲尋求其中的一個實根，并且達到一定的精度。定的精度。l畫出 的圖形；l程序如下：lx=-6:0.01:6;ly=x.5+2*x.2+4;lplot(x,y)lgrid on; %畫坐標(biāo)格52240 xx42)(25xxxf2022-3-166l我們可以看出方程在22范圍有一個實根。-6-4-20246-8000-6000-4000-2000020004000600080002022-3-167逐次縮小范圍得到較精確的根。逐次縮小范圍得到較精確的根。-2-1.5-1-0.500.511.52-20-1001020304050-1.6-1.58-1.5

4、6-1.54-1.52-1.5-1.48-1.46-1.5-1-0.500.51最后我們可以看出這個實根的值在-1.56-1.54之間。2022-3-1681.2 簡單迭代法簡單迭代法l迭代算法步驟：迭代算法步驟： l對方程f(x)=0求解l對方程經(jīng)過簡單變形得到 （不是唯一的），x 被稱之為不動點；l設(shè)置為迭代初值 ，迭代過程為 ，n=0,1,2當(dāng)兩次迭代結(jié)果之差小于某個設(shè)定的誤差值時，我們認(rèn)為迭代結(jié)果是收斂的，可得到結(jié)果的近似值 。)(xx)(1nnxx1nxx2022-3-169nx例1.2 求方程求方程 的非負(fù)實根。的非負(fù)實根。l解：由于函數(shù) 連續(xù)，并且在x=0和x=1處函數(shù)值符號相反

5、，可以判斷函數(shù)在區(qū)間（0，1）必有零點，即方程在（0，1）內(nèi)必然存在根。l先將函數(shù)變形為 ；l設(shè)置迭代初值為0，編程進行迭代。03xexxex 33/xex 2022-3-1610ln=1;lx=0;ly=exp(x)/3;lys=vpa(y,10); %ys取取y的的10位有效數(shù)字位有效數(shù)字lz=abs(y-x);lwhile z10(-5)l x=y;l y=exp(x)/3;l ys=vpa(y,10);l z=abs(y-x);l n=n+1;lendln,y,ys n = 21y = 0.6190ys =.6190471917 從該結(jié)果可以看出，迭代從該結(jié)果可以看出，迭代21次后次后

6、兩次迭代的結(jié)果誤差值滿足小于兩次迭代的結(jié)果誤差值滿足小于的條件，結(jié)果收斂，迭代結(jié)果為的條件，結(jié)果收斂，迭代結(jié)果為0.6190，若保留小數(shù)點后，若保留小數(shù)點后10位有位有效數(shù)字則結(jié)果為效數(shù)字則結(jié)果為0.6190471917。2022-3-1611例例1.3 用迭代方法求解方程用迭代方法求解方程解：（1）對方程變形為 ，有不同的形式，比如： ； (a) ； (b) ； (c)0123xxx)(xx123xxx321xxx2111xxx2022-3-1612（2）設(shè)定初始值為）設(shè)定初始值為1，編程迭代求解，編程迭代求解lx=1;y=1;z=1;lfor k=1:25l x=x3-x2-1;l y=(

7、y2+y+1)(1/3);l z=1+1/z+1/z2;lendlx,y,z x = -Infy = 1.8393z = 1.8393 在程序中，函數(shù)在程序中，函數(shù)x,y,z分別對應(yīng)方程（分別對應(yīng)方程（a）（b ）（）（c），從結(jié)果可以看出方程（），從結(jié)果可以看出方程（ a）不收斂，結(jié)果趨于負(fù)無窮大，方程不收斂，結(jié)果趨于負(fù)無窮大，方程(b)(c)收斂，結(jié)果為收斂，結(jié)果為1.8393。而且，還可證明。而且，還可證明（b）比（）比（c）收斂速度快。）收斂速度快。2022-3-16131.3 方程組的求解方法方程組的求解方法l1.3.1 線性方程組的求解線性方程組的求解l對于線性方程組對于線性方程組

8、 l可以寫成矩陣的形式可以寫成矩陣的形式l我們在線性代數(shù)中已經(jīng)學(xué)習(xí)了線性方程組的求我們在線性代數(shù)中已經(jīng)學(xué)習(xí)了線性方程組的求解方法。解方法。 11 111,1 1nnmmnnma xa xba xa xb. bAX 2022-3-1614l由線性代數(shù)的知識可知，線性方程組的解可能由線性代數(shù)的知識可知，線性方程組的解可能出現(xiàn)三種情形：無解、有唯一解、有無窮多組出現(xiàn)三種情形：無解、有唯一解、有無窮多組解。這主要取決于系數(shù)矩陣解。這主要取決于系數(shù)矩陣A 的秩與增廣矩陣的秩與增廣矩陣（Ab）的秩是否相等、秩與變量個數(shù)是否）的秩是否相等、秩與變量個數(shù)是否相等，具體地：相等，具體地：l若若R(A)R(Ab）

9、,則無解；則無解；l若若R(A)R(Ab）n（n為變量個數(shù)），則為變量個數(shù)），則有唯一一組解；有唯一一組解；l若若R(A)R(Ab） n，則有無窮多組解。，則有無窮多組解。l求矩陣求矩陣A的秩可以很方便的用的秩可以很方便的用MATLAB的的rank(A)函數(shù)求得。函數(shù)求得。2022-3-1615l求解線性方程組的方法大致可以劃分為兩類：直接消去法直接消去法、迭代數(shù)值解法迭代數(shù)值解法。l直接消去法在線性代數(shù)中已經(jīng)學(xué)過，這里不再贅述。l線性方程組可以看成是非線性方程組的特例，其迭代數(shù)值解法相同，在非線性方程組的迭代解法中介紹方法。2022-3-1616132 非線性方程組的迭代解法非線性方程組的迭

10、代解法l非線性方程組的一般形式為 l可以改寫為等價的方程組 l用這個方程組進行迭代求得精確解。 0),(, 0),(111nnnxxfxxf),(),(1111nnnnxxgxxxgx2022-3-1617l例14 求解方程組l解：對方程組進行變形，構(gòu)造如下的迭代函數(shù)：l或081008102122122121xxxxxxx8 . 01 . 01 . 08 . 01 . 01 . 01221222211xxxxxxx211222218101810 xxxxxx2022-3-1618求解方程組迭代產(chǎn)生的序列是數(shù)組：求解方程組迭代產(chǎn)生的序列是數(shù)組：l迭代程序如下：（初始值是（2,3），迭代次數(shù)為20

11、次）lx=2,3;y=2,3;lfor k=1:20la=0.1*x(1)2+0.1*x(2)2+0.8;lx(2)=0.1*x(1)*x(2)2+0.1*x(1)+0.8;lx(1)=a;lb=(10*y(2)-8)/(y(2)2+1);ly(2)=sqrt(10*y(1)-8-y(1)2);ly(1)= b;lendlx,y),( ,),(),2122211211nnxxxxxx（x = 1.0000 1.0000y = 2.1934 3.0205 2022-3-16191.3 MATLAB軟件直接求解法軟件直接求解法l1.3.1 任意函數(shù)方程與線性方程組可用命令任意函數(shù)方程與線性方程組可

12、用命令solve( )求解求解)( ,)( ,)(),(, 0)(, 0)(2111xfxfxfsolvexxxxfxfnnn2022-3-16201.3 MATLAB軟件直接求解法軟件直接求解法lsolve( )語句的用法：語句的用法：l1. 單變量方程：單變量方程：f(x)=0 l（1）符號解：）符號解：l例例13 求解方程：求解方程：l解：輸入：x=solve(a*x2+b*x+c) l輸出為：lx =l 1/2/a*(-b+(b2-4*a*c)(1/2)l 1/2/a*(-b-(b2-4*a*c)(1/2) 02cbxax2022-3-1621l（2）數(shù)字解：l如果不能求得精確的符號解

13、，可以計算可變精度的數(shù)值解。l例例1.4 解方程： l解：s=solve(x3-2*x2=x-1) ls = 1/6*(28+84*i*3(1/2)(1/3)+14/3/(28+84*i*3(1/2)(1/3)+2/3l -1/12*(28+84*i*3(1/2)(1/3)-7/3/(28+84*i*3(1/2)(1/3)+2/3+1/2*i*3(1/2)*(1/6*(28+84*i*3(1/2)(1/3)-14/3/(28+84*i*3(1/2)(1/3)l -1/12*(28+84*i*3(1/2)(1/3)-7/3/(28+84*i*3(1/2)(1/3)+2/3-1/2*i*3(1/2

14、)*(1/6*(28+84*i*3(1/2)(1/3)-14/3/(28+84*i*3(1/2)(1/3) 1223xxx2022-3-1622ldouble(s) lans =l 2.2470 l -0.8019 + 0.0000il 0.5550 - 0.0000i lvpa(s,10) lans =l2.246979604+.1e-9*il -.8019377358-.1866025404e-9*il .5549581322-.1339745960e-10*i l該方程無實根。2022-3-1623l（3）無窮解l例1.5 求解方程：l解：輸入：solve(tan(x)-sin(x)=0

15、) l輸出為：lans =l0 l該方程有無窮多個解，不能給出全部解，這里只得到其中的一個。2022-3-16242多變量方程組：多變量方程組：l例例1.6 解方程組 l解：輸入：lx,y=solve(x2*y2,x-y/2-b) 0)(, 0)(1xfxfnbyxyx2022x = 0 0 b by = -2*b -2*b 0 0 v = 0, -2*b 0, -2*b b, 0 b, 0 v=x,y2022-3-16251.3.2 非線性方程組非線性方程組l非線性方程組仍然可以用solve( )求解（如上例1.6），一般給出的是數(shù)值解。也可以用fsolve( )函數(shù)求解，格式是：l這里 為

16、變量的初始值；fun為m文件的文件名。m文件如下所示：110( )0,( )0,( ,)(,)nnf xfxxxxxfsolve fun x0 x1( )(1)( );( )( )nfunctionffun xff xf nfx2022-3-1626l例例1.7 求解方程組：l解：這是一個非線性方程組。l（1）首先建立關(guān)于該方程組的M文件nxxflfunction eq=nxxf(x)lglobal number;lnumber=number+1;leq(1)=sin(x(1)+x(2)2+log(x(3)-7;leq(2)=3*x(1)+2x(2)-x(3)3+1;leq(3)=x(1)+x

17、(2)+x(3)-5;05012307lnsin32zyxzxzyxy2022-3-1627l（2）執(zhí)行以下程序lglobal number;lnumber=0;lx=fsolve(nxxf,1,1,1)lnumber lx =l 0.5991 2.3959 2.0050lnumber =l 29 l這里迭代步驟為29次。2022-3-16281.3.3 多項式方程多項式方程l求解多項式方程除了可以用上面講述的solve( ) 函數(shù)外，還可以直接用求解多項式方程的函數(shù)roots(p)。 l對于多項式方程l用： 求解l該方法可以求出方程的全部根（包含重根）。00111axaxaxammmm)(;

18、,011prootsaaaapmm2022-3-1629l例例18 求解多項式方程：l解：p=1,1,0,0,0,0,0,0,0,1;lroots(p)l運行結(jié)果：運行結(jié)果： lans =l -1.2131 l -0.9017 + 0.5753il -0.9017 - 0.5753il -0.2694 + 0.9406il -0.2694 - 0.9406il 0.4168 + 0.8419il 0.4168 - 0.8419il 0.8608 + 0.3344il 0.8608 - 0.3344i 0189 xx2022-3-16301.3.4 線性方程組線性方程組l線性方程組除了可以使用s

19、olve( ) 求解外，還可以使用其他的MATLAB命令。l將線性方程組寫成矩陣形式：AX=b，就可以考慮用以下幾種形式之一求解。llinsolve(A,b)；6.5版無此函數(shù)lsym(A)sym(b)；lAb；linv(A)*b； inv(A)表示A的逆矩陣，這里A必須為方陣且可逆；lpinv(A)*b； pinv(A) 表示A的逆矩陣，A可以為任意矩陣。2022-3-1631l例例1.9 求解線性方程組：AXbl(1)lA=4 1 0 ;1 -1 5;2 2 -3;b=6;14;-3;lx=Ab lx =l 1l 2l 3 2022-3-16321.4 案例詳解：山崖高度案例詳解：山崖高度

20、l某人站在山崖頂且身上帶著一只具有跑表功能的手表，出于好奇心想用扔下一塊石頭聽回聲的方法來估計山崖的高度，假定他準(zhǔn)確地測定時間t=4s，那么應(yīng)該怎樣來推算山崖的高度呢？l模型一：模型一：假定空氣阻力不計，可以直接利用自由落體運動的公式來計算。l （1） l取g=9.8m/s2，則可求得 。這樣計算非常簡單，但略去了一些影響因素，誤差較大。221gth mh4 .782022-3-1633l模型二：模型二：l除去地球吸引力外，對石塊下落影響最大的當(dāng)屬空氣阻力。根據(jù)流體力學(xué)知識，可設(shè)空氣阻力正比于石塊下落的速度，阻力系數(shù)K為常數(shù)，因而，由牛頓第二定律可得：l令 ，解得l代入初始條件： ，得l對v再

21、積分一次，得KvmgdtdvmFmKk/kgecvkt10)0(vkgc/122cekgtkghkt2022-3-1634l代入初始條件： ，l得l得到計算山崖高度的公式：l （2）l若設(shè)k=0.05，將t=4s代入上式，則可求得0)0(h22/kgc22kgekgtkghktmh4 .732022-3-1635l模型三：模型三：l進一步考慮，聽到回聲按下跑表，進一步考慮，聽到回聲按下跑表，t=4s包含了包含了反應(yīng)時間，不妨假設(shè)平均反應(yīng)時間為反應(yīng)時間，不妨假設(shè)平均反應(yīng)時間為0.1s，則，則扣除反應(yīng)時間后扣除反應(yīng)時間后下下落時間為落時間為3.9s，代入式（，代入式（2）求得求得 。mh9 .69

22、2022-3-1636l模型四：模型四：l再深入考慮回聲傳回來所需要的時間，為此，再深入考慮回聲傳回來所需要的時間，為此，令石塊下落的的真正時間為令石塊下落的的真正時間為 ,聲音傳回來的聲音傳回來的時間記為時間記為 ，我們知道聲音在空氣中的傳播速，我們知道聲音在空氣中的傳播速度是度是340m/s,因此得到一個方程組因此得到一個方程組1t2t9 . 33402122211ttthkgekgtkghkt2022-3-1637l利用MATLAB的fsolve()函數(shù)求解，先建立m文件mx4.mlfunction f=mx4(x)lglobal number;lnumber=number+1;lg=9

23、.8;lk=0.05; lf(1)=x(1)-g*x(2)/k-g*exp(-k*x(2)/k2+g/k2;lf(2)=x(1)-340*x(3);lf(3)=x(2)+x(3)-3.9;2022-3-1638l再執(zhí)行以下程序：lglobal number;lnumber=0;lx=fsolve(mx4,0,0,0)lnumber l得到結(jié)果：lx = 63.5616 3.7131 0.1869lnumber = 125l即迭代125次后求出， l這樣計算的過程實際上就是一個數(shù)學(xué)建模并求解的過程，通過循序漸進的計算得到了越來越精確的解。ststmh1869. 0,7131. 3,5616.63

24、212022-3-1639范例：波音公司飛機最佳定價策略范例：波音公司飛機最佳定價策略l問題問題l全球最大的飛機制造商波音公司自1955年推出的波音707開始，成功地開發(fā)了一系列的噴氣式客機。問題：討論該公司對一種新型客機最優(yōu)定價策略的數(shù)學(xué)模型。l問題分析問題分析l定價策略涉及到諸多因素，這里考慮以下主要因素： l價格、競爭對手的行為、出售客機的數(shù)量、波音公司的客機制造量、制造成本、波音公司的市場占有率等等因素。2022-3-1640l假設(shè)及模型假設(shè)及模型l價格記為p，根據(jù)實際情況，對于民航飛機制造商，能夠與波音公司抗衡的競爭對手只有一個，因此他們可以在價格上達成一致，具體假設(shè)如下： 1. 型

25、號：為了研究方便，假設(shè)只有一種型號飛機；2. 銷售量：其銷售量只受飛機價格p的影響。預(yù)測以此價格出售，該型號飛機全球銷售量為N。N應(yīng)該受到諸多因素的影響，假設(shè)其中價格是最主要的因素。根據(jù)市場歷史的銷售規(guī)律和需求曲線，假設(shè)該公司銷售部門預(yù)測得到：12565578)2pppNN（2022-3-16413. 市場占有率：既然在價格上達成一致，即價格的變化是同步的，因此，不同定價不會影響波音公司的市場占有率，因此市場占有率是常數(shù)市場占有率是常數(shù)，記為h。4. 制造數(shù)量：假設(shè)制造量等于銷售量，記為x。既然可以預(yù)測該型號飛機全球銷售量，結(jié)合波音公司的市場占有率，可以得到)(*pNhx 2022-3-164

26、25. 制造成本：根據(jù)波音產(chǎn)品分析部門的估計，制造成本為：6. 利潤：假設(shè)利潤等于銷售收入去掉成本，并且公司的最優(yōu)策略原則為利潤R(p)最大。4385 . 150)(xxxC)()(xCpxpR2022-3-1643l由以上簡化的分析及假設(shè)得到波音公司飛機最佳定價策略的數(shù)學(xué)模型如下：l其中：l p，x，N0。 )()(xCpxpMaxR12565578)2pppNN（)(*pNhx 4385 . 150)(xxxC2022-3-1644模型求解模型求解l我們采用圖形放大圖形放大的方法求解。具體用Matlab 作出目標(biāo)函數(shù)曲線圖，得到一個直觀的印象：最優(yōu)定價策略下價格p大致在6到7之間；再用圖形放大方法，進一步估計出 lp 6.2859，R=1780.8336 6.2856.28526.28546.28566.28586.2866.28626.28646.28666.28686.2871780.83321780.83331780.83331780.83341780.83341780.83351780.83351780.83361780.8336價 格 p利 潤 曲 線 R(p)2022-3-1645lMatlab 程序如下(作函數(shù)曲線圖的基本程序)lh=0.5;a=6.2