數(shù)學(xué)實驗 課件 第9、10章 常微分方程實驗、優(yōu)化方法實驗

9.1常微分方程的解析解

定義9.1

含有未知函數(shù)的導(dǎo)數(shù)的方程稱為微分方程.

如果未知函數(shù)是一元函數(shù)，稱為常微分方程.

常微分方程的一般形式為.

由幾個微分方程聯(lián)立而成的方程組稱為微分方程組．

定義9.2

微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階．若方程中未知函數(shù)及其各階導(dǎo)數(shù)都是一次的，稱為線性常微分方程，一般表示為.

若上式中的系數(shù)

均與

x無關(guān)，稱之為常系數(shù).

二階及二階以上的微分方程稱為高階微分方程．

一階常微分方程與高階微分方程可以互化，已給一個n階方程設(shè)

，

可將上式化為一階方程組

反過來，在許多情況下，一階微分方程組也可化為高階方程.

所以一階微分方程組與高階常微分方程的理論與方法在許多方面是相通的，一階常系數(shù)線性微分方程組也可用特征根法求解．

MATLAB中主要用dsolve求符號解析解，調(diào)用格式如下：

S=dsolve(eqn,cond)

求帶初始條件的微分方程特解，若缺省cond，則求微分方程的通解.利用“diff”和“==”表示微分方程．例如diff(y,x)==y表示dy/dx=y．例9.1

求下列微分方程的解析解（1）（2）（3）解（1）>>symsy(t)ab

>>eqn=diff(y,t)==a*y+b;

>>dsolve(eqn)

ans=

-(b-C1*exp(a*t))/a可知微分方程的通解為

，其中

，C1是任意常數(shù)．（2）>>symsy(x)>>eqn=diff(y,x,2)==sin(2*x)-y;>>Dy=diff(y,x);>>cond=[y(0)==0,Dy(0)==1];>>dsolve(eqn,cond)ans=(5*sin(x))/3-sin(2*x)/3可知微分方程的特解為.（3）>>symsf(t)g(t)>>eqns=[diff(f,t)==f+g,diff(g,t)==g-f];>>cond=[f(0)==1,g(0)==1];>>S=dsolve(eqns,cond)S=包含以下字段的struct:g:[1×1sym]f:[1×1sym]>>S.fans=exp(t)*cos(t)+exp(t)*sin(t)>>S.gans=exp(t)*cos(t)-exp(t)*sin(t)可知微分方程組的特解為.9.2微分方程的數(shù)值解

除常系數(shù)線性微分方程可用特征根法求解，少數(shù)特殊方程可用初等積分法求解外，大部分微分方程無解析解，應(yīng)用中主要依靠數(shù)值解法．考慮一階常微分方程初值問題其中所謂數(shù)值解法，就是尋求y(x)在一系列離散節(jié)點

上的近似值yk

，

．稱

為步長，通常取為常量h．最簡單的數(shù)值解法是Euler法．9.2.1歐拉法

Euler法的思路極其簡單：在節(jié)點處用差商近似代替導(dǎo)數(shù)

這樣導(dǎo)出計算公式（稱為Euler格式）

它能求解各種形式的微分方程.Euler法也稱折線法．

Euler方法只有一階精度，改進方法有二階Runge-Kutta法、四階Runge-Kutta法、五階Runge-Kutta-Felhberg法和先行多步法等，這些方法可用于解高階常微分方程（組）初值問題.邊值問題采用不同方法，如差分法、有限元法等.

數(shù)值算法的主要缺點是它缺乏物理理解．9.2.2龍格-庫塔法

歐拉公式易于計算，但精度不高，收斂速度慢．在實際應(yīng)用中，我們采用龍格-庫塔（Runge-Kutta）方法，基本思想為進一步可得在MATLAB中，利用ode23，ode45求微分方程數(shù)值解．

[t,y]=ode23(odefun,tspan,y0)

[t,y]=ode45(odefun,tspan,y0)

求微分方程組y′=f(t,y)從t0到tf的積分，初始條件為y0.其中tspan=[t0tf].解數(shù)組y中的每一行都與列向量t中返回的值相對應(yīng)．

ode45是最常用的求解微分方程數(shù)值解的命令，采用四階和五階Runge-Kutta算法，是一種自適應(yīng)步長（變步長）的常微分方程數(shù)值解法，對于剛性方程組不宜采用.ode23與ode45類似，采用二階和三階Runge-Kutta算法，只是精度低一些．ode45是解決數(shù)值解問題的首選方法．例9.2

求解微分方程先求解析解，再求數(shù)值解，并進行比較.解>>clear;>>symsy(t);>>eqn=diff(y,t)==-y+t+1;>>cond=y(0)==1;>>dsolve(eqn,cond)ans=t+exp(-t)可得解析解為

．下面再求其數(shù)值解，先編寫M文件fun92.m．%M函數(shù)fun92.mfunctionf=fun92(t,y)f=-y+t+1;再用命令>>clear;t=0:0.1:1;>>y=t+exp(-t);plot(t,y);

%繪制解析解的圖形>>holdon;

%保留已經(jīng)畫好的圖形，如果下面再畫圖，兩個圖形和并在一起>>[t,y]=ode45(@fun92,[0,1],1);>>plot(t,y,'ro');

%繪制數(shù)值解圖形>>xlabel('t'),ylabel('y')結(jié)果見圖9-1．

由圖9-1所示，解析解和數(shù)值解吻合得很好．圖9-1解析解與數(shù)值解例9.3已知方程當

時，上面方程可化為求上面方程的解析解和數(shù)值解．解先求解析解，>>symsy(t)>>eqn=diff(y,t,2)==9.8*sin(t);>>Dy=diff(y,t);>>cond=[y(0)==15,Dy(0)==0];>>dsolve(eqn,cond)ans=(49*t)/5-(49*sin(t))/5+15可知方程的解析解為

.下面求數(shù)值解．令

可將原方程化為如下方程組建立函數(shù)文件fun93.m如下%M文件fun93.mfunctionf=fun93(t,y)f=[y(2),9.8*sin(y(1))]';

%f向量必須為一列向量運行MATLAB代碼>>clear;close;>>[t,y]=ode45(@fun93,[0,10],[15,0]);>>plot(t,y(:,1));

%畫

隨時間變化圖,y(:2)則表示

的值>>xlabel('t'),ylabel('y1')結(jié)果見圖9-2．

由圖9-2可見，

隨時間t周期變化．圖9-2數(shù)值解

例9.4

Lotka-Volterra方程，也即捕食者-獵物模型的一對一階常微分方程

（9-1）

變量

x

和

y

分別計算獵物和捕食者的數(shù)量.當沒有捕食者時，獵物數(shù)量將增加，當獵物匱乏時，捕食者數(shù)量將減少.使用初始條件x(0)=y(0)=20，使捕食者和獵物的數(shù)量相等．求當α=0.01和β=0.02時方程的解．解在MATLAB中，兩個變量x和y可以表示為向量y中的兩個值．同樣，導(dǎo)數(shù)是向量yp中的兩個值.當α=0.01和β=0.02時，方程組（9-1）可以表示為：

yp(1)=(1-alpha*y(2))*y(1)

yp(2)=(-1+beta*y(1))*y(2)MATLAB中的函數(shù)文件lotka.m來表示Lotka-Volterra方程，typelotkafunctionyp=lotka(t,y)%LOTKALotka-Volterrapredator-preymodel.%Copyright1984-2014TheMathWorks,Inc.yp=diag([1-.01*y(2),-1+.02*y(1)])*y;首先使用ode23在區(qū)間0<t<15中求解lotka中定義的微分方程．t0=0;tfinal=15;y0=[20;20];[t,y]=ode23(@lotka,[t0tfinal],y0);繪制兩個種群數(shù)量對時間的圖，見圖9-3a．plot(t,y(:,1),'--',t,y(:,2))title('捕食者-獵物模型')xlabel('時間')ylabel('數(shù)量')legend('獵物','捕食者','Location','North')再使用

ode45

求解該方程組，得到相軌線圖9-3b.[T,Y]=ode45(@lotka,[t0tfinal],y0);plot(y(:,1),y(:,2),'-',Y(:,1),Y(:,2),'--');title('相軌線')legend('ode23','ode45')

圖9-3捕食者-獵物模型a)捕食者—獵物模型b）相軌線

圖9.3a中實線表示的是微分方程的解，虛線表示的是解曲線的導(dǎo)數(shù)所描繪的曲線.可以看出，獵物數(shù)量減少時，捕食者的數(shù)量也在減少，捕食者的種群數(shù)量會隨著獵物種群數(shù)量的增加而不斷增加．獵物的種群數(shù)量增加時，捕食者數(shù)量也在增加，但是當捕食者達到一定的程度后，獵物又在不斷減少.即獵物種群數(shù)量增加→捕食者種群數(shù)量增加→獵物種群數(shù)量減少→捕食者種群數(shù)量減少→獵物種群數(shù)量增加→再次循環(huán)．這種變化趨勢反映了生態(tài)系統(tǒng)中普遍存在的負反饋調(diào)節(jié)．

圖9.3b中，使用不同的數(shù)值方法求解微分方程會產(chǎn)生略微不同的答案.可以看出利用ode45函數(shù)得到的圖形是平滑的，要比ode23函數(shù)得到的結(jié)果精度更高．

例9.5（羅倫茲吸引子與“蝴蝶效應(yīng)”）吸引子在1963年由美國麻省理工學(xué)院的氣象學(xué)家羅倫茲（E.N.Lorenz）發(fā)現(xiàn).羅倫茲教授在研究天氣的不可預(yù)測性時，通過簡化方程，獲得了具有三個自由度的系統(tǒng)．在計算機上用他所建立的微分方程模擬氣候變化，意外地發(fā)現(xiàn)，初始條件的極微小差別可以引起模擬結(jié)果的巨大變化，這表明天氣過程以及描述它們的非線性方程時如此的不穩(wěn)定，以至巴西熱帶雨林的一只蝴蝶偶然拍動一下翅膀，幾星期后可以在美國德克薩斯州引起一場龍卷風，這就是“蝴蝶效應(yīng)”．

羅倫茲根據(jù)牛頓定律建立的溫度、壓強和風速之間的微分方程組為給定初值條件求當

時，微分方程組的解．解當

時，得到微分方程組：將三個方程的右端函數(shù)寫成向量形式：新建函數(shù)文件flo.m，命令如下：functionz=flo(t,y)

A=[-5./30y(2);0-88;-y(2)25-1];

z=A*y;在命令行窗口調(diào)用flo.m，>>[t,y]=ode23(@flo,[0,80],[002.2204e-16]);

>>u=y(:,1);v=y(:,2);w=y(:,3);plot3(u,v,w)結(jié)果如圖9-4所示．圖9-4羅倫茲吸引子10.1線性規(guī)劃

線性規(guī)劃（Linearprogramming,簡稱LP），是運籌學(xué)中研究較早、發(fā)展較快、應(yīng)用廣泛、方法較成熟的一個重要分支，它是輔助人們進行科學(xué)管理的一種數(shù)學(xué)方法．研究線性約束條件下線性目標函數(shù)的極值問題的數(shù)學(xué)理論和方法．MATLAB解決線性規(guī)劃問題的標準形式為：其中c、x、b、beq、lb、ub均為列向量；A、Aeq為矩陣，求

z

的最大值就是求

–z

的最小值．

在MATLAB中利用函數(shù)linprog來解決這類問題.函數(shù)linprog的調(diào)用格式如下：

X=linprog(f,A,b)

[X,fval,exitflag,output,lamnda]=linprog(f,A,b,Aeq,Beq,LB,UB,X0,options)

這里，X是問題的解向量，f是由目標函數(shù)的系數(shù)構(gòu)成的向量，A是一個矩陣，b是一個向量，A，b和變量x={x1,x2,…,xn}一起，表示了線性規(guī)劃中不等式約束條件，A，b是系數(shù)矩陣和右端向量.Aeq和Beq表示了線性規(guī)劃中等式約束條件中的系數(shù)矩陣和右端向量.LB和UB是約束變量的下界和上界向量，X0是給定的變量的初始值，options為控制規(guī)劃過程的參數(shù)系列.返回值中fval是優(yōu)化結(jié)束后得到的目標函數(shù)值.exitflag=0表示優(yōu)化結(jié)果已經(jīng)超過了函數(shù)的估計值或者已聲明的最大迭代次數(shù)；exitflag>0表示優(yōu)化過程中變量收斂于解X，exitflag<0表示不收斂.output有3個分量，iterations表示優(yōu)化過程的迭代次數(shù)，cgiterations表示PCG迭代次數(shù)，algorithm表示優(yōu)化所采用的運算規(guī)則.lambda有4個分量，ineqlin是線性不等式約束條件，eqlin是線性等式約束條件，upper是變量的上界約束條件，lower是變量的下界約束條件.它們的返回值分別表示相應(yīng)的約束條件在約束條件在優(yōu)化過程中是否有效．例10.1求解線性規(guī)劃問題：

解MATLAB命令如下：clearf=-[5,4,6];A=[1,-2,1;3,2,4;3,2,0];b=[20,42,30];LB=[0;0;0];[X,fval]=linprog(f,A,b,[],[],LB)運行結(jié)果：Optimalsolutionfound.X=015.00003.0000fval=-78可知當x1=0，x2=15，x3=3時，得到最小值-78．例10.2求解線性規(guī)劃問題：解先將最大值問題轉(zhuǎn)化為標準形式：MATLAB命令如下：c=[-2,-3,5];A=[-2,5,-1];b=-10;Aeq=[111];beq=7;x0=[000];[x,fval]=linprog(c,A,b,Aeq,beq,x0)運行結(jié)果為：Optimalsolutionfound.x=6.42860.57140fval=-14.5714可知，當x1=6.4286，x2=0.5714，x3=0時，得到最大值z=14.5714．

例10.3某工廠生產(chǎn)A，B兩種產(chǎn)品，所用原料均為甲、乙、丙三種，生產(chǎn)一件產(chǎn)品所需原料和所獲利潤以及庫存原料情況如表10-1所示.

在該廠只有表中所列庫存原料的情況下，如何安排A，B兩種產(chǎn)品的生產(chǎn)數(shù)量可以獲得最大利潤？原料甲（公斤）原料乙（公斤）原料丙（公斤）利潤（元）產(chǎn)品A8447000產(chǎn)品B68610000庫存原料量380300220

表10-1利潤以及庫存原料情況表解設(shè)生產(chǎn)A產(chǎn)品x1件，生產(chǎn)B產(chǎn)品x2件，z為所獲利潤，我們將問題歸結(jié)為如下的線性規(guī)劃問題：

轉(zhuǎn)換成最小值問題

接著寫出MATLAB程序如下：clearf=-[7000,10000];A=[8,6;4,8;4,6];b=[380,300,220];[X,fval]=linprog(f,A,b)運行結(jié)果為：Optimalsolutionfound.X=40.000010.0000fval=-380000可知生產(chǎn)A產(chǎn)品40件，B產(chǎn)品10件時可獲得最大利潤380000元．10.2二次規(guī)劃類似于線性規(guī)劃，求解二次規(guī)劃之前需要先將其化為標準形式：其中f、x、b、beq、lb、ub均為列向量；A、Aeq為矩陣；H為二次型（對稱正定矩陣）．在MATLAB中利用函數(shù)quadprog求解，調(diào)用格式如下

[x,fval]=quadprog(H,f,A,b,Aeq,beq,lb,ub)求解約束條件下的二次規(guī)劃例10.4求解二次規(guī)劃：解將二次規(guī)劃寫成標準形式其中MATLAB程序如下：H=[1-1;-12];f=[-2;-6];A=[11;-12;21];b=[2;2;3];[x,fval]=quadprog(H,f,A,b)運行后得到Minimumfoundthatsatisfiestheconstraints.Optimizationcompletedbecausetheobjectivefunctionisnon-decreasinginfeasibledirections,towithinthevalueoftheoptimalitytolerance,andconstraintsaresatisfiedtowithinthevalueoftheconstrainttolerance.<stoppingcriteriadetails>x=0.66671.3333fval=-8.2222可得當x1=0.6667，x2=1.3333時，最小值是z=-8.2222．10.3無約束優(yōu)化10.3.1一元函數(shù)無約束優(yōu)化的最優(yōu)解求解

一元函數(shù)優(yōu)化一般要給定自變量的取值范圍，其標準形式為：

在MATLAB中利用fminbnd函數(shù)求解，調(diào)用格式如下：

[x,fval]=fminbnd(fun,x1,x2)返回一個值x，該值是fun中描述的標量值函數(shù)在區(qū)間x1<x<x2中的局部最小值；

[x,fval]=fminbnd(fun,x1,x2,options)使用options中指定的優(yōu)化選項求最小值.

在實際求解過程中如果需要給出求解的初值，要求能夠從問題本身進行初步分析后得到最優(yōu)解的大致位置．例10.5求sin(x)函數(shù)在0<x<2π范圍內(nèi)的最小值的點．解>>fun=@sin;>>x1=0;>>x2=2*pi;>>x=fminbnd(fun,x1,x2)x=4.7124此值與正確值

x=3π/2

相同>>3*pi/2ans=4.7124

例10.6對邊長為5m的正方形鐵板，在4個角處減去相等的正方形，以制成方形無蓋水槽，問如何剪才能使水槽的容積最大？

解假設(shè)剪去正方形的邊長為x，則水槽的容積為

水槽的容積最大的目標函數(shù)是

，將其轉(zhuǎn)化為標準形式為

首先建立函數(shù)文件：

functionv=myfun1(x)

v=-(5-2*x)^2*x;然后在命令行調(diào)用fminbnd函數(shù)，>>[x,v]=fminbnd('myfun1',0,2.5)x=0.8333v=-9.2593可知，當剪去正方形的邊長為0.8333m時，水槽的容積最大為9.2593m3．10.3.2多元函數(shù)無約束優(yōu)化的最優(yōu)解求解

多元函數(shù)無約束優(yōu)化的標準形式：

，其中

是向量

在MATLAB中利用函數(shù)fminunc和fminsearch求解，fminunc調(diào)用格式如下：

[x,fval]=fminunc(fun,x0)在點x0處開始并嘗試求fun中描述的函數(shù)的局部最小值點x和函數(shù)值fval

[x,fval]=fminunc(fun,x0,options)使用options中指定的優(yōu)化選項求最小值

fminsearch調(diào)用格式與fminunc相同，fminunc為無約束優(yōu)化提供了大型優(yōu)化和中型優(yōu)化算法，fminsearch采用簡單搜索法.當函數(shù)的階數(shù)大于2時，使用fminunc比fminsearch更有效，當函數(shù)高度不連續(xù)時，使用fminsearch更有效．

例10.7求

的最小值點.解為了對函數(shù)有一個直觀認識，先繪制出其三維圖形．