Matlab-求解化工常微分方程和偏微分方程課件

Matlab求解化工常微分方程和偏微分方程方利國(guó)Matlab求解化工常微分方程和偏微分方程方利國(guó)1Matlab求解化工常微分方程和偏微分方程1、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命令1.2初值問(wèn)題求解1.3邊值問(wèn)題求解1.4加權(quán)問(wèn)題求解（自學(xué)內(nèi)容）2、偏微分方程（組）求解2.1問(wèn)題描述及一維動(dòng)態(tài)PDE方程求解2

.2二維求解Matlab求解化工常微分方程和偏微分方程1、常微分方程（21、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命令

常微分方程：（初值問(wèn)題）常微分方程：（兩點(diǎn)邊值問(wèn)題）

1、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命31、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命令

微分方程組：

1、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命41、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命令

高級(jí)微分方程：

1、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命51、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命令Matlab調(diào)用命令：ODE45:4-5階龍格庫(kù)塔法（非剛性）ODE23:2-3階龍格庫(kù)塔法（非剛性）ODE113:可變D-B-M法（非剛性）ODE15S:基于數(shù)值差分的可變階方法法（剛性）ODE23S、ODE23t、ODE23tb（剛性）

1、常微分方程（組）求解1.1問(wèn)題描述及Matlab調(diào)用命61、常微分方程（組）求解通用調(diào)用格式：[x,y]=ode***(@odefun,xspan,y0,<option,p1,p2..>)X:自變量向量，在實(shí)際調(diào)用時(shí)取名不一定要用x，也可以用其他名稱(chēng)，只要前后一致即可。Y:應(yīng)變量向量，在實(shí)際調(diào)用時(shí)取名不一定要用Y，也可以用其他名稱(chēng)，只要前后一致即可。***:根據(jù)不同的問(wèn)題調(diào)用不同格式，如45，23s@odefun:自定義函數(shù)的函數(shù)名，該函數(shù)為Xspan：自變量的積分限，[xa,xb],也可以是離散點(diǎn)，[x0,x1,x2,…xf]y0：應(yīng)變量向量的初值<>:可以沒(méi)有該選項(xiàng)，如有，具體應(yīng)用見(jiàn)下面的實(shí)際例子

1、常微分方程（組）求解通用調(diào)用格式：71.2初值問(wèn)題求解例1：已知某高溫物體其溫降過(guò)程符合以下規(guī)律，其中溫度T的單位為K，時(shí)間

的單位為分鐘，零時(shí)刻高溫物體的溫度為2000K，以1分鐘作為時(shí)間步長(zhǎng)，請(qǐng)計(jì)算零時(shí)刻以后每隔1分鐘至170分鐘的溫度。單個(gè)微分方程1.2初值問(wèn)題求解例1：已知某高溫物體其溫降過(guò)程符合以下8functionxODEs%鐵球從2000K降溫曲線,在7.0版本上調(diào)試通過(guò)%由華南理工大學(xué)方利國(guó)編寫(xiě)，2012年2月29日%歡迎讀者調(diào)用，如有問(wèn)題請(qǐng)告知lgfang@clearall；clcy0=[2000];[x1,y1]=ode45(@f,[0:1:170],y0);%0到170分鐘，每分鐘一個(gè)計(jì)算點(diǎn)[x2,y2]=ode23(@f,[0:1:170],y0);plot(x1,y1,'r-')xlabel('時(shí)間，M')ylabel('溫度，K')holdondisp('Resultsbyusingode45()：')disp('xy(1)')disp([x1y1])disp('Resultsbyusingode23()：')disp('xy(2)')disp([x2y2])plot(x2,y2,'k:')%------------------------------------------------------------------functiondy=f(x,y)%定義降溫速率的微分方程%dy=0.04*y(1)-100;dy=-0.04*exp(0.001*(y(1)-300))*(-300+y(1));functionxODEs9Resultsbyusingode45()：

xy(1)1.0e+003*Columns1through1300.01000.02000.03000.04000.05000.06000.07000.08000.09000.10000.11000.12002.00000.87880.61330.48800.41810.37620.34980.33280.32180.31450.30970.30650.3043Columns14through180.13000.14000.15000.16000.17000.30290.30190.30130.30090.3006Resultsbyusingode23()：

xy(2)1.0e+003*Columns1through1300.01000.02000.03000.04000.05000.06000.07000.08000.09000.10000.11000.12002.00000.87790.61240.48730.41760.37560.34930.33230.32130.31400.30920.30610.3041Columns14through180.13000.14000.15000.16000.17000.30270.30180.30120.30080.3005Resultsbyusingode45()：10Matlab-求解化工常微分方程和偏微分方程課件111.2初值問(wèn)題求解該問(wèn)題相當(dāng)與一個(gè)自變量，兩個(gè)應(yīng)變量問(wèn)題，已知初值及微分表達(dá)式，可以利用ODE45求解。微分方程組求解1.2初值問(wèn)題求解該問(wèn)題相當(dāng)與一個(gè)自變量，兩個(gè)應(yīng)變量問(wèn)題12程序代碼functionuvDEs%微生物消亡問(wèn)題計(jì)算,在7.0版本上調(diào)試通過(guò)%由華南理工大學(xué)方利國(guó)編寫(xiě)，2012年3月12日%歡迎讀者調(diào)用，如有問(wèn)題請(qǐng)告知lgfang@clearall;Clcy0=[1.61.2];[x1,y1]=ode45(@f,[0:0.1:10],y0);%0到3分鐘，每0.1分鐘一個(gè)計(jì)算點(diǎn)u=y1(:,1);v=y1(:,2);plot(x1,u,'r-')xlabel('時(shí)間，M')ylabel('微生物濃度')holdonplot(x1,v,'k:')disp('Resultsbyusingode45()：')disp('xuv')disp([x1y1])%------------------------------------------------------------------functiondy=f(x,y)%定義降溫速率的微分方程f1=0.09*y(1)*(1-y(1)/20)-0.45*y(1)*y(2);f2=0.06*y(2)*(1-y(2)/15)-0.001*y(1)*y(2);dy=[f1;f2];程序代碼functionuvDEs131.2初值問(wèn)題求解

例3：當(dāng)X較大時(shí)，兩種方法計(jì)算結(jié)果有較大不同，為什么？單個(gè)微分方程有零點(diǎn)問(wèn)題？1.2初值問(wèn)題求解例3：當(dāng)X較大時(shí)，兩種方法計(jì)算結(jié)果有14functionL43ODEs%在7.0版本上調(diào)試通過(guò)%由華南理工大學(xué)方利國(guó)編寫(xiě)，2012年2月29日%歡迎讀者調(diào)用，如有問(wèn)題請(qǐng)告知lgfang@clearallclcy0=[1];[x1,y1]=ode45(@f,[0:0.05:10],y0);%0到10，每0.05間隔一個(gè)計(jì)算點(diǎn)[x2,y2]=ode23(@f,[0:0.05:10],y0);%0到10，每0.05間隔一個(gè)計(jì)算點(diǎn)plot(x1,y1,'r-')xlabel('x')ylabel('y')holdondisp('Resultsbyusingode45()：')disp('xy(1)')disp([x1y1])disp('Resultsbyusingode23()：')disp('xy(2)')disp([x2y2])plot(x2,y2,'b--')%------------------------------------------------------------------functiondy=f(x,y)%定義微分方程dy=y^2*cos(x);functionL43ODEs15計(jì)算值xy(1)Columns1through1300.05000.10000.15000.20000.25000.30000.35000.40000.45000.50000.55000.60001.00001.05261.11091.17571.24791.32871.41951.52181.63781.76981.92102.09512.2970Columns14through260.65000.70000.75000.80000.85000.90000.95001.00001.05001.10001.15001.20001.25002.53292.81073.14113.53814.02064.61525.35986.30787.54359.192811.461414.728319.5991Columns27through291.30001.35001.400027.428341.203068.6630計(jì)算值xy(1)16高階微分方程求解求解思路：將高階微分方程通過(guò)變量轉(zhuǎn)換，轉(zhuǎn)變成一級(jí)微分方程組進(jìn)行求解。例4：高階微分方程求解求解思路：將高階微分方程通過(guò)變量轉(zhuǎn)換，轉(zhuǎn)變成17高階微分方程求解程序及解Columns1through1300.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.200000.10020.20150.30520.41290.52630.64760.77890.92301.08271.26141.46271.69081.00001.00531.02271.05441.10261.16981.25881.37231.51361.68611.89352.13982.4296Columns14through211.30001.40001.50001.60001.70001.80001.90002.00001.95022.24612.58412.97053.41233.91714.49365.15132.76773.15953.61104.12864.71965.39186.15417.0159高階微分方程求解程序及解Columns1through18剛性方程求解有些微分組的系數(shù)變化很大，這時(shí)用ODE45就很難收斂求解，這時(shí)可用專(zhuān)門(mén)解決此類(lèi)微分方程的ODE23S來(lái)求解，需要注意的是在解的圖像繪制時(shí)，也需要考慮數(shù)值的波動(dòng)幅度很大，需要引入對(duì)數(shù)坐標(biāo)。例5：

dy1/dx=-0.03*y1+1e4*y2*y3dy2/dx=0.03*y1-2e4*y2*y3-5e7*y2^2dy3/dx=5e7*y2^2y1(0)=1,y2(0)=0,y3(0)=0剛性方程求解有些微分組的系數(shù)變化很大，這時(shí)用ODE45就很難19functiongangxinDEs%剛性問(wèn)題計(jì)算,在7.0版本上調(diào)試通過(guò)%由華南理工大學(xué)方利國(guó)編寫(xiě)，2012年3月13日%歡迎讀者調(diào)用，如有問(wèn)題請(qǐng)告知lgfang@%dy1=-0.03*y1+1e4*y2*y3%dy2=0.03*y1-2e4*y2*y3-5e7*y2^2%dy3=5e7*y2^2clearallclcfigurexspan=[06*logspace(-6,6)];y0=[100];[x1,y1]=ode15s(@f,xspan,y0);%用ode45計(jì)算剛性方程，可能有問(wèn)題u=y1(:,1);v=1e4*y1(:,2);w=y1(:,3);semilogx(x1,u,'r-','linewidth',2)xlabel('x')ylabel('1e4*v')holdonsemilogx(x1,v,'k:','linewidth',2)holdonsemilogx(x1,w,'g-','linewidth',2)gridaxis([10^(-10)10^10-0.21.2])legend('u','v','w')disp('Resultsbyusingode45()：')disp('xuvw')formatlongdisp([x1y1])%------------------------------------------------------------------functiondy=f(x,y)%定義微分方程f1=-0.03*y(1)+1e4*y(2)*y(3);f2=0.03*y(1)-2e4*y(2)*y(3)-5e7*y(2)^2;f3=5e7*y(2)^2;dy=[f1;f2;f3];functiongangxinDEsdisp('x201.3邊值問(wèn)題求解邊值問(wèn)題相對(duì)于初值問(wèn)題而言，多了一個(gè)端點(diǎn)的約束，如果在高階或微分方程組中端點(diǎn)約束過(guò)多，微分方程組可能無(wú)解，端點(diǎn)約束有一定限制?？梢酝ㄟ^(guò)建立離散的方程組，再利用ODE45進(jìn)行求解，但可以利用MATLAB的專(zhuān)用工具求解最好。下面介紹ODE-BVPs的求解器，主要有bvpinit,bvp4c,deval,solinit等。通過(guò)實(shí)際例子介紹這些內(nèi)部函數(shù)的功能。1.3邊值問(wèn)題求解邊值問(wèn)題相對(duì)于初值問(wèn)題而言，多了一個(gè)端點(diǎn)211.3邊值問(wèn)題求解solinit=bvpinit（x,yinit):產(chǎn)生在初始網(wǎng)格上的初始解，以便bvp4c調(diào)用，其中x為自變量網(wǎng)格，yinit為對(duì)應(yīng)函數(shù)的初值。sol=bvp4c(@odefun,@BCfun,solinit,<option,p1,p2..>)@Bcfun:為定義邊界條件方程，Bcfun(ya,yb),其中ya、yb分別表示左右邊界。其他符號(hào)意義同上。deval(sol,xint):計(jì)算任意點(diǎn)處的函數(shù)值。1.3邊值問(wèn)題求解solinit=bvpinit（x,yi221.3邊值問(wèn)題求解

例6：1.3邊值問(wèn)題求解例6：23程序functionBVP4c1%求解兩點(diǎn)邊值問(wèn)題的示例在7.0版本上調(diào)試通過(guò)%由華南理工大學(xué)方利國(guó)編寫(xiě)，2012年3月13日%歡迎讀者調(diào)用，如有問(wèn)題請(qǐng)告知lgfang@clearallclca=0;b=10;solinit=bvpinit(linspace(a,b,101),[00]);sol=bvp4c(@ODEfun,@BCfun,solinit);formatshortx=[0:0.1:10];y=deval(sol,x);y1=y(1,:)y2=y(2,:)plot(x,y1,'r-')xlabel('x')ylabel('y')holdongridplot(x,y2,'k:')legend('y','dy')disp('Resultsbyusingbvp4c：')disp('xydy')disp([x1y1])%------------------------------------------------------------------functiondy=ODEfun(x,y)f1=y(2);f2=0.05*(1+x^2)*y(1)+2;dy=[f1;f2];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-40;yb(1)-80];Resultsbyusingbvp4c：

xydyColumns1through1300.10000.20000.30000.40000.50000.60000.70000.80000.90001.00001.10001.200039.999838.172236.384034.634732.924431.253129.621428.029826.479124.970223.503822.081020.7025-18.4739-18.0779-17.6872-17.2985-16.9088-16.5158-16.1176-15.7125-15.2996-14.8781-14.4476-14.0080-13.5597程序functionBVP4c1disp('Results24計(jì)算結(jié)果計(jì)算結(jié)果251.3邊值問(wèn)題求解

例7：1.3邊值問(wèn)題求解例7：26主要程序

functiondy=ODEfun(x,y)f1=y(2);f2=4*y(1);dy=[f1;f2];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-1;yb(1)-exp(4)];主要程序

functiondy=ODEfun(x,y)27functionBVP4c2%求解兩點(diǎn)邊值問(wèn)題的示例在7.0版本上調(diào)試通過(guò)%由華南理工大學(xué)方利國(guó)編寫(xiě)，2012年3月13日%歡迎讀者調(diào)用，如有問(wèn)題請(qǐng)告知lgfang@clearallclca=0;b=2solinit=bvpinit(linspace(a,b,21),[00]);sol=bvp4c(@ODEfun,@BCfun,solinit);formatshortx=[0:0.1:2];y=deval(sol,x);y1=y(1,:);y2=y(2,:);y3=exp(2*x);plot(x,y1,'r-','linewidth',2)xlabel('x')ylabel('y')holdongridplot(x,y2,'k:','linewidth',2)holdonplot(x,y3,'b:','linewidth',2)legend('y','dy','real-y')disp('Resultsbyusingbvp4c：')disp('xydy')disp([x;y1;y2])%------------------------------------------------------------------functiondy=ODEfun(x,y)f1=y(2);f2=4*y(1);dy=[f1;f2];%------------------------------------------------------------------functionbc=BCfun(ya,yb)bc=[ya(1)-1;yb(1)-exp(4)];functionBVP4c2282、偏微分方程（組）求解

2.1問(wèn)題描述2

.2一維動(dòng)態(tài)PDE方程求解

2

.3二維穩(wěn)態(tài)PDE方程求解(自學(xué)）

2、偏微分方程（組）求解2.1問(wèn)題描述29問(wèn)題描述當(dāng)A，B，C為常數(shù)時(shí)，稱(chēng)為擬線性偏微分方程，當(dāng)A，B，C滿(mǎn)足不同條件時(shí)，分為三種不同的類(lèi)型：不同類(lèi)型的方程，MATLAB求解時(shí)，采用不同或相同的內(nèi)置函數(shù)或工具箱問(wèn)題描述當(dāng)A，B，C為常數(shù)時(shí)，稱(chēng)為擬線性偏微分方程，當(dāng)A，B30PDE求解數(shù)學(xué)模型通式m=0，表示平板，m=1表示圓柱，m=2表示球形，f項(xiàng)表示通量項(xiàng)，,s項(xiàng)表示源項(xiàng)，c項(xiàng)為對(duì)角陣，元素必須大于等于0才可以求解。PDE求解數(shù)學(xué)模型通式m=0，表示平板，m=1表示圓柱，31調(diào)用方法sol=pdepe(m,@pdefun,,@iCfun@BCfun,xspan,tspan,<option,p1,p2..>)@Bcfun:為定義邊界條件方程@iCfun:為定義初始條件方程。其他符號(hào)意義同上，注意邊界條件必須寫(xiě)成以上形式：調(diào)用方法sol=pdepe(m,@pdefun,,@iCf32應(yīng)用策略Pdepe內(nèi)部函數(shù)具體應(yīng)用時(shí)，需將實(shí)際的偏微分方程對(duì)照標(biāo)準(zhǔn)模型，確定c、f、s函數(shù)的具體形式及邊界條件p、q的具體形式及m值。蒸汽入口流體入口，u，t0冷凝液出口流體出口，u，tL應(yīng)用策略Pdepe內(nèi)部函數(shù)具體應(yīng)用時(shí)，需將實(shí)際的偏微分方程33套管動(dòng)態(tài)傳熱問(wèn)題原方程：代入具體數(shù)據(jù)，并對(duì)長(zhǎng)度歸一化無(wú)量綱處理后，得到以下方程：套管動(dòng)態(tài)傳熱問(wèn)題原方程：代入具體數(shù)據(jù)，并對(duì)長(zhǎng)度歸一化無(wú)量綱處34對(duì)照標(biāo)準(zhǔn)模型，T就是標(biāo)準(zhǔn)模型中的函數(shù)u，m=0,

a=0,b=1,t0=0,tf=1標(biāo)準(zhǔn)模型：初始條件：零時(shí)刻所有位置溫度為30，即u0=30

邊界條件：已知在零位置處任意時(shí)間溫度為30，在1位置處，偏導(dǎo)為0pa=u-30,qa=0;pb=0,qb=1

對(duì)照標(biāo)準(zhǔn)模型，T就是標(biāo)準(zhǔn)模型中的函數(shù)u，m=0,a=0,b35程序清單functionl51clcclearallglobaluaubalpha=1.0;%cm/sua=30;ub=30;m=0;a=0;