MATLAB教程【6】微分方程

1、1.4.8 解常微分方程解常微分方程 (ordinary differential equation, ode) 一、引言一、引言v微分方程求解的數(shù)值算法有多種，常用的有微分方程求解的數(shù)值算法有多種，常用的有euler（歐拉法）、（歐拉法）、runge kutta(龍格龍格-庫塔法）。庫塔法）。veuler法稱一步法，用于一階微分方程。法稱一步法，用于一階微分方程。00yxyyxfdxdy)(),(0 x0 x1 x2 xn xn+1v當(dāng)給定步長時(shí)：當(dāng)給定步長時(shí)：hyyxxyydxdyf010101(x)所以所以 y1 = y0 + hf (x0,y0) y2 = y1 + hf (x1,y1

2、) yi+1 = yi+ hf (xi,yi)vrunge kutta法法 ),(),()(21121211hkyxfkyxfkkkhyynnnnnn 龍格龍格-庫塔法：實(shí)際上取庫塔法：實(shí)際上取兩點(diǎn)斜率的平均斜率來計(jì)算兩點(diǎn)斜率的平均斜率來計(jì)算的，其精度高于歐拉算法的，其精度高于歐拉算法 。 )(4);););()(23121432113nnnnnnnnnnhkh,yxfkk2hy2hxfkk2hy2hxfk,yxfkk2k2kk6hyy,(,(下面是一個(gè)經(jīng)典下面是一個(gè)經(jīng)典 的四階龍格庫塔公式：的四階龍格庫塔公式：二、解題方法二、解題方法1.編寫編寫odefile文件格式：文件格式：functi

3、on ydot=odefile(t,y) ydot=待求的函數(shù)待求的函數(shù)2.選擇和調(diào)用解常微分方程的指令選擇和調(diào)用解常微分方程的指令( 常用的有常用的有ode23、ode45)指令格式：指令格式：t,y=ode23(f,tspan,y0,options,p1,p2,) f 求解的求解的odefile的文件名的文件名 tspan 單調(diào)遞增單調(diào)遞增(減減)的積分區(qū)間的積分區(qū)間 y0 初始條件矢量初始條件矢量 options 用用odeset建立的優(yōu)化選項(xiàng)，如用默認(rèn)值則不必輸入建立的優(yōu)化選項(xiàng)，如用默認(rèn)值則不必輸入 p1,p2 傳遞給傳遞給f的參數(shù)，的參數(shù)， t，yt是輸出的時(shí)間列矢量，矩陣是輸出的時(shí)

4、間列矢量，矩陣y每個(gè)列矢量是解的一個(gè)每個(gè)列矢量是解的一個(gè) 分量分量(1)建立建立ode函數(shù)文件函數(shù)文件% m-function, g1.m function dy=g1(x,y) dy=3*x.2;例例1：解一階微分方程在：解一階微分方程在區(qū)間區(qū)間2,4的數(shù)值解的數(shù)值解 dy/dx=3x2 y(2)=0.5(2)調(diào)用函數(shù)文件解微分方程調(diào)用函數(shù)文件解微分方程x,num_y=ode23(g1,2,4,0.5); % m-function, g3.m function dy=g3(x,y) dy=2*x*cos(y)2; x,num_y=ode23(g3,0,2,pi/4); 例例2：解一階微分方程

5、在區(qū)：解一階微分方程在區(qū)間間0 ，2的數(shù)值解的數(shù)值解 dy/dx=2xcos2y y(0)=pi/4x 的積分區(qū)間的積分區(qū)間y的初始條件的初始條件例例3 ：解二階微分方程：解二階微分方程12)4(0022 vxaadtxd解：解：1.先將方程寫成一階微分方程先將方程寫成一階微分方程 令令y(1)=x, y(2)=dx/dt,4)2()2()1( dtdyydtdy2.建立函數(shù)文件建立函數(shù)文件yjs.m并存盤并存盤 function ydot = yjs( t,y ) ydot=y(2); 4 ;3.調(diào)用函數(shù)文件解方程調(diào)用函數(shù)文件解方程 t,y=ode23(yjs,0:0.1:10,2,1);時(shí)

6、間時(shí)間t 的積的積分區(qū)間分區(qū)間初始條件初始條件為方便令為方便令x1=x，x2=x分別對分別對x1,x2求一求一階導(dǎo)數(shù)，整理后寫成一階微分方程組階導(dǎo)數(shù)，整理后寫成一階微分方程組形式形式 x1=x2 x2=x2(1-x12)-x1例：例：x+(x2-1)x+x=0 (t=0,20; x0=0,x0=0.25)1.建立m文件wf.mfunction xdot=wf(t,x)xdot=x(2); x(2)*(1-x(1)2)-x(1);1.建立建立m文件文件2.解微分方程解微分方程2.給定區(qū)間、初始值;求解微分方程t,x=ode23(wf, 0,20,0,0.25)plot(t,x), figure(

7、2),plot(x(:,1),x(:,2)例：研究有空氣阻力時(shí)拋體運(yùn)動(dòng)的特征。比較下面三種情況例：研究有空氣阻力時(shí)拋體運(yùn)動(dòng)的特征。比較下面三種情況下的拋體的軌道：沒有空氣阻力；空氣阻力與速度一次方成下的拋體的軌道：沒有空氣阻力；空氣阻力與速度一次方成正比；以及空氣阻力與速度二次方成正比。正比；以及空氣阻力與速度二次方成正比。vgr2p22vbmtddm 質(zhì)點(diǎn)運(yùn)動(dòng)微分方程為質(zhì)點(diǎn)運(yùn)動(dòng)微分方程為空氣阻力的三種情況分別對應(yīng)方程中參數(shù)值為空氣阻力的三種情況分別對應(yīng)方程中參數(shù)值為b=0,0.1,0.1,p=0,0,1令令y(1)=x, y(2)=dx/dt, y(3)=y , y(4)=dy/dt,將方程

8、寫成一階微分將方程寫成一階微分方程組方程組(4)d(3)d(2)d(1)dyyyytt222222(4)(2)(4)(4)(4)(2)(2)(2)ppyymbygdtdyyymbydtdy%program ddqxn.mm=1; b=0,0.1,0.1; p=0,0,1; %設(shè)置參數(shù)設(shè)置參數(shù)figureaxis(0 9 0 4); %設(shè)置坐標(biāo)軸的范圍設(shè)置坐標(biāo)軸的范圍hold onfor i=1:3 %解微分方程三次解微分方程三次 t,y=ode45(ddqxnfun,0:0.001:2,0,5,0,8, ,b(i),p(i),m); comet(y(:,1),y(:,3) %畫軌道的彗星軌跡畫

9、軌道的彗星軌跡end函數(shù)文件函數(shù)文件ddqxnfun.m如下：如下：function ydot=ddqxnfun(t,y,flag,b,p,m)ydot=y(2)； -b/m*y(2)*(y(2).2+y(4).2).(p/2); y(4); -9.8-b/m*y(4)*(y(2).2+y(4).2).(p/2);3.確定求解的條件和要求確定求解的條件和要求t,y=ode23(f,tspan,y0,options,p1,p2,)odeset 建立和修改優(yōu)化選項(xiàng)建立和修改優(yōu)化選項(xiàng)語句格式：語句格式：options=odeset(name1,value1, name2,value2,)指定各個(gè)參數(shù)

10、的取值，對不指定取值的參數(shù)，取默認(rèn)值指定各個(gè)參數(shù)的取值，對不指定取值的參數(shù)，取默認(rèn)值options=odeset(odeopts,name1,value1,)使用原來的優(yōu)化選項(xiàng)，但對其中部分參數(shù)指定新值使用原來的優(yōu)化選項(xiàng)，但對其中部分參數(shù)指定新值options=odeset(oldopts,newopts)合并了兩個(gè)優(yōu)化選項(xiàng)合并了兩個(gè)優(yōu)化選項(xiàng)oldopts和和newopts，重復(fù)部分取，重復(fù)部分取newopts的指定值的指定值 odeset abstol: positive scalar or vector 1e-6 絕對誤差。默認(rèn)絕對誤差。默認(rèn)1e6 bdf: on | off 在使用在使用

11、ode15s時(shí)是否使用反向差分公式時(shí)是否使用反向差分公式 reltol: positive scalar 1e-3 相對誤差。默認(rèn)相對誤差。默認(rèn)1e3 events: function 設(shè)置事件設(shè)置事件 initialstep: positive scalar 解方程時(shí)使用的初始步長解方程時(shí)使用的初始步長 jacobian: matrix | function 設(shè)定是否計(jì)算雅各比矩陣設(shè)定是否計(jì)算雅各比矩陣 jpattern: sparse matrix 若返回稀疏雅各比矩陣，為若返回稀疏雅各比矩陣，為on mass: matrix | function 返回質(zhì)量矩陣返回質(zhì)量矩陣 masssin

12、gular: yes | no | maybe 質(zhì)量矩陣是否為奇異矩陣質(zhì)量矩陣是否為奇異矩陣 maxorder: 1 | 2 | 3 | 4 | 5 ode15s解剛性方程的最高階數(shù)解剛性方程的最高階數(shù) maxstep: positive scalar 步長上界步長上界 normcontrol: on | off 解向量范數(shù)的誤差控制解向量范數(shù)的誤差控制outputsel: vector of integers 選擇輸出解向量哪個(gè)分量選擇輸出解向量哪個(gè)分量 refine: positive integer 細(xì)化輸出因子細(xì)化輸出因子 stats: on | off 顯示計(jì)算成本統(tǒng)計(jì)結(jié)果顯示計(jì)算成

13、本統(tǒng)計(jì)結(jié)果 vectorized: on | off 設(shè)定向量形式的設(shè)定向量形式的ode函數(shù)文件函數(shù)文件 outputfcn: function 輸出方式，其選項(xiàng)有：輸出方式，其選項(xiàng)有：odeplot 按時(shí)間順序畫出全部變量的解；按時(shí)間順序畫出全部變量的解；odephas2 二維相空間中前兩個(gè)變量的圖形；二維相空間中前兩個(gè)變量的圖形；odephas3 三維相空間中三個(gè)變量的圖形；三維相空間中三個(gè)變量的圖形；odeprint 打印輸出打印輸出 氣象學(xué)家氣象學(xué)家lorenz提出一篇論文，名叫一只蝴蝶拍一下翅膀會(huì)不會(huì)在提出一篇論文，名叫一只蝴蝶拍一下翅膀會(huì)不會(huì)在taxas州引起龍卷風(fēng)？論述某系統(tǒng)如果

14、初期條件差一點(diǎn)點(diǎn)，結(jié)果會(huì)很不州引起龍卷風(fēng)？論述某系統(tǒng)如果初期條件差一點(diǎn)點(diǎn)，結(jié)果會(huì)很不穩(wěn)定，他把這種現(xiàn)象戲稱做蝴蝶效應(yīng)。穩(wěn)定，他把這種現(xiàn)象戲稱做蝴蝶效應(yīng)。lorenz為何要寫這篇論文呢？為何要寫這篇論文呢？ 這故事發(fā)生在這故事發(fā)生在1961年的某個(gè)冬天，他如往常一般在辦公室操作氣象電年的某個(gè)冬天，他如往常一般在辦公室操作氣象電腦。平時(shí)，他只需要將溫度、濕度、壓力等氣象數(shù)據(jù)輸入，電腦就會(huì)依據(jù)腦。平時(shí)，他只需要將溫度、濕度、壓力等氣象數(shù)據(jù)輸入，電腦就會(huì)依據(jù)三個(gè)內(nèi)建的微分方程式，計(jì)算出下一刻可能的氣象數(shù)據(jù)，因此模擬出氣象三個(gè)內(nèi)建的微分方程式，計(jì)算出下一刻可能的氣象數(shù)據(jù)，因此模擬出氣象變化圖。變化圖。

15、這一天，這一天，lorenz想更進(jìn)一步了解某段紀(jì)錄的後續(xù)變化，他把某時(shí)刻的想更進(jìn)一步了解某段紀(jì)錄的後續(xù)變化，他把某時(shí)刻的氣象數(shù)據(jù)重新輸入電腦，讓電腦計(jì)算出更多的後續(xù)結(jié)果。當(dāng)時(shí)，電腦處理氣象數(shù)據(jù)重新輸入電腦，讓電腦計(jì)算出更多的後續(xù)結(jié)果。當(dāng)時(shí)，電腦處理數(shù)據(jù)資料的數(shù)度不快，在結(jié)果出來之前，足夠他喝杯咖啡并和友人閑聊一數(shù)據(jù)資料的數(shù)度不快，在結(jié)果出來之前，足夠他喝杯咖啡并和友人閑聊一陣。在一小時(shí)後，結(jié)果出來了，不過令他目瞪口呆。結(jié)果和原資訊兩相比陣。在一小時(shí)後，結(jié)果出來了，不過令他目瞪口呆。結(jié)果和原資訊兩相比較，初期數(shù)據(jù)還差不多，越到後期，數(shù)據(jù)差異就越大了，就像是不同的兩較，初期數(shù)據(jù)還差不多，越到後期，

16、數(shù)據(jù)差異就越大了，就像是不同的兩筆資訊。而問題并不出在電腦，問題是他輸入的數(shù)據(jù)差了筆資訊。而問題并不出在電腦，問題是他輸入的數(shù)據(jù)差了0.000127，而這，而這些微的差異卻造成天壤之別。所以長期的準(zhǔn)確預(yù)測天氣是不可能的。些微的差異卻造成天壤之別。所以長期的準(zhǔn)確預(yù)測天氣是不可能的。例：求解lorlenz方程例：求解lorlenz方程zyxydtdzzydtdyyzxdtdx35101038zyxydtdyzydtdyyzxdtdy35(3)1010(2)38(1)設(shè)設(shè)y(1)=x, y(2)=y, y(3)=z%program lor.maxis(10 50 -50 50 -50 50); %設(shè)置坐標(biāo)軸范圍設(shè)置坐標(biāo)軸范圍view(3) %設(shè)置觀察三維圖形的視角設(shè)置觀察三維圖形的視角hold ontitle(lonrenz attractor) %圖的標(biāo)題圖的標(biāo)題options =odeset(outputfcn,odephas3) ;%設(shè)置輸出方式為三維空間三變量圖形設(shè)置輸出方式為三維空間三變量圖形t,y=ode23(lorfun,0,20,0,0,eps,options);%odefi