常微分方程求解

常微分方程求解第一頁，共三十八頁，2022年，8月28日本章知識要點數(shù)值計算常微分方程初值問題常微分方程邊值問題MATLAB微分方程求解常微分方程的相關(guān)函數(shù)ode45ode23bvp4c第二頁，共三十八頁，2022年，8月28日微分方程在化工模型中的應用間歇反應器的計算活塞流反應器的計算全混流反應器的動態(tài)模擬定態(tài)一維熱傳導問題逆流壁冷式固定床反應器一維模型固定床反應器的分散模型第三頁，共三十八頁，2022年，8月28日Matlab常微分方程求解問題分類初值問題：定解附加條件在自變量的一端

一般形式為：初值問題的數(shù)值解法一般采用步進法，如Runge-Kutta法邊值問題：在自變量兩端均給定附加條件一般形式：邊值問題可能有解、也可能無解，可能有唯一解、也可能有無數(shù)解邊值問題有3種基本解法迭加法打靶法松弛法第四頁，共三十八頁，2022年，8月28日Matlab求解常微分方程初值問題方法將待求解轉(zhuǎn)化為標準形式，并“翻譯”成Matlab可以理解的語言，即編寫odefile文件選擇合適的解算指令求解問題根據(jù)求解問題的要求，設置解算指令的調(diào)用格式第五頁，共三十八頁，2022年，8月28日Matlab求解初值問題函數(shù)指令含義指令含義解算ode23普通2-3階法解ODEodefileODE文件格式ode45普通4-5階法解ODE選項odeset創(chuàng)建、更改ODE選項的設置ode113普通變階法解ODEodeget讀取ODE選項的設置ode23t解適度剛性ODE輸出odeplotODE的輸出時間序列圖ode15s變階法解剛性ODEodephas2ODE的二維相平面圖ode23s低階法解剛性ODEodephas3ODE的三維相空間圖ode23tb低階法解剛性ODEodeprint在Matlab指令窗顯示結(jié)果第六頁，共三十八頁，2022年，8月28日odefile所謂的odefile實際上是一個Matlab函數(shù)文件，一般作為整個求解程序的一個子函數(shù)，表示ode求解問題Matlab提供了odefile的模板，采用typeodefile命令顯示其詳細內(nèi)容，然后將其復制到腳本編輯窗口，在合適的位置填入所需內(nèi)容一般而言，對于程序通用性要求不高的場合，只需將原有模型寫成標準形式，然后“翻譯”成Matlab語言即可第七頁，共三十八頁，2022年，8月28日odefile的編寫規(guī)定ode文件的最簡單格式必須有一個自變量t和函數(shù)y作為輸入變量，一個y的導函數(shù)作為輸出變量。其中自變量t不論在ode文件中是否使用都必須作為第一輸入變量，y則必須作為第二輸入變量，位置不能顛倒?？梢韵騩de文件中傳遞參數(shù)，數(shù)目不受限制第八頁，共三十八頁，2022年，8月28日odefile的編寫functionf=fun(x,y)f=y-2*x/y;求解初值問題：

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

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

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

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

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

第十四頁，共三十八頁，2022年，8月28日ode45和ode23的比較functionCha5demo1%Comparisonofresultsobtainedfromode45andode23solverformatlongy0=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('ODE45','ODE23','location','Northwest')

disp('ComparativeResultsatx=1:');fprintf('\nODE45\t\t\ty=%.8f\nODE23\t\t\ty=%.8f\nPrecisiveResult...=%.8f\n',y1(end),y2(end),1.7320508)%------------------------------------------------------------------------functionf=fun(x,y)f=y-2*x/y;第十五頁，共三十八頁，2022年，8月28日ode解算指令的選擇(2)2.根據(jù)常微分方程組是否為剛性方程

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

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

第十六頁，共三十八頁，2022年，8月28日ode解算指令的選擇(2)剛性方程的物理意義：常微分方程組所描述的物理化學變化過程中包含了多個子過程，其變化速度相差非常大的數(shù)量級，于是常微分方程組含有快變和慢變分量。常微分方程組數(shù)值積分的穩(wěn)定步長受模值最大的特征值控制，即受快變量分量約束，特征值大則允許步長??；而過程趨于穩(wěn)定的時間又由模值最小的特征值控制，特征值小則積分到穩(wěn)定的時間則長。Matlab提供了不同種類的剛性方程求解指令：ode15sode23sode23tode23tb，可根據(jù)實際情況選用第十七頁，共三十八頁，2022年，8月28日functionCha5demo3figureode23s(@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];剛性常微分方程組求解第十八頁，共三十八頁，2022年，8月28日解算指令的輸出控制[T,Y]=ode45(@fun,TSPAN,Y0,options,P1,P2,…)sol=ode45(@fun,TSPAN,Y0)將解輸出指結(jié)構(gòu)體sol中；SXINT=deval(SOL,XINT)用于計算解在XINT處的值，XINT必須位于區(qū)間[SOL.x(1)SOL.x(end)]內(nèi)。在無輸出變量時，將調(diào)用默認的odeplot輸出解的圖形。除了以odeplot形式輸出外，還可以以odephas2，和odephas3的形式輸出解向量的二維和三維相平面圖。采用以下語句options=odeset(‘outputfcn’,‘odephas2’)可以將輸出方法改變?yōu)槠矫孑敵鰋deprint輸出求解過程每一步的解第十九頁，共三十八頁，2022年，8月28日解算指令的options選項RelTol－相對誤差，它應用于解向量的所有分量。在每一步積分過程中，第i個分量誤差e(i)滿足：e(i)<=max(RelTol*abs(y(i),AbsTol(i))。AbsTol－絕對誤差，若是實數(shù)，則應用于解向量的所有分量，若是向量，則它的每一個元素應用于對應位置解向量元素。OutputFcn－可調(diào)用的輸出函數(shù)名。每一步計算完后，這個函數(shù)將被調(diào)用輸出結(jié)果，可以選擇的值為：odeplot，odephas2，odephas3，odeprint。OutputSel－輸出序列選擇。指定解向量的哪個分量被傳遞給OutputFcn。MaxSetp－步長上界，缺省值為求解區(qū)間的1/10。InitialStep－初始步長，缺省時自動設置。采用odeset改變原有選項的值第二十頁，共三十八頁，2022年，8月28日在間歇反應器中進行液相反應制備產(chǎn)物B，反應網(wǎng)絡如圖5-1所示。反應可在180~260℃的溫度范圍內(nèi)進行，反應物X大量過剩，而C,D和E為副產(chǎn)物。各反應均為一級動力學關(guān)系：r＝－kC，式中

間歇反應器中串聯(lián)－平行復雜反應系統(tǒng)已知參數(shù)：k01＝5.78052×1010，k02＝3.92317×1012，k03＝1.64254×104，k04＝6.264×108，Ea1=124670，Ea2=150386，Ea3=77954，Ea4=111528。初始濃度：CA=1kmol/m3，其余物質(zhì)濃度為0。已知是產(chǎn)物B收率最大的最優(yōu)反應溫度為224.6℃試計算1)在最優(yōu)反應溫度下各組分濃度隨時間的動態(tài)變化；2)最優(yōu)反應時間；3)輸出產(chǎn)物D對反應物濃度A的關(guān)系圖。第二十一頁，共三十八頁，2022年，8月28日間歇反應器中串聯(lián)－平行復雜反應系統(tǒng)數(shù)學模型

第二十二頁，共三十八頁，2022年，8月28日間歇反應器中串聯(lián)－平行復雜反應系統(tǒng)functionCha5demo4T=224.6+273.15;R=8.31434;k0=[5.78052E+103.92317E+121.64254E+46.264E+8];Ea=[12467015038677954111528];%Initialconcentration:C0(i),kmol/m^3C0=[10000];tspan=[01e4];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/m^3)');legend('A','B','C','D')CBmax=max(C(:,2));

%CBmaxyBmax=CBmax/C0(1)%yBmax:index=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];第二十三頁，共三十八頁，2022年，8月28日Matlab求解邊值問題方法：bvp4c函數(shù)把待解的問題轉(zhuǎn)化為標準邊值問題因為邊值問題可以多解，所以需要為期望解指定一個初始猜測解。該猜測解網(wǎng)（Mesh）包括區(qū)間[a,b]內(nèi)的一組網(wǎng)點（Meshpoints）和網(wǎng)點上的解S(x)根據(jù)原微分方程構(gòu)造殘差函數(shù)利用原微分方程和邊界條件，借助迭代不斷產(chǎn)生新的S(x)，使殘差不斷減小，從而獲得滿足精度要求的解第二十四頁，共三十八頁，2022年，8月28日Matlab求解邊值問題的基本指令solinit＝bvpinit（x,v,parameters） 生成bvp4c調(diào)用指令所必須的“解猜測網(wǎng)”sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…） 給出微分方程邊值問題的近似解

sxint＝deval（sol，xint） 計算微分方程積分區(qū)間內(nèi)任何一點的解值第二十五頁，共三十八頁，2022年，8月28日初始解生成函數(shù)：bvpinit()solinit＝bvpinit（x,v,parameters）x指定邊界區(qū)間[a,b]上的初始網(wǎng)絡，通常是等距排列的（1×M）一維數(shù)組。注意：使x(1)=a，x(end)=b；格點要單調(diào)排列。v是對解的初始猜測solinit（可以取別的任意名）是“解猜測網(wǎng)（Mesh）”。它是一個結(jié)構(gòu)體，帶如下兩個域：solinit.x是表示初始網(wǎng)格有序節(jié)點的（1×M）一維數(shù)組，并且solinit.x(1)一定是a，solinit.x(end)一定是b。M不宜取得太大，10數(shù)量級左右即可。solinit.y是表示網(wǎng)點上微分方程解的猜測值的（N×M）二維數(shù)組。solinit.y(:,i)表示節(jié)點solinit.x(i)處的解的猜測值。第二十六頁，共三十八頁，2022年，8月28日sol＝bvp4c（odefun，bcfun，solinit，options，p1，p2，…）輸入?yún)?shù)：odefun是計算導數(shù)的M函數(shù)文件。該函數(shù)的基本形式為：dydx＝odefun(x,y,parameters,p1,p2，…)，在此，自變量x是標量，y，dydx是列向量。其基本形式與初值問題的odefile形式相同bcfun是計算邊界條件下殘數(shù)的M函數(shù)文件。其基本形式為：res＝bcfun(ya,yb,parameters,p1,p2,……)，文件輸入宗量ya,yb是邊界條件列向量,分別代表y在a和b處的值。res是邊界條件滿足時的殘數(shù)列向量。注意：例如odefun函數(shù)的輸入宗量中包含若干“未知”和“已知”參數(shù)，那么不管在邊界條件計算中是否用到，它們都應作為bcfun的輸入宗量。輸入宗量options是用來改變bvp4c算法的控制參數(shù)的。在最基本用法中，它可以缺省，此時一般可以獲得比較滿意的邊值問題解。如需更改可采用bvpset函數(shù)，使用方法同odeset函數(shù)。輸入宗量p1，p2等表示希望向被解微分方程傳遞的已知參數(shù)。如果無須向微分方程傳遞參數(shù)，它們可以缺省。邊值問題求解指令：bvp4c()第二十七頁，共三十八頁，2022年，8月28日輸出參數(shù)：輸出變量sol是一個結(jié)構(gòu)體sol.x是指令bvp4c所采用的網(wǎng)格節(jié)點；sol.y是y(x)在sol.x網(wǎng)點上的近似解值；sol.yp是y'(x)在sol.x網(wǎng)點上的近似解值；sol.parameters是微分方程所包含的未知參數(shù)的近似解值。當被解微分方程包含未知參數(shù)時，該域存在。邊值問題求解指令：bvp4c()第二十八頁，共三十八頁，2022年，8月28日邊值問題的求解原方程組等價于以下標準形式的方程組：functionCha5demo5solinit=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)];求解兩點邊值問題：

令：

邊界條件為：第二十九頁，共三十八頁，2022年，8月28日邊值問題的求解原方程組等價于以下標準形式的方程組：functionCha5demo6solinit=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];求解：

令：

邊界條件為：第三十頁，共三十八頁，2022年，8月28日邊值問題的求解functionCha5demo7c=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)格點解')%------------------------------------------------------------------functiondydx=ODEfun(x,y,c)dydx=[y(2);-c*abs(y(1))];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1);yb(1)+2];求解：

令：

c＝1此程序有什么錯誤？第三十一頁，共三十八頁，2022年，8月28日邊值問題的求解functionCha5demo8lmb=15;solinit=bvpinit(linspace(0,pi,10),[1;1],lmb);opts=bvpset('Stats','on');sol=bvp4c(@ODEfun,@BCfun,solinit,opts);lambda=sol.parametersx=[0:pi/60:pi];y=deval(sol,x);plot(x,y(1,:),'b-',sol.x,sol.y(1,:),'ro')legend('解曲線','初始網(wǎng)格點解')%------------------------------------------------------------------functiondydx=ODEfun(x,y,lmb)q=5;dydx=[y(2);-(lmb-2*q*cos(2*x))*y(1)];%------------------------------------------------------------------functionbc=BCfun(ya,yb,lmb)bc=[ya(1)-1;ya(2);yb(2)];求解：

邊界條件：本例中，微分方程與參數(shù)λ的數(shù)值有關(guān)。一般而言，對于任意的λ值，該問題無解，但對于特殊的λ值（特征值），它存在一個解，這也稱為微分方程的特征值問題。對于此問題，可在bvpinit中提供參數(shù)的猜測值，然后重復求解BVP得到所需的參數(shù)，返回參數(shù)為sol.parameters第三十二頁，共三十八頁，2022年，8月28日邊值問題的求解functionCha5demo9solinit=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;........%--------------------------------------------------------------------------functiondydx=ODEfun(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))];%-------------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-1ya(2)-1ya(3)-1ya(4)+10yb(3)-yb(5)];%-------------------------------------------------------------------------functionv=ex1init(x)v=[1；1；-4.5*x^2+8.91*x+1；-10；-4.5*x^2+9*x+0.91];求解：

邊界條件：初始猜測解為：第三十三頁，共三十八頁，2022年，8月28日邊值問題的求解functionCha5demo11infinity=6;solinit=bvpinit(linspace(0,infinity,5),[001]);options=bvpset('stats','on');sol=bvp4c(@ex5ode,@ex5bc,solinit,options);eta=sol.x;f=sol.y;fprintf('\n');fprintf('Cebeci&Kellerreportf''''(0)=0.92768.\n')fprintf('Valuecomputedhereisf''''(0)=%7.5f.\n',f(3,1))clfresetplot(eta,f(2,:));axis([0infinity01.4]);title('Falkner-Skanequation,positivewallshear,\beta=0.5.')xlabel(‘\eta’),ylabel('df/d\eta'),shg%--------------------------------------------------------------------------functiondfdeta=ex5ode(eta,f)beta=0.5;dfdeta=[f(2)；f(3)；-f(1)*f(3)-beta*(1-f(2)^2)];%--------------------------------------------------------------------------functionres=ex5bc(f0,finf)res=[f0(1)；f0(2)；finf(2)-1];求解：

邊界條件：β＝0.5如果取1，計算結(jié)果如何第三十四頁，共三十八頁，2022年，8月28日常微分方程邊值問題的應用已知在球形氧化鋁催化劑進行環(huán)己烷脫氫反應，催化劑直徑為5mm，操作溫度為700K，在此溫度下，k＝4s-1，D=0.05cm2/s。計算催化劑顆粒內(nèi)環(huán)己烷的濃度分布及催化劑的有效因子。數(shù)學模型：質(zhì)量傳遞方程