重要：MATLAB常微分方程(組)數(shù)值解法

重要：MATLAB常微分方程(組)數(shù)值解法第一頁(yè)，共31頁(yè)。知識(shí)要點(diǎn)常微分方程初值問(wèn)題---ode45,0de23常微分方程邊值問(wèn)題---bvp4c第二頁(yè)，共31頁(yè)。微分方程在化工模型中的應(yīng)用間歇反應(yīng)器的計(jì)算活塞流反應(yīng)器的計(jì)算全混流反應(yīng)器的動(dòng)態(tài)模擬定態(tài)一維熱傳導(dǎo)問(wèn)題逆流壁冷式固定床反應(yīng)器一維模型固定床反應(yīng)器的分散模型第三頁(yè)，共31頁(yè)。Matlab常微分方程求解問(wèn)題分類(lèi)初值問(wèn)題：定解附加條件在自變量的一端

一般形式為：初值問(wèn)題的數(shù)值解法一般采用步進(jìn)法，如Runge-Kutta法邊值問(wèn)題：在自變量?jī)啥司o定附加條件一般形式：邊值問(wèn)題可能有解、也可能無(wú)解，可能有唯一解、也可能有無(wú)數(shù)解邊值問(wèn)題有3種基本解法迭加法打靶法松弛法第四頁(yè)，共31頁(yè)。Matlab求解常微分方程初值問(wèn)題方法將待求解轉(zhuǎn)化為標(biāo)準(zhǔn)形式，并“翻譯”成Matlab可以理解的語(yǔ)言，即編寫(xiě)odefile文件選擇合適的解算指令求解問(wèn)題根據(jù)求解問(wèn)題的要求，設(shè)置解算指令的調(diào)用格式第五頁(yè)，共31頁(yè)。Matlab求解初值問(wèn)題函數(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è)DEode15s變階法解剛性O(shè)DEode23s低階法解剛性O(shè)DEode23tb低階法解剛性O(shè)DE第六頁(yè)，共31頁(yè)。odefile所謂的odefile實(shí)際上是一個(gè)Matlab函數(shù)文件，一般作為整個(gè)求解程序的一個(gè)子函數(shù)，表示ode求解問(wèn)題一般而言，對(duì)于程序通用性要求不高的場(chǎng)合，只需將原有模型寫(xiě)成標(biāo)準(zhǔn)形式，然后“翻譯”成Matlab語(yǔ)言即可第七頁(yè)，共31頁(yè)。odefile的編寫(xiě)---常微分方程functionf=fun(x,y)f=y-2*x/y;求解初值問(wèn)題：

自變量在前，因變量在后ode輸入函數(shù)輸出變量為因變量導(dǎo)數(shù)的表達(dá)式初值問(wèn)題：

functionf=fun(x,y)f=y＋y^2;第八頁(yè)，共31頁(yè)。odefile的編寫(xiě)---常微分方程組常微分方程組與單個(gè)常微分方程求解方法相同，只需在編寫(xiě)odefile時(shí)將整個(gè)方程組作為一個(gè)向量輸出。functionf=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];第九頁(yè)，共31頁(yè)。高階微分方程odefile的編寫(xiě)本例的難度：求解：

y(0)=0，y'(0)=1，方程系數(shù)非線(xiàn)性可在odefile中定義方程高階，非標(biāo)準(zhǔn)形式方程變形：令y1＝y(tǒng)；y2＝y(tǒng)’則原方程等價(jià)于：functionf=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];第十頁(yè)，共31頁(yè)。解算指令的使用方法[T,Y]=ode45(@fun,TSPAN,Y0)[T,Y]=ode45(@fun,TSPAN,Y0,options)[T,Y]=ode45(@fun,TSPAN,Y0,options,P1,P2,…)[T,Y,TE,YE,IE]=ode45(@fun,TSPAN,Y0,options,P1,P2,…)輸出變量T為返回時(shí)間列向量；解矩陣Y的每一行對(duì)應(yīng)于T的一個(gè)元素，列數(shù)與求解變量數(shù)相等。@fun為函數(shù)句柄，為根據(jù)待求解的ODE方程所編寫(xiě)的ode文件（odefile）；TSPAN＝[T0TFINAL]是微分系統(tǒng)y'＝F(t,y)的積分區(qū)間；Y0為初始條件options用于設(shè)置一些可選的參數(shù)值，缺省時(shí)，相對(duì)于第一種調(diào)用格式。options中可以設(shè)置的參數(shù)參見(jiàn)odesetP1，P2，…的作用是傳遞附加參數(shù)P1，P2，…到ode文件。當(dāng)options缺省時(shí)，應(yīng)在相應(yīng)位置保留[]，以便正確傳遞參數(shù)。第十一頁(yè)，共31頁(yè)。常微分方程初值問(wèn)題解算指令比較解算指令算法精度ode45四五階Runge-Kutta法較高ode23二三階Runge-Kutta法低ode113可變階Adams-Bashforth-Moulton法ode15s基于數(shù)值差分的可變階方法（BDFs，Gear）低～中ode23s二階改進(jìn)的Rosenbrock法低ode23t使用梯形規(guī)則適中ode23tbTR-BDF2（隱式Runge-Kutta法）低第十二頁(yè)，共31頁(yè)。ode解算指令的選擇(1)求解初值問(wèn)題：

比較ode45和ode23的求解精度和速度

1.根據(jù)常微分方程要求的求解精度與速度要求

第十三頁(yè)，共31頁(yè)。ode45和ode23的比較-1functionxODEclearallclcformatlongy0=1;[x1,y1]=ode45(@f,[0,1],y0);[x2,y2]=ode23(@f,[0,1],y0);plot(x1,y1,'k-',x2,y2,'b--')xlabel('x')ylabel('y')%------------------------------------------------------------------functiondydx=f(x,y)dydx=y-2*x/y;第十四頁(yè)，共31頁(yè)。ode45和ode23的比較-2functionxODEsclearallclcy0=[01];[x1,y1]=ode45(@f,[0:0.2:1],y0);[x2,y2]=ode23(@f,[0:0.2:1],y0);disp('Resultsbyusingode45()：')disp('xy(1)y(2)')disp([x1y1]);disp('Resultsbyusingode23()：')disp('xy(1)y(2)')disp([x2y2]);%------------------------------------------------------------------functiondydx=f(x,y)%DefinesimultaneousODEequationsdydx=[3*y(1)+2*y(2);4*y(1)+y(2)];第十五頁(yè)，共31頁(yè)。ode解算指令的選擇(2)2.根據(jù)常微分方程組是否為剛性方程

如果系數(shù)矩陣A的特征值連乘積小于零，且絕對(duì)值最大和最小的特征值之比（剛性比）很大，則稱(chēng)此類(lèi)方程為剛性方程

剛性方程在化學(xué)工程和自動(dòng)控制領(lǐng)域的模型中比較常見(jiàn)。剛性比：100/0.01＝10000

第十六頁(yè)，共31頁(yè)。ode解算指令的選擇(2)剛性方程的物理意義：常微分方程組所描述的物理化學(xué)變化過(guò)程中包含了多個(gè)子過(guò)程，其變化速度相差非常大的數(shù)量級(jí)，于是常微分方程組含有快變和慢變分量。Matlab提供了不同種類(lèi)的剛性方程求解指令：ode15sode23sode23tode23tb，可根據(jù)實(shí)際情況選用第十七頁(yè)，共31頁(yè)。functiondemo1figureode23s(@fun,[0,100],[0;1])holdon,ode45(@fun,[0,100],[0;1])%--------------------------------------------------------------------------functionf=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];剛性常微分方程組求解第十八頁(yè)，共31頁(yè)。間歇反應(yīng)器中串聯(lián)－平行復(fù)雜反應(yīng)系統(tǒng)（例4-1）functionBatchReactorclearallclcT=224.6+273.15;%Reactortemperature,KelvinR=8.31434; %Gasconstant,kJ/kmolK%Arrheniusconstant,1/sk0=[5.78052E+103.92317E+121.64254E+46.264E+8];%Activationenergy,kJ/kmolEa=[12467015038677954111528];%初始濃度C0(i),kmol/m^3C0=[10000];tspan=[01e4];[t,C]=ode45(@MassEquations,tspan,C0,[],k0,Ea,R,T);%繪圖plot(t,C(:,1),'r-',t,C(:,2),'k:',t,C(:,3),'b-.',t,C(:,4),'k--');xlabel('Time(s)');ylabel('Concentration(kmol/m^3)');legend('A','B','C','D')CBmax=max(C(:,2));%CBmax:themaximumconcentrationofB,kmol/m^3yBmax=CBmax/C0(1)%yBmax:themaximumyieldofBindex=find(C(:,2)==CBmax);t_opt=t(index)%t_opt:theoptimumbatchtime,s%------------------------------------------------------------------functiondCdt=MassEquations(t,C,k0,Ea,R,T)%Reactionrateconstants,1/sk=k0.*exp(-Ea/(R*T));k(5)=2.16667E-04;%Reactionrates,kmoles/m3srA=-(k(1)+k(2))*C(1);rB=k(1)*C(1)-k(3)*C(2);rC=k(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);%MassbalancesdCdt=[rA;rB;rC;rD;rE];

第十九頁(yè)，共31頁(yè)。三個(gè)串聯(lián)的CSTR等溫反應(yīng)器（例4-3）functionIsothermCSTRsclearallclcCA0=1.8;%kmol/m^3CA10=0.4;%kmol/m^3CA20=0.2;%kmol/m^3CA30=0.1;%kmol/m^3k=0.5;%1/mintau=2;stoptime=2.9;%min[t,y]=ode45(@Equations,[0stoptime],[CA10CA20CA30],[],k,CA0,tau);disp('Results:')disp('tCA1CA2CA3')disp([t,y])plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-')legend('CA_1','CA_2','CA_3')xlabel('Time(min)')ylabel('Concentration')%------------------------------------------------------------------functiondydt=Equations(t,y,k,CA0,tau)CA1=y(1);CA2=y(2);CA3=y(3);dCA1dt=(CA0-CA1)/tau-k*CA1;dCA2dt=(CA1-CA2)/tau-k*CA2;dCA3dt=(CA2-CA3)/tau-k*CA3;dydt=[dCA1dt;dCA2dt;dCA3dt];第二十頁(yè)，共31頁(yè)。Matlab求解邊值問(wèn)題方法：bvp4c函數(shù)把待解的問(wèn)題轉(zhuǎn)化為標(biāo)準(zhǔn)邊值問(wèn)題因?yàn)檫呏祮?wèn)題可以多解，所以需要為期望解指定一個(gè)初始猜測(cè)解。該猜測(cè)解網(wǎng)（Mesh）包括區(qū)間[a,b]內(nèi)的一組網(wǎng)點(diǎn)（Meshpoints）和網(wǎng)點(diǎn)上的解S(x)根據(jù)原微分方程構(gòu)造殘差函數(shù)利用原微分方程和邊界條件，借助迭代不斷產(chǎn)生新的S(x)，使殘差不斷減小，從而獲得滿(mǎn)足精度要求的解第二十一頁(yè)，共31頁(yè)。Matlab求解邊值問(wèn)題的基本指令solinit＝bvpinit（x,v,parameters） 生成bvp4c調(diào)用指令所必須的“解猜測(cè)網(wǎng)”sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…） 給出微分方程邊值問(wèn)題的近似解sxint＝deval（sol，xint） 計(jì)算微分方程積分區(qū)間內(nèi)任何一點(diǎn)的解值第二十二頁(yè)，共31頁(yè)。初始解生成函數(shù)：bvpinit()solinit＝bvpinit（x,v,parameters）x指定邊界區(qū)間[a,b]上的初始網(wǎng)絡(luò)，通常是等距排列的（1×M）一維數(shù)組。注意：使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)的（1×M）一維數(shù)組，并且solinit.x(1)一定是a，solinit.x(end)一定是b。M不宜取得太大，10數(shù)量級(jí)左右即可。solinit.y是表示網(wǎng)點(diǎn)上微分方程解的猜測(cè)值的（N×M）二維數(shù)組。solinit.y(:,i)表示節(jié)點(diǎn)solinit.x(i)處的解的猜測(cè)值。第二十三頁(yè)，共31頁(yè)。sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…）輸入?yún)?shù)：odefun是計(jì)算導(dǎo)數(shù)的M函數(shù)文件。該函數(shù)的基本形式為：dydx＝odefun(x,y,parameters,p1,p2，…)，在此，自變量x是標(biāo)量，y，dydx是列向量。其基本形式與初值問(wèn)題的odefile形式相同bcfun是計(jì)算邊界條件下殘數(shù)的M函數(shù)文件。其基本形式為：res＝bcfun(ya,yb,parameters,p1,p2,……)，文件輸入?yún)⒘縴a,yb是邊界條件列向量,分別代表y在a和b處的值。res是邊界條件滿(mǎn)足時(shí)的殘數(shù)列向量。注意：例如odefun函數(shù)的輸入?yún)⒘恐邪舾伞拔粗焙汀耙阎眳?shù)，那么不管在邊界條件計(jì)算中是否用到，它們都應(yīng)作為bcfun的輸入?yún)⒘俊＿呏祮?wèn)題求解指令：bvp4c()第二十四頁(yè)，共31頁(yè)。sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…）輸入?yún)?shù)：輸入?yún)⒘縪ptions是用來(lái)改變bvp4c算法的控制參數(shù)的。在最基本用法中，它可以缺省，此時(shí)一般可以獲得比較滿(mǎn)意的邊值問(wèn)題解。如需更改可采用bvpset函數(shù)，使用方法同odeset函數(shù)。輸入?yún)⒘縫1，p2等表示希望向被解微分方程傳遞的已知參數(shù)。如果無(wú)須向微分方程傳遞參數(shù)，它們可以缺省。邊值問(wèn)題求解指令：bvp4c()第二十五頁(yè)，共31頁(yè)。sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…）輸出參數(shù)：輸出變量sol是一個(gè)結(jié)構(gòu)體sol.x是指令bvp4c所采用的網(wǎng)格節(jié)點(diǎn)；sol.y是y(x)在sol.x網(wǎng)點(diǎn)上的近似解值；sol.yp是y'(x)在sol.x網(wǎng)點(diǎn)上的近似解值；sol.parameters是微分方程所包含的未知參數(shù)的近似解值。邊值問(wèn)題求解指令：bvp4c()第二十六頁(yè)，共31頁(yè)。邊值問(wèn)題的求解-1原方程組等價(jià)于以下標(biāo)準(zhǔn)形式的方程組：functiondemo2solinit=bvpinit(linspace(0,1,10),[10]);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('AnalyticalSolution','NumericalSolution')%-----------------------------------------------------------functiondydx=ODEfun(x,y)dydx=[y(2);y(1)+10];%-----------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1);yb(1)];求解兩點(diǎn)邊值問(wèn)題：

令：

邊界條件為：第二十七頁(yè)，共31頁(yè)。邊值問(wèn)題的求解-2原方程組等價(jià)于以下標(biāo)準(zhǔn)形式的方程組：functiondemo3solinit=bvpinit(linspace(0,1,10),[01]);sol=bvp4c(@ODEfun,@BCfun,solinit);x=[0:0.1:1];y=deval(sol,x);yP=[11.07431.16951.28691.4284...1.59651.79472.02742.30042.62143];xP=[0:0.1:1.0];plot(xP,yP,'o',x,y(1,:),'r-')legend('AnalyticalSolution','NumericalSolution','location','Northwest')legendboxoff%------------------------------------------------------------------functiondydx=ODEfun(x,y)dydx=[y(2);(1+x^2)*y(1)+1];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-1;yb(1)-3];求解：

令