MATLAB教程【6】微分方程課件

1、1.4.8 解常微分方程 (ordinary differential equation, ODE) 一、引言微分方程求解的數(shù)值算法有多種，常用的有Euler（歐拉法）、Runge Kutta(龍格-庫(kù)塔法）。Euler法稱一步法，用于一階微分方程。0 x0 x1 x2 xn xn+1當(dāng)給定步長(zhǎng)時(shí)：所以 y1 = y0 + hf (x0,y0) y2 = y1 + hf (x1,y1) yi+1 = yi+ hf (xi,yi)第1頁(yè)，共16頁(yè)。Runge Kutta法 龍格-庫(kù)塔法：實(shí)際上取兩點(diǎn)斜率的平均斜率來(lái)計(jì)算的，其精度高于歐拉算法 。 下面是一個(gè)經(jīng)典 的四階龍格庫(kù)塔公式：第2頁(yè)，共16

2、頁(yè)。二、解題方法1.編寫(xiě)odefile文件格式：function ydot=odefile(t,y) ydot=待求的函數(shù)2.選擇和調(diào)用解常微分方程的指令( 常用的有ode23、ode45)指令格式：T,Y=ode23(F,tspan,y0,options,p1,p2,) F 求解的odefile的文件名 tspan 單調(diào)遞增(減)的積分區(qū)間 y0 初始條件矢量 options 用odeset建立的優(yōu)化選項(xiàng)，如用默認(rèn)值則不必輸入 p1,p2 傳遞給F的參數(shù)， T，YT是輸出的時(shí)間列矢量，矩陣Y每個(gè)列矢量是解的一個(gè) 分量第3頁(yè)，共16頁(yè)。(1)建立ode函數(shù)文件% m-function, g1.

3、m function dy=g1(x,y) dy=3*x.2;例1：解一階微分方程在區(qū)間2,4的數(shù)值解 dy/dx=3x2 y(2)=0.5(2)調(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：解一階微分方程在區(qū)間0 ，2的數(shù)值解 dy/dx=2xcos2y y(0)=pi/4x 的積分區(qū)間y的初始條件第4頁(yè)，共16頁(yè)。例3 ：解二階微分方程解：1.先將方程寫(xiě)成一階微分方程 令y(1)=x, y(

4、2)=dx/dt,2.建立函數(shù)文件yjs.m并存盤(pán) function ydot = yjs( t,y ) ydot=y(2); 4 ;3.調(diào)用函數(shù)文件解方程 T,Y=ode23(yjs,0:0.1:10,2,1);時(shí)間t 的積分區(qū)間初始條件第5頁(yè)，共16頁(yè)。為方便令x1=x，x2=x分別對(duì)x1,x2求一階導(dǎo)數(shù)，整理后寫(xiě)成一階微分方程組形式 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);建立m文件解微

5、分方程2.給定區(qū)間、初始值;求解微分方程t,x=ode23(wf, 0,20,0,0.25)plot(t,x), figure(2),plot(x(:,1),x(:,2)第6頁(yè)，共16頁(yè)。例：研究有空氣阻力時(shí)拋體運(yùn)動(dòng)的特征。比較下面三種情況下的拋體的軌道：沒(méi)有空氣阻力；空氣阻力與速度一次方成正比；以及空氣阻力與速度二次方成正比。質(zhì)點(diǎn)運(yùn)動(dòng)微分方程為空氣阻力的三種情況分別對(duì)應(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,將方程寫(xiě)成一階微分方程組第7頁(yè)，共16頁(yè)。%program ddqxn.mm=1; b=0,

6、0.1,0.1; p=0,0,1; %設(shè)置參數(shù)figureaxis(0 9 0 4); %設(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) %畫(huà)軌道的彗星軌跡end函數(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).(

7、p/2);第8頁(yè)，共16頁(yè)。3.確定求解的條件和要求T,Y=ode23(F,tspan,y0,options,p1,p2,)Odeset 建立和修改優(yōu)化選項(xiàng)語(yǔ)句格式：options=odeset(name1,value1, name2,value2,)指定各個(gè)參數(shù)的取值，對(duì)不指定取值的參數(shù)，取默認(rèn)值options=odeset(odeopts,name1,value1,)使用原來(lái)的優(yōu)化選項(xiàng)，但對(duì)其中部分參數(shù)指定新值options=odeset(oldopts,newopts)合并了兩個(gè)優(yōu)化選項(xiàng)oldopts和newopts，重復(fù)部分取newopts的指定值第9頁(yè)，共16頁(yè)。 odeset Ab

8、sTol: positive scalar or vector 1e-6 絕對(duì)誤差。默認(rèn)1e6 BDF: on | off 在使用ODE15s時(shí)是否使用反向差分公式 RelTol: positive scalar 1e-3 相對(duì)誤差。默認(rèn)1e3 Events: function 設(shè)置事件 InitialStep: positive scalar 解方程時(shí)使用的初始步長(zhǎng) Jacobian: matrix | function 設(shè)定是否計(jì)算雅各比矩陣 JPattern: sparse matrix 若返回稀疏雅各比矩陣，為on Mass: matrix | function 返回質(zhì)量矩陣 Mass

9、Singular: yes | no | maybe 質(zhì)量矩陣是否為奇異矩陣 MaxOrder: 1 | 2 | 3 | 4 | 5 ODE15s解剛性方程的最高階數(shù) MaxStep: positive scalar 步長(zhǎng)上界 第10頁(yè)，共16頁(yè)。NormControl: on | off 解向量范數(shù)的誤差控制OutputSel: vector of integers 選擇輸出解向量哪個(gè)分量 Refine: positive integer 細(xì)化輸出因子 Stats: on | off 顯示計(jì)算成本統(tǒng)計(jì)結(jié)果 Vectorized: on | off 設(shè)定向量形式的ODE函數(shù)文件 OutputF

10、cn: function 輸出方式，其選項(xiàng)有：odeplot 按時(shí)間順序畫(huà)出全部變量的解；odephas2 二維相空間中前兩個(gè)變量的圖形；odephas3 三維相空間中三個(gè)變量的圖形；Odeprint 打印輸出第11頁(yè)，共16頁(yè)。 氣象學(xué)家Lorenz提出一篇論文，名叫一只蝴蝶拍一下翅膀會(huì)不會(huì)在Taxas州引起龍卷風(fēng)？論述某系統(tǒng)如果初期條件差一點(diǎn)點(diǎn)，結(jié)果會(huì)很不穩(wěn)定，他把這種現(xiàn)象戲稱做蝴蝶效應(yīng)。Lorenz為何要寫(xiě)這篇論文呢？ 這故事發(fā)生在1961年的某個(gè)冬天，他如往常一般在辦公室操作氣象電腦。平時(shí)，他只需要將溫度、濕度、壓力等氣象數(shù)據(jù)輸入，電腦就會(huì)依據(jù)三個(gè)內(nèi)建的微分方程式，計(jì)算出下一刻可能的氣

11、象數(shù)據(jù)，因此模擬出氣象變化圖。 這一天，Lorenz想更進(jìn)一步了解某段紀(jì)錄的後續(xù)變化，他把某時(shí)刻的氣象數(shù)據(jù)重新輸入電腦，讓電腦計(jì)算出更多的後續(xù)結(jié)果。當(dāng)時(shí)，電腦處理數(shù)據(jù)資料的數(shù)度不快，在結(jié)果出來(lái)之前，足夠他喝杯咖啡并和友人閑聊一陣。在一小時(shí)後，結(jié)果出來(lái)了，不過(guò)令他目瞪口呆。結(jié)果和原資訊兩相比較，初期數(shù)據(jù)還差不多，越到後期，數(shù)據(jù)差異就越大了，就像是不同的兩筆資訊。而問(wèn)題并不出在電腦，問(wèn)題是他輸入的數(shù)據(jù)差了0.000127，而這些微的差異卻造成天壤之別。所以長(zhǎng)期的準(zhǔn)確預(yù)測(cè)天氣是不可能的。例：求解lorlenz方程第12頁(yè)，共16頁(yè)。例：求解lorlenz方程設(shè)y(1)=x, y(2)=y, y(3)=z第13頁(yè)，共16頁(yè)。%program lor.maxis(10 50 -50 50 -50 50); %設(shè)置坐標(biāo)軸范圍view(3) %設(shè)置觀察三維圖形的視角hold ontitle(Lonrenz Attractor) %圖的標(biāo)題options =odeset(OutputFcn,odephas3) ;%設(shè)置輸出方式為三維空間三變量圖形t,y=ode23(lorfun,0,20,0,0,eps,options);%odefile lorfun.mfunction ydot=l