MATLAB課件(常微分方程數(shù)值解)

1、常微分方程數(shù)值解化工學(xué)院軟件應(yīng)用教科組化工學(xué)院軟件應(yīng)用教科組2006-10本章知識(shí)要點(diǎn) 數(shù)值計(jì)算 常微分方程初值問題 常微分方程邊值問題 MATLAB 微分方程求解常微分方程的相關(guān)函數(shù) ode45 ode23 bvp4c微分方程在化工模型中的應(yīng)用間歇反應(yīng)器的計(jì)算活塞流反應(yīng)器的計(jì)算全混流反應(yīng)器的動(dòng)態(tài)模擬定態(tài)一維熱傳導(dǎo)問題逆流壁冷式固定床反應(yīng)器一維模型固定床反應(yīng)器的分散模型Matlab常微分方程求解問題分類初值問題： 定解附加條件在自變量的一端 一般形式為： 初值問題的數(shù)值解法一般采用步進(jìn)法，如Runge-Kutta法邊值問題：在自變量?jī)啥司o定附加條件一般形式：邊值問題可能有解、也可能無(wú)解，可能

2、有唯一解、也可能有無(wú)數(shù)解邊值問題有3種基本解法 迭加法 打靶法 松弛法0)(),(yayyxfy21)(,)(),(ybyyayyxfyMatlab求解常微分方程初值問題方法將待求解轉(zhuǎn)化為標(biāo)準(zhǔn)形式，并“翻譯”成Matlab可以理解的語(yǔ)言，即編寫odefile文件選擇合適的解算指令求解問題根據(jù)求解問題的要求，設(shè)置解算指令的調(diào)用格式Matlab求解初值問題函數(shù)指令含義指令含義解算ode23普通2-3階法解ODEodefileODE文件格式ode45普通4-5階法解ODE選項(xiàng)odeset創(chuàng)建、更改ODE選項(xiàng)的設(shè)置ode113普通變階法解ODEodeget讀取ODE選項(xiàng)的設(shè)置ode23t解適度剛性O(shè)D

3、E輸出odeplotODE的輸出時(shí)間序列圖ode15s變階法解剛性O(shè)DEodephas2ODE的二維相平面圖ode23s低階法解剛性O(shè)DEodephas3ODE的三維相空間圖ode23tb低階法解剛性O(shè)DEodeprint在Matlab指令窗顯示結(jié)果odefileZ 所謂的odefile實(shí)際上是一個(gè)Matlab函數(shù)文件，一般作為整個(gè)求解程序的一個(gè)子函數(shù)，表示ode求解問題Z Matlab提供了odefile的模板，采用type odefile命令顯示其詳細(xì)內(nèi)容，然后將其復(fù)制到腳本編輯窗口，在合適的位置填入所需內(nèi)容Z 一般而言，對(duì)于程序通用性要求不高的場(chǎng)合，只需將原有模型寫成標(biāo)準(zhǔn)形式，然后“翻譯

4、”成Matlab語(yǔ)言即可odefile的編寫規(guī)定1.ode文件的最簡(jiǎn)單格式必須有一個(gè)自變量t和函數(shù)y作為輸入變量，一個(gè)y的導(dǎo)函數(shù)作為輸出變量。其中自變量t不論在ode文件中是否使用都必須作為第一輸入變量，y則必須作為第二輸入變量，位置不能顛倒。2.可以向ode文件中傳遞參數(shù)，數(shù)目不受限制odefile的編寫function f=fun(x,y) f=y-2*x/y; 求解初值問題： 1)0(2yyxyy) 10( x自變量在前，因變量在后ode輸入函數(shù)輸出變量為因變量導(dǎo)數(shù)的表達(dá)式0)(),(yayyxfy初值問題： 1)0(2yyyy) 10( xfunction f=fun(x,y) f=y

5、y2; 常微分方程組odefile的編寫1000 , 1)0(, 0)0()1 (300010)1 (0001. 0)1 ()1 (04. 02122142221211xyyyyyyyyyy常微分方程組與單個(gè)常微分方程求解方法相同，只需在編寫odefile時(shí)將整個(gè)方程組作為一個(gè)向量輸出。function f=fun(x,y)dy1dx = 0.04*(1-y(1)-(1-y(2).*y(1)+0.0001*(1-y(2).2; dy2dx = -1e4*dy1dx + 3000*(1-y(2).2;f = dy1dx; dy2dx;高階微分方程高階微分方程odefile的編寫的編寫本例的難度：

6、teytbytayt2cos)() )(2)2cos()(,2cos)(2ttbteetatt求解： y(0)=0，y(0)=1，方程系數(shù)非線性可在odefile中定義方程高階，非標(biāo)準(zhǔn)形式方程變形：令y1y；y2y則原方程等價(jià)于：tebyayyyyt2cos122221function f=fun(t,y)a=-exp(-t)+cos(2*pi*t)*exp(-2*t);b=cos(2*pi*t);f=y(2) -a*y(2)2-b*y(1)+exp(t)*b;解算指令的使用方法解算指令的使用方法調(diào)用格式：1. T,Y=ode45(fun, TSPAN,Y0)2. T,Y=ode45(fun,

7、 TSPAN,Y0,options)3. T,Y= ode45(fun, TSPAN,Y0,options,P1,P2,)4. T,Y,TE,YE,IE= ode45(fun, TSPAN,Y0,options,P1,P2,) 說(shuō)明：v 輸出變量T為返回時(shí)間列向量；解矩陣Y的每一行對(duì)應(yīng)于T的一個(gè)元素，列數(shù)與求解變量數(shù)相等。v fun為函數(shù)句柄，為根據(jù)待求解的ODE方程所編寫的ode文件（odefile）；v TSPANT0 TFINAL是微分系統(tǒng)yF(t,y)的積分區(qū)間；Y0為初始條件v options用于設(shè)置一些可選的參數(shù)值，缺省時(shí)，相對(duì)于第一種調(diào)用格式。options中可以設(shè)置的參數(shù)參見o

8、desetv P1，P2，的作用是傳遞附加參數(shù)P1，P2，到ode文件。當(dāng)options缺省時(shí)，應(yīng)在相應(yīng)位置保留，以便正確傳遞參數(shù)。常微分方程初值問題解算指令比較常微分方程初值問題解算指令比較解算指令算法精度ode45四五階Runge-Kutta法較高ode23二三階Runge-Kutta法低ode113可變階Adams-Bashforth-Moulton法ode15s基于數(shù)值差分的可變階方法（BDFs，Gear）低中ode23s二階改進(jìn)的Rosenbrock法低ode23t使用梯形規(guī)則適中ode23tbTR-BDF2（隱式Runge-Kutta法）低ode解算指令的選擇(1)求解初值問題：

9、1)0(2yyxyy) 10( x比較ode45和ode23的求解精度和速度 1.根據(jù)常微分方程要求的求解精度與速度要求 ode45和ode23的比較function Cha5demo1%Comparison of results obtained from ode45 and ode23 solverformat longy0=1;tic,x1,y1=ode45(fun,0,1,y0);t_ode45=toctic,x2,y2=ode23(fun,0,1,y0);t_ode23=tocplot(x1,y1,b-,x2,y2,m-.),xlabel(x),ylabel(y),legend(OD

10、E45,ODE23,location,Northwest) disp(Comparative Results at x=1:); fprintf(nODE45ttt y=%.8fnODE23ttt y=%.8fnPrecisive Result . =%.8fn,y1(end),y2(end),1.7320508) %-function f=fun(x,y) f=y-2*x/y;ode解算指令的選擇(2)2.根據(jù)常微分方程組是否為剛性方程 如果系數(shù)矩陣A的特征值連乘積小于零，且絕對(duì)值最大和最小的特征值之比（剛性比）很大，則稱此類方程為剛性方程 1)0(, 2)0(10099.9901. 021

11、22211yyyyyyy00)()(yxyxbAyy 剛性方程在化學(xué)工程和自動(dòng)控制領(lǐng)域的模型中比較常見。剛性比：100/0.0110000 ode解算指令的選擇(2) 剛性方程的物理意義：常微分方程組所描述的物理化學(xué)變化過程中包含了多個(gè)子過程，其變化速度相差非常大的數(shù)量級(jí)，于是常微分方程組含有快變和慢變分量。 常微分方程組數(shù)值積分的穩(wěn)定步長(zhǎng)受模值最大的特征值控制，即受快變量分量約束，特征值大則允許步長(zhǎng)??；而過程趨于穩(wěn)定的時(shí)間又由模值最小的特征值控制，特征值小則積分到穩(wěn)定的時(shí)間則長(zhǎng)。 Matlab提供了不同種類的剛性方程求解指令：ode15s ode23s ode23t ode23tb，可根據(jù)實(shí)

12、際情況選用function Cha5demo3figureode23s(fun,0,100,0;1)hold on,ode45(fun,0,100,0;1)%-function f=fun(x,y)dy1dx = 0.04*(1-y(1)-(1-y(2).*y(1)+0.0001*(1-y(2).2; dy2dx = -1e4*dy1dx + 3000*(1-y(2).2;f = dy1dx; dy2dx;剛性常微分方程組求解解算指令的輸出控制1. T,Y= ode45(fun, TSPAN,Y0,options,P1,P2,)2. sol=ode45(fun,TSPAN,Y0)將解輸出指結(jié)構(gòu)

13、體sol中；SXINT = deval(SOL,XINT)用于計(jì)算解在XINT處的值，XINT必須位于區(qū)間SOL.x(1) SOL.x(end)內(nèi)。3. 在無(wú)輸出變量時(shí)，將調(diào)用默認(rèn)的odeplot輸出解的圖形。4. 除了以odeplot形式輸出外，還可以以odephas2，和odephas3的形式輸出解向量的二維和三維相平面圖。5. 采用以下語(yǔ)句options=odeset(outputfcn,odephas2)可以將輸出方法改變?yōu)槠矫孑敵?. odeprint輸出求解過程每一步的解解算指令的options選項(xiàng)1. RelTol相對(duì)誤差，它應(yīng)用于解向量的所有分量。在每一步積分過程中，第i個(gè)分量

14、誤差e(i)滿足：e(i)=max(RelTol*abs(y(i),AbsTol(i)。2. AbsTol絕對(duì)誤差，若是實(shí)數(shù)，則應(yīng)用于解向量的所有分量，若是向量，則它的每一個(gè)元素應(yīng)用于對(duì)應(yīng)位置解向量元素。3. OutputFcn可調(diào)用的輸出函數(shù)名。每一步計(jì)算完后，這個(gè)函數(shù)將被調(diào)用輸出結(jié)果，可以選擇的值為：odeplot，odephas2，odephas3，odeprint。4. OutputSel輸出序列選擇。指定解向量的哪個(gè)分量被傳遞給OutputFcn。5. MaxSetp步長(zhǎng)上界，缺省值為求解區(qū)間的1/10。6. InitialStep初始步長(zhǎng)，缺省時(shí)自動(dòng)設(shè)置。7. 采用odeset改變

15、原有選項(xiàng)的值在間歇反應(yīng)器中進(jìn)行液相反應(yīng)制備產(chǎn)物B，反應(yīng)網(wǎng)絡(luò)如圖5-1所示。反應(yīng)可在180260的溫度范圍內(nèi)進(jìn)行，反應(yīng)物X大量過剩，而C, D和E為副產(chǎn)物。各反應(yīng)均為一級(jí)動(dòng)力學(xué)關(guān)系：rkC，式中 間歇反應(yīng)器中串聯(lián)平行復(fù)雜反應(yīng)系統(tǒng)已知參數(shù)：k015.780521010，k023.923171012，k031.64254104，k046.264108，Ea1=124670，Ea2=150386，Ea3=77954，Ea4=111528。初始濃度：CA=1kmol/m3，其余物質(zhì)濃度為0。已知是產(chǎn)物B收率最大的最優(yōu)反應(yīng)溫度為224.6試計(jì)算1)在最優(yōu)反應(yīng)溫度下各組分濃度隨時(shí)間的動(dòng)態(tài)變化；2)最優(yōu)反應(yīng)時(shí)

16、間；3)輸出產(chǎn)物D對(duì)反應(yīng)物濃度A的關(guān)系圖。 RTEaekk0間歇反應(yīng)器中串聯(lián)平行復(fù)雜反應(yīng)系統(tǒng)數(shù)學(xué)模型 AACkkdtdC)(21BABCkCkdtdC31CACCkCkdtdC42DBDCkCkdtdC53DCECkCkdtdC54間歇反應(yīng)器中串聯(lián)平行復(fù)雜反應(yīng)系統(tǒng)function Cha5demo4T = 224.6 + 273.15; R = 8.31434; k0 = 5.78052E+10 3.92317E+12 1.64254E+4 6.264E+8;Ea = 124670 150386 77954 111528;% Initial concentration:C0(i), kmol/

17、m3C0 = 1 0 0 0 0; tspan = 0 1e4;opt=odeset(reltol,1e-4,outputfcn,odephas2,outputsel,1;4)t,C = ode45(MassEquations, tspan, C0,opt,k0,Ea,R,T)% plotplot(t,C(:,1),r-,t,C(:,2),k:,t,C(:,3),b-.,t,C(:,4),k-);xlabel(Time (s);ylabel(Concentration (kmol/m3);legend(A,B,C,D)CBmax = max(C(:,2); % CBmaxyBmax = CBm

18、ax/C0(1) % yBmax: index = find(C(:,2)=CBmax);t_opt = t(index) % t_opt: the optimum batch time, s% -function dCdt = MassEquations(t,C,k0,Ea,R,T)% Reaction rate constants, 1/sk = k0.*exp(-Ea/(R*T); k(5) = 2.16667E-04; % Reaction rates, kmoles/m3 srA = -(k(1)+k(2)*C(1); rB = k(1)*C(1)-k(3)*C(2);rC = k(

19、2)*C(1)-k(4)*C(3); rD = k(3)*C(2)-k(5)*C(4);rE = k(4)*C(3)+k(5)*C(4);% Mass balancesdCdt = rA; rB; rC; rD; rE; Matlab求解邊值問題方法：bvp4c函數(shù)1. 把待解的問題轉(zhuǎn)化為標(biāo)準(zhǔn)邊值問題 2. 因?yàn)檫呏祮栴}可以多解，所以需要為期望解指定一個(gè)初始猜測(cè)解。該猜測(cè)解網(wǎng)（Mesh）包括區(qū)間a, b內(nèi)的一組網(wǎng)點(diǎn)（Mesh points）和網(wǎng)點(diǎn)上的解S(x) 3. 根據(jù)原微分方程構(gòu)造殘差函數(shù)4. 利用原微分方程和邊界條件，借助迭代不斷產(chǎn)生新的S(x)，使殘差不斷減小，從而獲得滿足精度要求的解

20、 0)(),(),(byaygyxfy)(,()()(xSxfxSxrMatlab求解邊值問題求解邊值問題的基本指令的基本指令solinitbvpinit（x,v,parameters）生成bvp4c調(diào)用指令所必須的“解猜測(cè)網(wǎng)”solbvp4c（odefun，bcfun，solinit，options，p1，p2，）給出微分方程邊值問題的近似解 sxintdeval（sol，xint）計(jì)算微分方程積分區(qū)間內(nèi)任何一點(diǎn)的解值初始解生成函數(shù)：初始解生成函數(shù)：bvpinit()solinitbvpinit（x,v,parameters） x指定邊界區(qū)間a,b上的初始網(wǎng)絡(luò)，通常是等距排列的（1M）一維數(shù)

21、組。注意：使x(1)=a，x(end)=b；格點(diǎn)要單調(diào)排列。 v是對(duì)解的初始猜測(cè) solinit（可以取別的任意名）是“解猜測(cè)網(wǎng)（Mesh）”。它是一個(gè)結(jié)構(gòu)體，帶如下兩個(gè)域：solinit.x是表示初始網(wǎng)格有序節(jié)點(diǎn)的（1M）一維數(shù)組，并且solinit.x(1)一定是a，solinit.x(end)一定是b。M不宜取得太大，10數(shù)量級(jí)左右即可。solinit.y是表示網(wǎng)點(diǎn)上微分方程解的猜測(cè)值的（NM）二維數(shù)組。solinit.y(:,i)表示節(jié)點(diǎn)solinit.x(i)處的解的猜測(cè)值。solbvp4c（odefun，bcfun，solinit，options，p1，p2，）輸入?yún)?shù)： odef

22、un是計(jì)算導(dǎo)數(shù)的M函數(shù)文件。該函數(shù)的基本形式為：dydxodefun(x,y,parameters,p1,p2，)，在此，自變量x是標(biāo)量，y，dydx是列向量。其基本形式與初值問題的odefile形式相同 bcfun是計(jì)算邊界條件下殘數(shù)的M函數(shù)文件。其基本形式為：resbcfun(ya,yb,parameters,p1,p2,)，文件輸入宗量ya,yb是邊界條件列向量,分別代表y在a和b處的值。res是邊界條件滿足時(shí)的殘數(shù)列向量。注意：例如odefun函數(shù)的輸入宗量中包含若干“未知”和“已知”參數(shù)，那么不管在邊界條件計(jì)算中是否用到，它們都應(yīng)作為bcfun的輸入宗量。 輸入宗量options是用

23、來(lái)改變bvp4c算法的控制參數(shù)的。在最基本用法中，它可以缺省，此時(shí)一般可以獲得比較滿意的邊值問題解。如需更改可采用bvpset函數(shù)，使用方法同odeset函數(shù)。 輸入宗量p1，p2等表示希望向被解微分方程傳遞的已知參數(shù)。如果無(wú)須向微分方程傳遞參數(shù)，它們可以缺省。邊值問題求解指令：bvp4c ()輸出參數(shù)：輸出變量sol是一個(gè)結(jié)構(gòu)體usol.x是指令bvp4c所采用的網(wǎng)格節(jié)點(diǎn)；usol.y是y(x)在sol.x網(wǎng)點(diǎn)上的近似解值；usol.yp是y(x)在sol.x網(wǎng)點(diǎn)上的近似解值；usol.parameters是微分方程所包含的未知參數(shù)的近似解值。當(dāng)被解微分方程包含未知參數(shù)時(shí)，該域存在。邊值問題

24、求解指令：bvp4c ()邊值問題的求解邊值問題的求解原方程組等價(jià)于以下標(biāo)準(zhǔn)形式的方程組：function Cha5demo5solinit=bvpinit(linspace(0,1,10),1 0);sol=bvp4c(ODEfun,BCfun,solinit);x=0:0.05:0.5;y=deval(sol,x);yP=0 -0.41286057 -0.729740656. -0.95385538 -1.08743325 -1.13181116;xP=0:0.1:0.5;plot(xP,yP,o,x,y(1,:),r-),legend(Analytical Solution,Numeri

25、cal Solution)% - function dydx=ODEfun(x,y)dydx=y(2);y(1)+10;% -function bc=BCfun(ya,yb)bc=ya(1);yb(1);0) 1 ()0(10yyyy求解兩點(diǎn)邊值問題： ,121yyyy令： 101221yyyy0) 1 (, 0)0(11yy邊界條件為： 邊值問題的求解邊值問題的求解原方程組等價(jià)于以下標(biāo)準(zhǔn)形式的方程組：function Cha5demo6solinit = bvpinit(linspace(0,1,10),0 1);sol = bvp4c(ODEfun,BCfun,solinit);x = 0

26、:0.1:1;y = deval(sol,x);yP=1 1.0743 1.1695 1.2869 1.4284. 1.5965 1.7947 2.0274 2.3004 2.6214 3;xP=0:0.1:1.0;plot(xP,yP,o,x,y(1,:),r-)legend(Analytical Solution,Numerical Solution,location,Northwest)legend boxoff% -function dydx = ODEfun(x,y)dydx = y(2); (1+x2)*y(1)+1;% -function bc = BCfun(ya,yb)bc

27、= ya(1)-1;yb(1)-3;求解： ,121yyyy令： 邊界條件為： 3) 1 (, 1)0(1)1 (2yyyxy1)1 (12221yxyyy03) 1 (01)0(11yy邊值問題的求解邊值問題的求解function Cha5demo7c=1;solinit = bvpinit(linspace(0,4,10),1;1);sol = bvp4c(ODEfun,BCfun,solinit,c);x = 0:0.1:4;y = deval(sol,x);plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro)legend(解曲線,初始網(wǎng)格點(diǎn)解)% -functio

28、n dydx = ODEfun(x,y,c)dydx = y(2); -c*abs(y(1);% -function bc = BCfun(ya,yb)bc = ya(1);yb(1)+2;求解： ,121zzzz令： 0 zcz2)4(, 0)0(zzc1 此程序有什么錯(cuò)誤？邊值問題的求解邊值問題的求解function Cha5demo8lmb=15;solinit = bvpinit(linspace(0,pi,10),1;1,lmb);opts=bvpset(Stats,on);sol = bvp4c(ODEfun,BCfun,solinit,opts);lambda=sol.param

29、etersx = 0:pi/60:pi;y = deval(sol,x);plot(x,y(1,:),b-,sol.x,sol.y(1,:),ro)legend(解曲線,初始網(wǎng)格點(diǎn)解)% -function dydx = ODEfun(x,y,lmb)q=5;dydx = y(2); -(lmb-2*q*cos(2*x)*y(1);% -function bc = BCfun(ya,yb,lmb)bc = ya(1)-1;ya(2);yb(2);求解： 邊界條件： 本例中，微分方程與參數(shù)的數(shù)值有關(guān)。一般而言，對(duì)于任意的值，該問題無(wú)解，但對(duì)于特殊的值（特征值），它存在一個(gè)解，這也稱為微分方程的特

30、征值問題。對(duì)于此問題，可在bvpinit中提供參數(shù)的猜測(cè)值，然后重復(fù)求解BVP得到所需的參數(shù)，返回參數(shù)為sol.parameters 0)2cos2( zxqz15q0)(, 0)0(, 1)0(zzz邊值問題的求解邊值問題的求解function Cha5demo9solinit = bvpinit(linspace(0,1,5),ex1init);options = bvpset(Stats,on,RelTol,1e-5);sol = bvp4c(ODEfun,BCfun,solinit,options);x = sol.x;y = sol.y;.% -function dydx = ODE

31、fun(x,y)dydx = 0.5*y(1)*(y(3) - y(1)/y(2) -0.5*(y(3) - y(1) (0.9 - 1000*(y(3) - y(5) - 0.5*y(3)*(y(3) - y(1)/y(4) 0.5*(y(3) - y(1) 100*(y(3) - y(5) ;%-function bc = BCfun(ya,yb)bc = ya(1) - 1 ya(2) - 1 ya(3) - 1 ya(4) + 10 yb(3) - yb(5);%- function v = ex1init(x)v = 1；1； -4.5*x2+8.91*x+1；-10； -4.5*x

32、2+9*x+0.91 ;求解： 邊界條件： )(100)(5 . 0/)(5 . 0)(10009 . 0()(5 . 0/ )(5 . 0wyyuwzzuwwywwuwvvuwuu) 1 () 1 (,10)0(, 1)0()0()0(ywzwvu初始猜測(cè)解為： 91. 095 . 4)(,10)(, 191. 85 . 4)(, 1)(, 1)(22xxxyxzxxxwxvxu邊值問題的求解邊值問題的求解function Cha5demo11infinity = 6;solinit = bvpinit(linspace(0,infinity,5),0 0 1);options = bvps

33、et(stats,on);sol = bvp4c(ex5ode,ex5bc,solinit,options);eta = sol.x;f = sol.y;fprintf(n);fprintf(Cebeci & Keller report f(0) = 0.92768.n)fprintf(Value computed here is f(0) = %7.5f.n,f(3,1)clf resetplot(eta,f(2,:);axis(0 infinity 0 1.4);title(Falkner-Skan equation, positive wall shear, beta = 0.5.)xlabel(eta),ylabel(df/deta),shg% -function dfdeta = ex5ode(eta,f)beta = 0.5;dfdeta = f(2)； f(3)； -f(1)*f(3) - beta*(1 - f(2)2) ;% -function res = ex5bc(f0,finf)res = f0(1)； f0(2)； finf(2) - 1;求解： 邊界條件： 0) (1 ( 2ffff0.5 1)( , 0)0( , 0)0(fff如果取1，計(jì)算結(jié)果如何常微分方程邊值問題的應(yīng)用已