部分的講函數(shù)的等高線

1、第二講 函數(shù)的等高線、梯度線及有關(guān)的作圖問題鯊魚襲擊目標(biāo)的前進(jìn)途徑等高線和梯度線有廣泛的實(shí)際應(yīng)用，例如在地理學(xué)中繪制地貌圖，在氣象學(xué)中繪制氣象圖等等.本實(shí)驗(yàn)通過鯊魚襲擊目標(biāo)這一例子介紹二元函數(shù)的等高線和梯度線的繪制，最后介紹用等高線來做一元隱函數(shù)的圖形及微分方程的積分曲線.2.1 等高線的繪制二元函數(shù)在空間表示的是一張曲面，這個(gè)曲面與平面的交線在面上的投影曲線稱為函數(shù)的一條等高線，我們可以用Matlab作出等高線的圖形.等高線的作圖命令為“Contour”，最基本格式為Contour二元函數(shù)，自變量1，自變量1最小值，自變量1最大值，自變量2，自變量2最小值，自變量2最大值例1 作出在區(qū)間上的

2、等高線.解 X,Y = meshgrid(-2:.2:2,-2:.2:3);Z = X.*exp(-X.2-Y.2);C,h = contour(X,Y,Z);set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)colormap cool運(yùn)行后見圖（2.1）.2.2 矢量場圖矢量場圖（又稱速度圖）是指由指令quiver實(shí)現(xiàn)的.它主要用于描寫函數(shù)在點(diǎn)的梯度大小和方向。其一般的調(diào)用格式為： quiver(X,Y,DZX,DZY)例2 作出函數(shù)的等高線和矢量場.解 X,Y = me

3、shgrid(-2:.2:2,-1:.2:2);Z = X.*exp(-X.2 - Y.2);請預(yù)覽后下載！DX,DY = gradient(Z,.2,.2);% 求二元函數(shù)矩陣Z的梯度指令，0.2 為x、y方向上的計(jì)算步長. DX,DY是.C,h = contour(X,Y,Z);hold onquiver(X,Y,DX,DY)colormap hsvhold off運(yùn)行后見圖（2.2）. 圖（2.1）等高線及其標(biāo)注圖（2.2）等高線和矢量場2.3 梯度線的描繪請預(yù)覽后下載！設(shè)為平面曲線，如果上任意一點(diǎn)處的切線與函數(shù)在該店處的梯度位于同一直線上，則稱為的梯度線?，F(xiàn)在來討論如何作出函數(shù)的梯度線

4、。下面我們一等步長的折線段來近似模擬函數(shù)的梯度線。設(shè)步長為，從點(diǎn)()出發(fā)，沿梯度方向前進(jìn)得到點(diǎn)，即,再從出發(fā)沿梯度線向前進(jìn)得到點(diǎn),依次得到一列點(diǎn)，利用"plot"做出此點(diǎn)集的圖形，即得梯度線的圖形. 例3.作出函數(shù)的梯度線.clear allt=cputime;syms x yS= sym(x2 -y2);Sx=diff(S,'x');Sy=diff(S,'y');x0=1;y0=1;lamda=0.01;i=1;sx(1)=x0;sy(1)=y0;for i=2:400 fx=subs(Sx,x,y,sx(i-1),sy(i-1); fy=

5、subs(Sy,x,y,sx(i-1),sy(i-1); sx(i)=sx(i-1)+lamda*fx./sqrt(fx.2+fy.2); sy(i)=sy(i-1)+lamda*fy./sqrt(fx.2+fy.2);請預(yù)覽后下載！endplot(sx,sy)cputime-t運(yùn)行后見圖（2.3）. 圖（2.3）梯度線2.4 鯊魚襲擊目標(biāo)的前進(jìn)途徑海洋生物學(xué)家發(fā)現(xiàn)，當(dāng)鯊魚在海水中察覺到血液的存在時(shí)，就會(huì)沿著血液濃度增加得最快的方向前進(jìn)去尋找目標(biāo).根據(jù)在海水中世紀(jì)測試的結(jié)果，如果以流血目標(biāo)處作為原點(diǎn)在海面上建立直角坐標(biāo)系，則在海面上點(diǎn)P(x,y)處的血液濃度近似等于（x,y的單位為m,f(x,

6、y)單位的百萬分之一）鍵入x,y = meshgrid(-1:.05:1,-1:.05:1);z =exp(-x.2-2*y.2)/104);C,h = contour(x,y,z);hold on 運(yùn)行后，作為函數(shù)的等高線，得到圖（2.4）.由題設(shè)條件和梯度的性質(zhì)可知，鯊魚襲擊目標(biāo)的前進(jìn)途徑即為的梯度線，下面作出的梯度線，有前面梯度線的繪制可知， 請預(yù)覽后下載！ 圖（2.4）的等高線syms x yS= sym(exp(-x.2-2*y.2)/104);Sx=diff(S,'x');Sy=diff(S,'y');x0=1;y0=1;lamda=0.01;i=1

7、;sx(1)=x0;sy(1)=y0;for i=2:400 fx=subs(Sx,x,y,sx(i-1),sy(i-1); fy=subs(Sy,x,y,sx(i-1),sy(i-1); sx(i)=sx(i-1)+lamda*fx./sqrt(fx.2+fy.2); sy(i)=sy(i-1)+lamda*fy./sqrt(fx.2+fy.2);endplot(sx,sy)hold on得到圖（2.5）：請預(yù)覽后下載！圖（2.5）鯊魚襲擊目標(biāo)的前進(jìn)途徑再運(yùn)用hold on 命令把等高線和梯度線在同一坐標(biāo)系中顯示，得圖（6）. 圖（2.6）等高線和梯度線函數(shù)的梯度線在實(shí)際中有廣泛的應(yīng)用，例如

8、在溫度場總熱量的流動(dòng)也是沿著梯度線方向的. 習(xí)題1.一塊長方形的金屬板，四個(gè)頂點(diǎn)的坐標(biāo)是(1,1)，(5,1)，(1,3)，(5,3)在坐標(biāo)原點(diǎn)處有一個(gè)火焰，它使金屬板受熱假定板上任意一點(diǎn)處的溫度與該點(diǎn)到原點(diǎn)的距離成反比在(3,2)處有一個(gè)青蛙，問這只青蛙應(yīng)沿什么方向爬行才能最快到達(dá)較涼快的地點(diǎn)？(應(yīng)沿由熱變冷變化最驟烈的方向（即梯度方向）爬行) 請預(yù)覽后下載！第三講 函數(shù)的極值、最值及有關(guān)的最優(yōu)化問題水輪機(jī)最優(yōu)化問題最優(yōu)化概念反映了人類實(shí)踐活動(dòng)中十分普遍的現(xiàn)象，即要在盡可能節(jié)省人力、物力和時(shí)間前提下，爭取獲得在可能范圍內(nèi)的最佳效果，因此，最優(yōu)化問題成為現(xiàn)代數(shù)學(xué)的一個(gè)重要課題，涉及統(tǒng)籌、線性規(guī)

9、劃一排序不等式等內(nèi)容。3.1多元函數(shù)的偏導(dǎo)數(shù)1調(diào)用格式一：diff('多元函數(shù)','自變量',n)其中，n 為所求偏導(dǎo)數(shù)的階數(shù)例 1 已知，求和.解 打開文件編輯窗口，在其中輸入下面命令集：pzpx=diff('x2*cos(2*y)','x')p2zpypx=diff(pzpx,'y')p2zpy2=diff('x2*cos(2*y)','y',2)取名為exa9 保存，再在命令窗口中輸入命令exa9，程序運(yùn)行結(jié)果如下：pzpx =2*x*cos(2*y)p2zpypx =-4*x

10、*sin(2*y)p2zpy2 =-4*x2*cos(2*y)即,.2調(diào)用格式二：syms x y z diff(f,自變量,n)例2 已知，求解 在命令行中依次輸入：請預(yù)覽后下載！syms x y zu=sin(x2-y3+5*z);ux=diff(u,x);uxy=diff(ux,y);uxyz=diff(uxy,z);uz3=diff(u,z,3);ux,uxyz,uz3運(yùn)行結(jié)果如下：ux =2*cos(x2-y3+5*z)*xuxyz =30*cos(x2-y3+5*z)*y2*xuz3 =-125*cos(x2-y3+5*z)3.2 隱函數(shù)的導(dǎo)數(shù)在 Matlab 中沒有直接求隱函數(shù)導(dǎo)

11、數(shù)的命令，但可調(diào)用Maple 中求隱函數(shù)導(dǎo)數(shù)的命令，調(diào)用格式如下：maple('implicitdiff(f(u,x,y,z,，)=0,u,x)')例3 求由多元方程x2 + y2 + z2 = xyz所確定的隱函數(shù)解 在命令行中輸入：pzpx=maple('implicitdiff(x2+y2+z2-x*y*z=0,z,x)')運(yùn)行結(jié)果是：pzpx =(2*x-y*z)/(-2*z+x*y)即xy zx yz3.3一元函數(shù)的極（或最）值例1 求 在中的最小值與最大值主程序?yàn)? f='2*exp(-x).*sin(x)' fplot(f,0,8)

12、; %作圖語句請預(yù)覽后下載！ xmin,ymin=fminbnd (f, 0,8) f1='-2*exp(-x).*sin(x)' xmax,zmin=fminbnd (f1, 0,8)ymax=-zmin運(yùn)行結(jié)果： xmin = 3.9270 ymin = -0.0279 xmax = 0.7854 ymax = 0.6448 3.4 多元函數(shù)的極（或最）值在 Matlab 中同樣有求多元函數(shù)的極（或最）小值的函數(shù)，但由于多元函數(shù)的形式比較復(fù)雜，不同情況用到不同的Matlab 函數(shù)若要求多元函數(shù)在某一區(qū)域的極（或最）大值，可轉(zhuǎn)化為求在該區(qū)域內(nèi)的極（或最）小值1非線性無約束情形

13、求極（或最）小值點(diǎn)或極（或最）小值的調(diào)用格式是：x,fval=fminsearch(f,x0) %fminsearch是不能設(shè)定約束范圍的f是被最小化的目標(biāo)函數(shù)名，x0 是求解的初始值向量例 求二元函數(shù)的最值點(diǎn)和最值解 打開文件編輯窗口，在其中輸入下面命令集：%必須對(duì)自變量進(jìn)行轉(zhuǎn)化x=x(1)，y=x(2)Xmin,fmin=fminsearch('2*x(1)3+4*x(1)*x(2)3-10*x(1)*x(2)+x(2)2',0,0);Xmax,Fmin=fminsearch('-2*x(1)3-4*x(1)*x(2)3+10*x(1)*x(2)-x(2)2'

14、;,0,0);fmax=-Fmin;Xmin,fminXmax,fmax取名為exa10保存，再在命令窗口中輸入命令exa10，程序運(yùn)行結(jié)果如下：Xmin =1.0016 0.8335fmin =-3.3241Xmax =-1.0000 1.0000請預(yù)覽后下載！fmax =5.00002非線性有約束情形非線性有約束優(yōu)化問題的數(shù)學(xué)模型如下： 式中，和是向量，和是矩陣，和為函數(shù)，返回量和可以是非線性函數(shù)求極（或最）小值點(diǎn)或極（或最）小值的調(diào)用格式如下：x,fval=fmincon('fun',x0,A,b,Aeq,beq,lb,ub,nonlcon)nonlcon參數(shù)計(jì)算非線性不

15、等式約束c(x)<=0 和非線性等式約束ceq(x)=0例 5 求表面積為的體積最大的長方體體積解 設(shè)長方體的長、寬、高分別為x1、x2、x3，則f(x)=-x(1)*x(2)*x(3)，S.t x(1)*x(2)+x(2)*x(3)+x(3)*x(1)-3=0，x(i)>0,i=1,2,3 建立函數(shù)文件fun1打開文件編輯窗口，在其中輸入下面命令集：function F=fun1(x) %函數(shù)文件必須是function開頭F=-x(1)*x(2)*x(3);單擊“保存”按鈕，自動(dòng)取名為fun1，再擊保存 建立非線性約束函數(shù)文件yceqfunction c,ceq1=yceq1(x

16、)c=x(1)*x(2)+x(2)*x(3)+x(3)*x(1)-3;ceq1=;保存方法同上，自動(dòng)取名為yceq1，再擊保存 編制主程序：打開文件編輯窗口，在其中輸入下面命令集：x0=3;3;3; %給長寬高一個(gè)初值A(chǔ)=;b=;Aeq=;beq=;請預(yù)覽后下載！lb=0,0,0;ub=;xmax,fmin=fmincon('fun1',x0,A,b,Aeq,beq,lb,ub,'yceq1'); %函數(shù)要加單引號(hào)Vmax=-fmin;xmax,Vmax取名為exa11保存，再在命令窗口中輸入命令exa11，程序運(yùn)行結(jié)果如下：xmax =1.00001.0000

17、1.0000Vmax =1.0000例6求三元函數(shù)滿足約束條件的最小值。 建立函數(shù)文件fun3打開文件編輯窗口，在其中輸入下面命令集：function F=fun3(x) %函數(shù)文件必須是function開頭F=-x(1)*x(2)*x(3);單擊“保存”按鈕，自動(dòng)取名為fun3，再擊保存 建立非線性約束函數(shù)文件yceq3把約束條件改寫為因?yàn)閮蓚€(gè)約束均為線性，所以符合矩陣不等式的形式，其中，.(2) 編制主程序：打開文件編輯窗口，在其中輸入下面命令集：x0=10;10;10; %在此點(diǎn)處開始尋找A=-1 -2 -2; 1 2 2;b=0; 72; xmax,fmin=fmincon('

18、fun3',x0,A,b); %函數(shù)要加單引號(hào)xmax,fmin請預(yù)覽后下載！取名為zuizhi3保存，再在命令窗口中輸入命令zuizhi3，程序運(yùn)行結(jié)果如下： xmin = 24.0000 12.0000 12.0000fmin = -3.4560e+003例 7水輪機(jī)最優(yōu)化問題眾所周知，三峽水電站是中國最大的水電站。已知，水是通過輸水管從水壩送到發(fā)電所，而輸水管中水流速度的變化依賴于外界條件。若發(fā)電所有三個(gè)不同的水電渦輪，每個(gè)渦輪都有一個(gè)已知的、特定的功率輸出函數(shù)來提供發(fā)電所需的功率。而發(fā)電功率又是輸送到渦輪的水流速度的函數(shù)。將進(jìn)入發(fā)電所的水分成體積不同的三部分，分別到達(dá)每個(gè)渦輪，

19、因此，我們的目標(biāo)就是，若給定任意一個(gè)總的水流速度，如何分流才能使總的輸出功率最大。根據(jù)試驗(yàn)證據(jù)和伯努利方程，每個(gè)渦輪的輸出功率可用二次方程式模型描述，其中，表示通過第個(gè)渦輪機(jī)的水流流速，單位為立方英尺/秒；表示第個(gè)渦輪機(jī)的輸出功率，單位為千瓦；表示進(jìn)入發(fā)電所的總的水流速度，單位為立方英尺/秒。1、如果三個(gè)渦輪機(jī)一起工作，需要確定通過每個(gè)渦輪機(jī)的水流速度，使得總輸出功率最大。當(dāng)來流速度一定時(shí)，可以用拉格朗日乘子解出總的輸出功率在滿足約束條件下的最大值，這個(gè)最大值一定是來流速度請預(yù)覽后下載！的函數(shù)。2、當(dāng)取何值時(shí)問題1中的結(jié)果才是有效的？3、當(dāng)進(jìn)入發(fā)電所的水流速度為2500ft3/s（1立方英尺(

20、ft3)=0.0283立方米(m3)）,該如何分流才能使總輸出功率最大？（1）用1的結(jié)果計(jì)算；（2）直接用Matlab中的非線性有約束優(yōu)化問題的fmincon函數(shù)求解4、現(xiàn)在我們是讓三個(gè)渦輪一起工作，那有沒有可能在某種情況下一個(gè)渦輪工作能產(chǎn)生更大的電量呢？做出三個(gè)功率輸出函數(shù)的圖像，并討論當(dāng)來流速度為1000 ft3/s時(shí)應(yīng)如何分配水流，是讓三個(gè)渦輪一起工作還是只用一個(gè)？（如果只用其中一個(gè)，那么該用哪一個(gè)）當(dāng)來流速度為600 ft3/s時(shí)又該如何分配呢？5、對(duì)于某個(gè)來流速度而言，或許選擇兩個(gè)渦輪工作最好。若來流速度為1500 ft3/s，那么應(yīng)該選擇哪兩個(gè)渦輪一起工作最好？用拉格朗日乘子計(jì)算，

21、如何給兩個(gè)渦輪分配來流才能使輸出功率最大，這樣分配是否比同時(shí)啟動(dòng)三個(gè)渦輪更有效？6、如果來流速度為3400 ft3/s，那么該如何分配？解答：1、當(dāng)來流速度一定時(shí)，若同時(shí)啟動(dòng)三個(gè)渦輪機(jī)，那么，根據(jù)拉格朗日乘子法計(jì)算總輸出功率滿足約束條件的極大值。設(shè)拉格朗日函數(shù)為.由此得到微分方程組求解命令如：請預(yù)覽后下載！eq1=sym('(0.1277-8.16*10(-5)*Q1)*(170-1.6*10(-6)*QT2)+Q4=0');eq2=sym('(0.1358-9.38*10(-5)*Q2)*(170-1.6*10(-6)*QT2)+Q4=0');eq3=sym(

22、'(0.1380-7.68*(10(-5)*Q3)*(170-1.6*(10(-6)*QT2)+Q4=0');eq4=sym('Q1+Q2+Q3-QT=0');Q1,Q2,Q3,Q4=solve(eq1,eq2,eq3,eq4)運(yùn)行得：Q1 = .34101340604408089070665757782322*QT-75.182723623418919942437324850413 Q2 = .29665985003408316291751874573960*QT+20.949784139968189047943649170643 Q3 = 54.232939

23、483450730894493675679770+.36232674392183594637582367643717*QT取四位有效數(shù)字得： (3.1)每個(gè)分水流速度都是總水流速度的函數(shù)。顯然，總輸出功率也是的函數(shù)。2、根據(jù)(3.1)式，及（）的取值范圍，可以解得， 請預(yù)覽后下載！因此，當(dāng)同時(shí)啟動(dòng)三個(gè)渦輪機(jī)時(shí)，總水流速度必須滿足，才能應(yīng)用(3.1)式計(jì)算。3、（1）當(dāng)總流速度 ft3/s時(shí)，在問題2中算出的取值范圍之內(nèi)，應(yīng)用(3.1)式算出當(dāng) ft3/s，ft3/s， ft3/s時(shí)總輸出功率取得最大值，且最大值為28412千瓦。（2）計(jì)算機(jī)求解首先建立函數(shù) 建立函數(shù)文件fun2打開文件編輯窗口

24、，在其中輸入下面命令集：function F=fun2(x) %函數(shù)文件必須是function開頭y1=(-18.89+0.1277*x(1)-4.08*10(-5)*x(1)2)*(170-1.6*10(-6)*(x(1)+x(2)+x(3)2);y2=(-24.51+0.1358*x(2)-4.69*10(-5)*x(2)2)*(170-1.6*10(-6)*(x(1)+x(2)+x(3)2);y3=(-27.02+0.1380*x(3)-3.84*(10(-5)*x(3)2)*(170-1.6*10(-6)*(x(1)+x(2)+x(3)2);F=-(y1+y2+y3); 建立非線性約束

25、函數(shù)文件yceq2function c,ceq2=yceq2(x)c=x(1)+x(2)+x(3)-2500;ceq2=; 編制主程序：打開文件編輯窗口，在其中輸入下面命令集：x0=20;30;10; %從該點(diǎn)處開始解決問題 A=;b=;Aeq=;beq=; %如果沒有不等式存在,則令A(yù)=,b= lb=250,250,250;ub=1110,1110,1225;Qmax,fmin=fmincon('fun2',x0,A,b,Aeq,beq,lb,ub,'yceq2');%函數(shù)要加單引號(hào)Wmax=-fmin;Qmax,Wmax請預(yù)覽后下載！運(yùn)行得Qmax = 77

26、7.3508 762.5994 960.0498Wmax = 2.8412e+0044、做出每個(gè)渦輪機(jī)單獨(dú)工作時(shí)其輸出功率關(guān)于來流速度的函數(shù)圖像，程序如下:x1=250:0.01:1110;x2=250:0.01:1110;x3=250:0.01:1250;q=2500;y1=(-18.89+0.1277.*x1-4.08.*10.(-5).*x1.*x1).*(170-1.6.*10.(-6).*q.*q);y2=(-24.51+0.1358.*x2-4.69.*10.(-5).*x2.*x2).*(170-1.6.*10.(-6).*q.*q);y3=(-27.02+0.1380.*x3-

27、4.08.*10.(-5).*x3.*x3).*(170-1.6.*10.(-6).*q.*q);plot(x1,y1,':') hold onplot(x2,y2,'-')plot(x3,y3,'-')xlabel('流速');ylabel('輸出功率');運(yùn)行得圖（3.1）請預(yù)覽后下載！圖（3.1）三個(gè)渦輪機(jī)單獨(dú)工作時(shí)其輸出功率關(guān)于來流速度的函數(shù)圖像如圖（3.1）所示，“:”線表示第一個(gè)渦輪機(jī)隨水流速度的輸出功率，“-”線表示第二個(gè)渦輪機(jī)隨水流速度的輸出功率，而“”線則表示第三個(gè)渦輪機(jī)的輸出功率。觀察上圖可知，

28、當(dāng)來流速為1000ft3/s時(shí)，“”線在最上方，如果只用一個(gè)渦輪機(jī)工作，應(yīng)選擇第三個(gè)渦輪機(jī)才能使輸出功率最大。而當(dāng)來流速度為600 ft3/s時(shí)，“:”線在最上方，即第一個(gè)渦輪機(jī)輸出功率最大。當(dāng)來流速度為1000 ft3/s時(shí)，若只啟動(dòng)第三個(gè)渦輪機(jī)，根據(jù)(3.1)式，及的取值范圍解得 ft3/s， ft3/s， ft3/s時(shí)算得總的輸出功率最大，最大值為千瓦。而如果只用第三個(gè)渦輪機(jī)工作，則總的輸出功率為千瓦。所以應(yīng)只選擇第三個(gè)渦輪機(jī)單獨(dú)工作最好。當(dāng)來流速度為600 ft3/s時(shí)，根據(jù)的取值范圍，不能同時(shí)啟動(dòng)三個(gè)渦輪機(jī)，若只啟動(dòng)一個(gè)，則選擇第一個(gè)渦輪機(jī)單獨(dú)工作輸出功率最大，最大值為7292千瓦。

29、5、觀察圖（3.1）可以看出，當(dāng)來流速度超過500 ft3/s后，“-”一直處于最下方，也就是說，在同樣的水流速度下，第二個(gè)渦輪機(jī)的輸出功率最小。當(dāng)來流速度為1500 ft3/s時(shí)，若只用兩個(gè)渦輪機(jī)工作，應(yīng)該選擇第一和第三個(gè)渦輪機(jī)一起工作最好。此時(shí)用拉格朗日乘子法計(jì)算總輸出功率請預(yù)覽后下載！滿足約束條件的極大值。設(shè)拉格朗日函數(shù)由方程此時(shí)解得 (3.2)根據(jù)和的取值范圍，從(3.2)式解得適合此式的的取值范圍為。當(dāng)=1500 ft3/s時(shí)，用Matlab求解的解得 ft3/s， ft3/s，總的輸出功率最大，最大值為18208 ft3/s。而此時(shí)若同時(shí)啟動(dòng)三個(gè)渦輪機(jī)一起工作，根據(jù)（3.1）式計(jì)算

30、出 ft3/s， ft3/s， ft3/s，總輸出功率為16539千瓦。顯然，此時(shí)只啟動(dòng)兩個(gè)渦輪機(jī)工作更有效。根據(jù)(3.2)式和(3.1)式，繪出總輸出功率關(guān)于總來流速度的函數(shù)圖像。程序如下：l=0.01;%步長q0=953;qe=3636;q1=q0:l:qe;x1=-65.0253+0.4848.*q1;x3=65.0253+0.5152.*q1;y1=(-18.89+0.1277.*x1-4.08.*10.(-5).*x1.*x1).*(170-1.6.*10.(-6).*q1.*q1);y3=(-27.02+0.1380.*x3-3.84.*10.(-5).*x3.*x3).*(170

31、-1.6.*10.(-6).*q1.*q1);z=y1+y3;請預(yù)覽后下載！plot(q1,z,':','LineWidth',2)hold onq1=q0:l:qe;x1=-75.1827+0.341.*q1;x2=-20.9498+00.2967.*q1;x3=54.233+0.3623.*q1;y1=(-18.89+0.1277.*x1-4.08.*10.(-5).*x1.*x1).*(170-1.6.*10.(-6).*q1.*q1);y2=(-24.51+0.1358.*x2-4.69.*10.(-5).*x2.*x2).*(170-1.6.*10.(

32、-6).*q1.*q1);y3=(-27.02+0.1380.*x3-3.84.*10.(-5).*x3.*x3).*(170-1.6.*10.(-6).*q1.*q1);y=y1+y2+y3;h=z-y;aa=find(h=min(abs(h);%找交點(diǎn)所對(duì)應(yīng)的迭代次數(shù)q2=l*aa+q0 %計(jì)算交點(diǎn)所對(duì)應(yīng)的流速y(aa) %計(jì)算交點(diǎn)所對(duì)應(yīng)的輸出功率text(q2,y(aa),'*(2.1038e+003,2.3898e+004)')plot(q1,y,'-','LineWidth',2)xlabel('流速');ylabel(

33、'輸出功率');運(yùn)行得圖（3. 2）。虛線表示同時(shí)啟動(dòng)第一和第三個(gè)渦輪機(jī)時(shí)總輸出功率隨總水流速度的變化情況，而實(shí)線則表示三個(gè)渦輪機(jī)同時(shí)啟動(dòng)時(shí)總輸出功率隨總水流速度的變化情況。觀察圖（3. 2），兩條曲線大約在處有個(gè)交點(diǎn)，說明，當(dāng)時(shí)，同時(shí)啟動(dòng)兩個(gè)渦輪機(jī)比同時(shí)啟動(dòng)三個(gè)能獲得更大的總輸出功率，而當(dāng)后，三個(gè)渦輪同時(shí)啟動(dòng)較好。請預(yù)覽后下載！圖（3. 2）總輸出功率關(guān)于總來流速度的函數(shù)圖像6、請同學(xué)們自己解答。習(xí)題1. 求下列函數(shù)的極小點(diǎn) 1) ；2) ；3) . 第1），2）題的初始點(diǎn)可任意選取，第3）題的初始點(diǎn)取為.2. 梯子長度問題一樓房的后面是一個(gè)很大的花園. 在花園中緊靠著樓房有

34、一個(gè)溫室,溫室伸入花園2m,高3m,溫室正上方是樓房的窗臺(tái). 清潔工打掃窗臺(tái)周圍,他得用梯子越過溫室,一頭放在花園中,一頭靠在樓房的墻上. 因?yàn)闇厥沂遣荒艹惺芴葑訅毫Φ?所以梯子太短是不行的.現(xiàn)清潔工只有一架7m長的梯子,你認(rèn)為它能達(dá)到要求嗎? 能滿足要求的梯子的最小長度為多少?請預(yù)覽后下載！3. 陳酒出售的最佳時(shí)機(jī)問題某酒廠有批新釀的好酒,如果現(xiàn)在就出售,可得總收入=50萬元(人民幣),如果窖藏起來待來日(第年)按陳酒價(jià)格出售,第n年末可得總收入(萬元),而銀行利率為=0.05,試分析這批好酒窖藏多少年后出售可使總收入的現(xiàn)值最大. (假設(shè)現(xiàn)有資金萬元,將其存入銀行,到第年時(shí)增值為萬元,則稱X

35、為的現(xiàn)值.)并填下表：第一種方案：將酒現(xiàn)在出售，所獲50萬元本金存入銀行；第二種方案：將酒窖藏起來，待第年出售。（1） 計(jì)算15年內(nèi)采用兩種方案，50萬元增值的數(shù)目并填入表1，2中；（2） 計(jì)算15年內(nèi)陳酒出售后總收入的現(xiàn)值填入表3中。表1 第一種方案第1年第2年第3年第4年第5年第6年第7年第8年第9年第10年第11年第12年第13年第14年第15年表2 第二種方案第1年第2年第3年第4年第5年第6年第7年第8年第9年第10年第11年第12年第13年第14年第15年表3 陳酒出售后的現(xiàn)值第1年第2年第3年第4年第5年請預(yù)覽后下載！第6年第7年第8年第9年第10年第11年第12年第13年第14

36、年第15年4. 請解答例7中的第6個(gè)問題。第四講 微分方程模型 導(dǎo)彈系統(tǒng)改進(jìn)、交通管理色燈及自動(dòng)噴水滅火裝置設(shè)計(jì)實(shí)際應(yīng)用問題通過數(shù)學(xué)建模所歸納而得到的方程，絕大多數(shù)都是微分方程. 另一方面，能夠求解的微分方程也是十分有限的，特別是高階方程和偏微分方程（組）.這就要求我們必須研究微分方程組的解法，既要研究微分方程組的解析解法（精確解），更要研究微分方程組的數(shù)值解法（近似解）.請預(yù)覽后下載！ 對(duì)微分方程組的解析解法，Matlab有專門的函數(shù)可以用，本實(shí)驗(yàn)將作一定的介紹.然后建模實(shí)驗(yàn)研究導(dǎo)彈系統(tǒng)的改進(jìn)等幾個(gè)實(shí)際問題的微分方程建模及計(jì)算機(jī)求解.4.1 預(yù)備知識(shí)微分方程、通解、特解的概念補(bǔ)充一般說來，微

37、分方程就是聯(lián)系自變量、未知函數(shù)以及未知函數(shù)的某些導(dǎo)數(shù)之間的關(guān)系式.如果其中的未知函數(shù)只與一個(gè)自變量有關(guān)，則稱為常微分方程；如果未知函數(shù)是兩個(gè)或兩個(gè)以上自變量的函數(shù)，并且在方程中出現(xiàn)偏導(dǎo)數(shù)，則稱為偏微分方程. 本書所介紹的都是常微分方程，有時(shí)就簡稱微分方程或方程. 微分方程的解中可以包含任意常數(shù)，其中任意常數(shù)的個(gè)數(shù)可以多到與方程的階數(shù)相等，也可以不含任意常數(shù). 我們把n階常微分方程的含有n個(gè)獨(dú)立的任意常數(shù)的解 ，稱為該方程的通解，如果方程的解不包含任意常數(shù)，則稱它為特解.4.2相關(guān)函數(shù)命令及簡介1.dsolve( ) 求微分方程的解析解 Matlab語言的符號(hào)運(yùn)算工具箱提供了一個(gè)常系數(shù)線性微分方

38、程求解的實(shí)用函數(shù)dsolve( ),該函數(shù)允許用字符串的形式描述微分方程及初值、邊值條件，最終得出解析解。調(diào)用格式指明自變量，其中 可以描述微分方程，又可以描述初始條件或邊界條件。D4y表示，D2y(2)=3表示。2.T,Y=solver(odefun,tspan,y0),求微分方程的數(shù)值解.說明:(1)其中solver為命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb之一.請預(yù)覽后下載！(2)odefun是顯式常微分方程：（3）在積分區(qū)間上，從到，用初始條件求解。（4）要獲得問題在其他指定時(shí)間點(diǎn)上的解，則令（要求是單調(diào)的）。3. inline

39、（）：建立一個(gè)內(nèi)聯(lián)函數(shù).格式：inline（expr, var1, var2,）,注意括號(hào)里的表達(dá)式要加引號(hào).4.3計(jì)算實(shí)驗(yàn)：（自學(xué)）直接可以用Matlab求微分方程精確解的例子。 例1假設(shè)輸入信號(hào)為 ,求下面微分方程的通解 （1） syms t;u=exp(-5*t)*cos(2*t+1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*usyms t y;y=dsolve('D4y+10*D3y+35*D2y+50*Dy+24*y=','87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10')

40、latex(y)例2 假設(shè) 方程（1）滿足初始條件 求特解>>syms t;u=exp(-5*t)*cos(2*t+1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*usyms t y;y=dsolve('D4y+10*D3y+35*D2y+50*Dy+24*y=','87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10','y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')請預(yù)覽后下載！

41、例3 求解下面線性微分方程組 x,y=dsolve('D2x+2*Dx=x+2*y-exp(-t)','Dy=4*x+3*y+4*exp(-t)')例4Lorenz模型其中設(shè). 令初值為.試求解該微分方程組.求解：編程一：function lo()t_final=100;x0=0;0;1e-10;t,x=ode45(lorenzeq,0,t_final,x0);plot(t,x),figure;plot3(x(:,1),x(:,2),x(:,3);axis(10 42 -20 20 -20 25);comet3(x(:,1),x(:,2),x(:,3);func

42、tion xdot=lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3)-10*x(2)+10*x(3) -x(1)*x(2)+28*x(2)-x(3);編程二：f1=inline('-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)','t','x');請預(yù)覽后下載！t_final=100;x0=0;0;1e-10;t,x=ode45(f1,0,t_final,x0);plot(t,x),figure;plot3(x(:,1),x(:,2),x(:,3)

43、;axis(10 42 -20 20 -20 25);comet3(x(:,1),x(:,2),x(:,3);例1 求解微分方程例2求微分方程在初始條件下的特解，并畫出解函數(shù)的圖形。例3求微分方程組在初始條件下的特解，并畫出解函數(shù)的圖形。例4求解微分方程初值問題的數(shù)值解，求解范圍為0,0.5.例5求解描述振蕩器的經(jīng)典的VerderPol微分方程4.4建模實(shí)驗(yàn)：例1 導(dǎo)彈系統(tǒng)的改進(jìn)海軍方面要求改進(jìn)現(xiàn)有的艦對(duì)艦系統(tǒng).目前的電子系統(tǒng)能迅速測出敵艦的種類、位置以及敵艦行駛的速度和方向，且導(dǎo)彈自動(dòng)制導(dǎo)系統(tǒng)能保證在發(fā)射后的任意時(shí)刻都能對(duì)準(zhǔn)目標(biāo).根據(jù)情報(bào)，這種敵艦?zāi)茉谖臆娕灠l(fā)射導(dǎo)彈后T分鐘作出反應(yīng)并摧毀導(dǎo)彈

44、.現(xiàn)在要求改進(jìn)電子導(dǎo)彈系統(tǒng)使能自動(dòng)計(jì)算出敵艦是否在有效打擊范圍之內(nèi).設(shè)我艦發(fā)射導(dǎo)彈時(shí)位置在坐標(biāo)原點(diǎn)，敵艦在軸正向km處，其行駛速度為km/h,方向與軸夾角為，導(dǎo)彈水平飛行線速度為km/h. 問題的關(guān)鍵是求出導(dǎo)彈擊中敵艦的時(shí)間.特別地，在導(dǎo)彈系統(tǒng)中設(shè)km/h, 請預(yù)覽后下載！ km/h,h. 現(xiàn)在求的有效范圍.求解：設(shè)時(shí)刻導(dǎo)彈位置為.則 （4.1）易知時(shí)刻敵艦位置為，為了保持對(duì)準(zhǔn)目標(biāo)，導(dǎo)彈軌跡切線方向應(yīng)為 （4.2）由（4.1）和（4.2）得下列微分方程， (4.3) （4.4）初始條件對(duì)給定的進(jìn)行計(jì)算.當(dāng)滿足， （4.5）則認(rèn)為已擊中目標(biāo).如果,則敵艦在打擊范圍內(nèi),可以發(fā)射. 有兩個(gè)極端情形容

45、易算出. 若即敵艦正好背向行駛,即軸正向.那么導(dǎo)彈直線飛行,擊中時(shí)間 得. 若即迎面駛來,類似地.一般地,應(yīng)有請預(yù)覽后下載！%M函數(shù)missilefun.mfunction dy=missilefun(t,y,a,b,d,theta)dydx=(a*t*sin(theta)-y(2) +1e-8)/.(abs(d+a*t*cos(theta)-y(1)+1e-8);dy(1)=b/(1+dydx2)0.5;dy(2)=b/(1+dydx(-2)0.5;dy=dy(:);在指令窗口執(zhí)行下列腳本文件得-5.7410. %M腳本missilefunªclear;close;a=90;b=4

46、50;d=50;theta=pi/2;t,y=ode45(missilefun,0,0.1,0 0,a,b,d,theta);plot(y(:,1),y(:,2);max(y(:,1)-d-a*t*cos(theta)由于在內(nèi)，（5）式不成立，所以敵艦不在有效打擊范圍，應(yīng)等近一些再發(fā)射.例2. 交通管理色燈中黃燈應(yīng)亮多長時(shí)間？讓我們來考慮這樣一個(gè)問題：紅綠燈在亮紅燈之前黃燈應(yīng)亮多長時(shí)間？在交通管理中，定期的亮一段時(shí)間黃燈是為了讓那些正行駛在交叉路口上或距交叉路口太近以致無法停下的車輛通過路口。這樣，紅綠燈應(yīng)保持足夠長時(shí)間的黃燈，便那些無法停下的駕駛員有機(jī)會(huì)在黃燈亮的時(shí)候通過路口。對(duì)于一位駛近交

47、叉路口的駕駛員來說，萬萬不可處于這樣的進(jìn)退兩難的境地：要安全停車，離路口太近，而要在紅燈亮之前通過路口又顯得太遠(yuǎn)了。關(guān)于這個(gè)時(shí)間的計(jì)算，有一個(gè)粗糙的“經(jīng)驗(yàn)方法”：對(duì)法定迫近速度的每個(gè)10mi/h亮1s黃燈。我們來看理論計(jì)算是否證實(shí)這個(gè)經(jīng)驗(yàn)方法。駛近交叉路口的駕駛員，在看到黃色信號(hào)后要作出決定：是停車還是通過路口。如果他以法定速度（或低于法定速度）行駛，當(dāng)決定停車時(shí)，他必須有足夠的停車距離。當(dāng)決定通過路口時(shí)，他必須有足夠的時(shí)間使他能夠完全通過路口，這包括作出停車決定的時(shí)間（反應(yīng)時(shí)間）以及通過停車所需的最短距離的駕駛時(shí)間。能夠很快看到黃燈的駕駛員可以利用剎車距離將車停下。于是，黃燈狀態(tài)應(yīng)持續(xù)的時(shí)間

48、包括駕駛員的反應(yīng)時(shí)間、他通過交叉路口的時(shí)間以及停車所需的時(shí)間（剎車距離）。有了這么多的時(shí)間，駕駛員就能夠在剎車距離內(nèi)安全停車。請預(yù)覽后下載！如果法定速度為，交叉路口的寬度為,典型的車身長度為，那么通過路口的時(shí)間為。（注意車的尾部必須通過路口，這樣，路口的實(shí)際長度就是.)現(xiàn)在，我們來計(jì)算剎車距離。注意到實(shí)際的剎車和停車過程是相當(dāng)復(fù)雜的，駕駛員首先開加速器，然后不同程度地用勁踩住剎車踏板（也許用了氣動(dòng)剎車）使車子減速，直到停止。這些過程的大部分是很難用精確地模型來描述的，因而我們將回避它們。通過在剎車過程引入一個(gè)抵抗摩擦力，我們將假定剎車效果可用模型來反映。設(shè)為汽車重量，為摩擦系數(shù)。根據(jù)定義，對(duì)汽

49、車的制動(dòng)力為,其方向與運(yùn)動(dòng)方向相反. 汽車在停車過程中，行駛的距離可通過解下面的微分方程求得，這個(gè)方程反映的是在常力作用下的直線運(yùn)動(dòng)： (4.6) 其中是重力加速度。賦予的正確條件是時(shí)，以及. 于是，剎車距離就是直到時(shí)汽車駛過的距離。 求解：在的條件下對(duì)（4.6）積分，得 (4.7)因此，當(dāng)時(shí)，速度為零。在條件下對(duì)（4.7）積分，得請預(yù)覽后下載！ (4.8)時(shí)，的值是 （4.9）上述求解的過程可以用計(jì)算機(jī)求解：（huangdeng1.m）syms t x f g v0 y;x=dsolve('D2x=','-f*g','Dx(0)=v0',

50、9;x(0)=0')x =-1/2*f*g*t2+v0*t注意，對(duì)必須用到的單位要謹(jǐn)慎：對(duì)我們所涉及的距離，交通工程師是用來度量的（ft為英尺的簡寫，1英尺=0.3048米），而速度則通常用（mi/h=0.44704m/s）來度量。重力加速度= 32.16在作任何計(jì)算之前，讀者應(yīng)將轉(zhuǎn)換成. 我們將計(jì)算一下黃燈狀態(tài)其中是駕駛員的反應(yīng)時(shí)間。于是 假設(shè) 另外，我們將接受公路工程師提出的具有代表性的(參看練習(xí)1）當(dāng)=30，40以及50時(shí)，黃燈時(shí)間如表1所示。表中同時(shí)也給出了經(jīng)驗(yàn)法的值。請預(yù)覽后下載！黃燈周期與法定速度的關(guān)系編寫程序如下（huangdeng2.m）v0=1:0.1:50;f=0.

51、2;g=9.8/0.3048;I=30;L=15;T=1;A=(v0*0.44704/0.3048)./(2*f*g)+(I+L)./(v0*0.44704/0.3048)+T;text(30,5.46,'*')plot(v0,A)xlabel('v_0(mi/h)');ylabel('A(s)');關(guān)于的圖像描繪出來，如圖 4.1所示。圖4.1 黃燈周期與法定速度的關(guān)系 表 1（mi/h) A 經(jīng)驗(yàn)法30 5.46s 3s40 6.35s 4s50 7.34s 5s 我們注意到，經(jīng)驗(yàn)法的結(jié)果一律比我們預(yù)測的黃燈狀態(tài)短些。這使人想起，許多交叉路口

52、的設(shè)計(jì)很可能使車輛在紅綠燈轉(zhuǎn)為紅燈時(shí)正處于交叉路口上。請預(yù)覽后下載！ 即使給了充分的停車時(shí)間，仍有許多汽車司機(jī)企圖加速想搶在紅燈亮之前沖過十字路口。這當(dāng)中的大多數(shù)駕駛員不知道（有些人甚至不予注意）什么時(shí)候綠燈轉(zhuǎn)紅。有一種“倒數(shù)”型的紅綠燈可以部分地解決這個(gè)問題，在黃燈亮的最后幾秒鐘內(nèi)，黃燈上顯示出一串倒著數(shù)的數(shù)字，它們準(zhǔn)確地向駕駛員警告紅綠燈何時(shí)將變色。這種系統(tǒng)幾年前就在休斯頓試用了，它成功地降低了事故發(fā)生率。參看練習(xí)6. 習(xí)題1 (廣告效應(yīng))某公司生產(chǎn)一種耐用消費(fèi)品,市場占有率為時(shí)開始做廣告,一段時(shí)間的市場跟蹤調(diào)查后,該公司發(fā)現(xiàn):單位時(shí)間內(nèi)購買人口百分比的相對(duì)增長率與當(dāng)時(shí)還沒有買的百分比成正比,且估得此比例系數(shù)為0.5. (1)建立該問題的數(shù)學(xué)模型,分別求其解析解和數(shù)值解,并作比較. (2)廠家問:要做多少時(shí)間廣告,可使市場占有率達(dá)到