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

第一章方程（組）模型2023/2/41本章學(xué)習(xí)目的：本章我們就是要學(xué)習(xí)求解線性方程組、非線性方程（組）的方法，以及利用數(shù)學(xué)軟件和計算機(jī)對方程和方程組進(jìn)行求解。考慮求方程f(x)=0的解，我們通常采用這樣的幾種方法：因式分解法、圖形放大法、數(shù)值迭代逼近法。2023/2/42本章學(xué)習(xí)目的：復(fù)習(xí)求解方程的基本原理和方法，掌握解方程的圖形放大法和迭代算法；熟練掌握用MATLAB軟件的函數(shù)來求解方程和方程組；通過范例展現(xiàn)求解實(shí)際問題的初步建模過程和MATLAB程序設(shè)計。2023/2/431.1圖形放大法：由于計算機(jī)的廣泛應(yīng)用，可以非常方便地作出函數(shù)f(x)的圖形（曲線），找出曲線與x軸的交點(diǎn)的橫坐標(biāo)值，就可求出f(x)＝0的近似根。這些值盡管不精確，但是直觀，方程有多少個根、在什么范圍，一目了然。并且可以借助于計算機(jī)使用圖形局部放大功能，將根定位得更加準(zhǔn)確。2023/2/44用圖形放大法求解方程f(x)＝0的步驟：建立坐標(biāo)系，作曲線f(x)；觀察f(x)與x軸的交點(diǎn)；將其中的一個交點(diǎn)進(jìn)行局部放大；該交點(diǎn)的橫坐標(biāo)值就是方程的一個根；對所有的交點(diǎn)進(jìn)行相同的處理，就得到方程的所有解。2023/2/45例1.1求方程所有的根及大致分布范圍，欲尋求其中的一個實(shí)根，并且達(dá)到一定的精度。畫出的圖形；程序如下：x=-6:0.01:6;y=x.^5+2*x.^2+4;plot(x,y)gridon;%畫坐標(biāo)格2023/2/46我們可以看出方程在－2～＋2范圍有一個實(shí)根。2023/2/47逐次縮小范圍得到較精確的根。最后我們可以看出這個實(shí)根的值在-1.56~-1.54之間。2023/2/481.2簡單迭代法迭代算法步驟：對方程f(x)=0求解對方程經(jīng)過簡單變形得到（不是唯一的），x被稱之為不動點(diǎn)；設(shè)置為迭代初值，迭代過程為，n=0,1,2……當(dāng)兩次迭代結(jié)果之差小于某個設(shè)定的誤差值時，我們認(rèn)為迭代結(jié)果是收斂的，可得到結(jié)果的近似值。2023/2/49例1.2求方程的非負(fù)實(shí)根。解：由于函數(shù)連續(xù)，并且在x=0和x=1處函數(shù)值符號相反，可以判斷函數(shù)在區(qū)間（0，1）必有零點(diǎn)，即方程在（0，1）內(nèi)必然存在根。先將函數(shù)變形為；設(shè)置迭代初值為0，編程進(jìn)行迭代。2023/2/410n=1;x=0;y=exp(x)/3;ys=vpa(y,10);%ys取y的10位有效數(shù)字z=abs(y-x);whilez>10^(-5)

x=y;y=exp(x)/3;ys=vpa(y,10);z=abs(y-x);n=n+1;endn,y,ys

n=21y=0.6190ys=.6190471917從該結(jié)果可以看出，迭代21次后兩次迭代的結(jié)果誤差值滿足小于的條件，結(jié)果收斂，迭代結(jié)果為0.6190，若保留小數(shù)點(diǎn)后10位有效數(shù)字則結(jié)果為0.6190471917。2023/2/411例1.3用迭代方法求解方程解：（1）對方程變形為，有不同的形式，比如：；(a)

；(b)

；……(c)2023/2/412（2）設(shè)定初始值為1，編程迭代求解x=1;y=1;z=1;fork=1:25x=x^3-x^2-1;y=(y^2+y+1)^(1/3);z=1+1/z+1/z^2;endx,y,z

x=-Infy=1.8393z=1.8393在程序中，函數(shù)x,y,z分別對應(yīng)方程（a）（b）（c），從結(jié)果可以看出方程（a）不收斂，結(jié)果趨于負(fù)無窮大，方程(b)(c)收斂，結(jié)果為1.8393。而且，還可證明（b）比（c）收斂速度快。2023/2/4131.3方程組的求解方法1.3.1線性方程組的求解對于線性方程組

可以寫成矩陣的形式我們在線性代數(shù)中已經(jīng)學(xué)習(xí)了線性方程組的求解方法。

2023/2/414由線性代數(shù)的知識可知，線性方程組的解可能出現(xiàn)三種情形：無解、有唯一解、有無窮多組解。這主要取決于系數(shù)矩陣A的秩與增廣矩陣（A︱b）的秩是否相等、秩與變量個數(shù)是否相等，具體地：若R(A)≠R(A︱b）,則無解；若R(A)＝R(A︱b）＝n（n為變量個數(shù)），則有唯一一組解；若R(A)＝R(A︱b）<n，則有無窮多組解。求矩陣A的秩可以很方便的用MATLAB的rank(A)函數(shù)求得。2023/2/415求解線性方程組的方法大致可以劃分為兩類：直接消去法、迭代數(shù)值解法。直接消去法在線性代數(shù)中已經(jīng)學(xué)過，這里不再贅述。線性方程組可以看成是非線性方程組的特例，其迭代數(shù)值解法相同，在非線性方程組的迭代解法中介紹方法。2023/2/4161．3．2非線性方程組的迭代解法非線性方程組的一般形式為可以改寫為等價的方程組用這個方程組進(jìn)行迭代求得精確解。2023/2/417例1．4求解方程組解：對方程組進(jìn)行變形，構(gòu)造如下的迭代函數(shù)：或2023/2/418求解方程組迭代產(chǎn)生的序列是數(shù)組：迭代程序如下：（初始值是（2,3），迭代次數(shù)為20次）x=[2,3];y=[2,3];fork=1:20a=0.1*x(1)^2+0.1*x(2)^2+0.8;x(2)=0.1*x(1)*x(2)^2+0.1*x(1)+0.8;x(1)=a;b=(10*y(2)-8)/(y(2)^2+1);y(2)=sqrt(10*y(1)-8-y(1)^2);y(1)=b;endx,y x=1.00001.0000y=2.19343.02052023/2/4191.3MATLAB軟件直接求解法1.3.1任意函數(shù)方程與線性方程組可用命令solve()求解2023/2/4201.3MATLAB軟件直接求解法solve()語句的用法：1.單變量方程：f(x)=0（1）符號解：例1．3求解方程：解：輸入：x=solve('a*x^2+b*x+c')

輸出為：x=[1/2/a*(-b+(b^2-4*a*c)^(1/2))][1/2/a*(-b-(b^2-4*a*c)^(1/2))]2023/2/421（2）數(shù)字解：如果不能求得精確的符號解，可以計算可變精度的數(shù)值解。例1.4

解方程：解：s=solve('x^3-2*x^2=x-1')

s=1/6*(28+84*i*3^(1/2))^(1/3)+14/3/(28+84*i*3^(1/2))^(1/3)+2/3-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))-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))2023/2/422double(s)

ans=2.2470-0.8019+0.0000i0.5550-0.0000ivpa(s,10)

ans=2.246979604+.1e-9*i-.8019377358-.1866025404e-9*i.5549581322-.1339745960e-10*i該方程無實(shí)根。2023/2/423（3）無窮解例1.5求解方程：解：輸入：solve('tan(x)-sin(x)=0')輸出為：ans=0 該方程有無窮多個解，不能給出全部解，這里只得到其中的一個。2023/2/4242．多變量方程組：例1.6

解方程組解：輸入：[x,y]=solve('x^2*y^2','x-y/2-b')

x=00bby=-2*b-2*b00v=[0,-2*b][0,-2*b][b,0][b,0]v=[x,y]2023/2/4251.3.2非線性方程組非線性方程組仍然可以用solve()求解（如上例1.6），一般給出的是數(shù)值解。也可以用fsolve()函數(shù)求解，格式是：這里

為變量的初始值；fun為m文件的文件名。m文件如下所示：2023/2/426例1.7

求解方程組：解：這是一個非線性方程組。（1）首先建立關(guān)于該方程組的M文件nxxffunctioneq=nxxf(x)globalnumber;number=number+1;eq(1)=sin(x(1))+x(2)^2+log(x(3))-7;eq(2)=3*x(1)+2^x(2)-x(3)^3+1;eq(3)=x(1)+x(2)+x(3)-5;2023/2/427（2）執(zhí)行以下程序globalnumber;number=0;x=fsolve('nxxf',[1,1,1])number

x=0.59912.39592.0050number=29這里迭代步驟為29次。2023/2/4281.3.3多項式方程求解多項式方程除了可以用上面講述的solve()函數(shù)外，還可以直接用求解多項式方程的函數(shù)roots(p)。

對于多項式方程用：求解該方法可以求出方程的全部根（包含重根）。2023/2/429例1．8

求解多項式方程：解：p=[1,1,0,0,0,0,0,0,0,1] ;roots(p)運(yùn)行結(jié)果：

ans=-1.2131-0.9017+0.5753i-0.9017-0.5753i-0.2694+0.9406i-0.2694-0.9406i0.4168+0.8419i0.4168-0.8419i0.8608+0.3344i0.8608-0.3344i2023/2/4301.3.4線性方程組線性方程組除了可以使用solve()求解外，還可以使用其他的MATLAB命令。 將線性方程組寫成矩陣形式：AX=b，就可以考慮用以下幾種形式之一求解。linsolve(A,b)；6.5版無此函數(shù)sym(A)\sym(b)；A\b；inv(A)*b；inv(A)表示A的逆矩陣，這里A必須為方陣且可逆；pinv(A)*b；pinv(A)表示A的逆矩陣，A可以為任意矩陣。2023/2/431例1.9

求解線性方程組：AX＝b(1)A=[410;1-15;22-3];b=[6;14;-3];x=A\b

x=1232023/2/4321.4案例詳解：山崖高度某人站在山崖頂且身上帶著一只具有跑表功能的手表，出于好奇心想用扔下一塊石頭聽回聲的方法來估計山崖的高度，假定他準(zhǔn)確地測定時間t=4s，那么應(yīng)該怎樣來推算山崖的高度呢？模型一：假定空氣阻力不計，可以直接利用自由落體運(yùn)動的公式來計算。（1）取g=9.8m/s2，則可求得。這樣計算非常簡單，但略去了一些影響因素，誤差較大。2023/2/433模型二：除去地球吸引力外，對石塊下落影響最大的當(dāng)屬空氣阻力。根據(jù)流體力學(xué)知識，可設(shè)空氣阻力正比于石塊下落的速度，阻力系數(shù)K為常數(shù)，因而，由牛頓第二定律可得：令，解得代入初始條件：，得對v再積分一次，得2023/2/434代入初始條件：，得得到計算山崖高度的公式：（2）若設(shè)k=0.05，將t=4s代入上式，則可求得2023/2/435模型三：進(jìn)一步考慮，聽到回聲按下跑表，t=4s包含了反應(yīng)時間，不妨假設(shè)平均反應(yīng)時間為0.1s，則扣除反應(yīng)時間后下落時間為3.9s，代入式（2）求得。2023/2/436模型四：再深入考慮回聲傳回來所需要的時間，為此，令石塊下落的的真正時間為,聲音傳回來的時間記為，我們知道聲音在空氣中的傳播速度是340m/s,因此得到一個方程組2023/2/437利用MATLAB的fsolve()函數(shù)求解，先建立m文件mx4.mfunctionf=mx4(x)globalnumber;number=number+1;g=9.8;k=0.05;

f(1)=x(1)-g*x(2)/k-g*exp(-k*x(2))/k^2+g/k^2;f(2)=x(1)-340*x(3);f(3)=x(2)+x(3)-3.9;2023/2/438再執(zhí)行以下程序：globalnumber;number=0;x=fsolve('mx4',[0,0,0])number

得到結(jié)果：x=63.56163.71310.1869number=125即迭代125次后求出，這樣計算的過程實(shí)際上就是一個數(shù)學(xué)建模并求解的過程，通過循序漸進(jìn)的計算得到了越來越精確的解。2023/2/439范例：波音公司飛機(jī)最佳定價策略問題全球最大的飛機(jī)制造商——波音公司自1955年推出的波音707開始，成功地開發(fā)了一系列的噴氣式客機(jī)。問題：討論該公司對一種新型客機(jī)最優(yōu)定價策略的數(shù)學(xué)模型。問題分析定價策略涉及到諸多因素，這里考慮以下主要因素：價格、競爭對手的行為、出售客機(jī)的數(shù)量、波音公司的客機(jī)制造量、制造成本、波音公司的市場占有率等等因素。2023/2/440假設(shè)及模型價格記為p，根據(jù)實(shí)際情況，對于民航飛機(jī)制造商，能夠與波音公司抗衡的競爭對手只有一個，因此他們可以在價格上達(dá)成一致，具體假設(shè)如下：1.型號：為了研究方便，假設(shè)只有一種型號飛機(jī)；2.銷售量：其銷售量只受飛機(jī)價格p的影響。預(yù)測以此價格出售，該型號飛機(jī)全球銷售量為N。N應(yīng)該受到諸多因素的影響，假設(shè)其中價格是最主要的因素。根據(jù)市場歷史的銷售規(guī)律和需求曲線，假設(shè)該公司銷售部門預(yù)測得到：2023/2/4413.市場占有率：既然在價格上達(dá)成一致，即價格的變化是同步的，因此，不同定價不會影響波音公司的市場占有率，因此市場占有率是常數(shù)，記為h。4.制造數(shù)量：假設(shè)制造量等于銷售量，記為x。既然可以預(yù)測該型號飛機(jī)全球銷售量，結(jié)合波音公司的市場占有率，可以得到2023/2/4425.制造成本：根據(jù)波音產(chǎn)品分析部門的估計，制造成本為：6.利潤：假設(shè)利潤等于銷售收入去掉成本，并且公司的最優(yōu)策略原則為利潤R(p)最大。2023/2/443由以上簡化的分析及假設(shè)得到波音公司飛機(jī)最佳定價策略的數(shù)學(xué)模型如下：其中：

p，x，N≥0。2023/2/444模型求解我們采用圖形放大的方法求解。具體用Matlab作出目標(biāo)函數(shù)曲線圖，得到一個直觀的印象：最優(yōu)定價策略下價格p大致在6到7之間；再用圖形放大方法，進(jìn)一步估計出p≈6.2859，R=1780.83362023/2/445Matlab程序如下(作函數(shù)曲線圖的基本程序)h=0.5;a=6.285;b=6.287;n=80;d=(b-a)/n;fori=1:n+1pr(i)=a+(i-1)*d;p=pr(i);n=-78*p^2+655*p+125;x=h*n;r=p*x;c=50+