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

1、化工常微分方程和偏微分方程Matlab求解Matlab 求解化工常微分方程和偏微分方程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 二維求解1、常微分方程組求解1.1 問(wèn)題描繪及Matlab調(diào)用命令 常微分方程：初值問(wèn)題 常微分方程：兩點(diǎn)邊值問(wèn)題 1、常微分方程組求解1.1 問(wèn)題描繪及Matlab調(diào)用命令 微分方程組： 1、常微分方程組求解1.1 問(wèn)題描繪及Matlab調(diào)用命令 高級(jí)微分方程： 1、常微分方程組求解1.1 問(wèn)題描繪及Matlab調(diào)

2、用命令 Matlab 調(diào)用命令：ODE45:4-5 階龍格庫(kù)塔法非剛性O(shè)DE23:2-3 階龍格庫(kù)塔法非剛性O(shè)DE113:可變D-B-M法非剛性O(shè)DE15S:基于數(shù)值差分的可變階方法法剛性O(shè)DE23S、 ODE23t、 ODE23tb 剛性 1、常微分方程組求解通用調(diào)用格式： x,y=ode*(odefun,xspan,y0,)X:自變量向量，在實(shí)際調(diào)用時(shí)取名不一定要用x，也可以用其他名稱，只要前后一致即可。Y:應(yīng)變量向量，在實(shí)際調(diào)用時(shí)取名不一定要用Y，也可以用其他名稱，只要前后一致即可。 *:根據(jù)不同的問(wèn)題調(diào)用不同格式，如45，23sodefun: 自定義函數(shù)的函數(shù)名，該函數(shù)為Xspan：自

3、變量的積分限，xa,xb,也可以是離散點(diǎn)，x0,x1,x2,xfy0 ：應(yīng)變量向量的初值: 可以沒(méi)有該選項(xiàng)，如有，詳細(xì)應(yīng)用見(jiàn)下面的實(shí)際例子 1.2 初值問(wèn)題求解例1：某高溫物體其溫降過(guò)程符合以下規(guī)律，其中溫度T的單位為K，時(shí)間 的單位為分鐘，零時(shí)刻高溫物體的溫度為2000K，以1分鐘作為時(shí)間步長(zhǎng)，請(qǐng)計(jì)算零時(shí)刻以后每隔1分鐘至170分鐘的溫度。單個(gè)微分方程function xODEs% 鐵球從2000K降溫曲線,在7.0 版本上調(diào)試通過(guò)% 由華南理工大學(xué)方利國(guó)編寫(xiě)，年2月29日% 歡送讀者調(diào)用，如有問(wèn)題請(qǐng)告知 clear all；clcy0 = 2000;x1,y1 = ode45(f,0:1:

4、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)hold ondisp(Results by using ode45()：)disp( x y(1) )disp(x1 y1)disp(Results by using ode23()：)disp( x y(2) )disp(x2 y2)plot(x2,y2,k:)% -function dy = f(x,y) % 定義降溫速率的微分方程%dy = 0.04*y(1)-100;dy=-0.04*exp(0.001

5、*(y(1)-300)*(-300+y(1);Results by using ode45()： x y(1) 1.0e+003 * Columns 1 through 13 Columns 14 through 18 Results by using ode23()： x y(2) 1.0e+003 * Columns 1 through 13 Columns 14 through 18 1.2 初值問(wèn)題求解該問(wèn)題相當(dāng)與一個(gè)自變量，兩個(gè)應(yīng)變量問(wèn)題，初值及微分表達(dá)式，可以利用ODE45求解。微分方程組求解程序代碼function uvDEs% 微生物消亡問(wèn)題計(jì)算,在7.0 版本上調(diào)試通過(guò)% 由

6、華南理工大學(xué)方利國(guó)編寫(xiě)，年3月12日% 歡送讀者調(diào)用，如有問(wèn)題請(qǐng)告知 clear all;Clcy0 = 1.6 1.2;x1,y1 = ode45(f,0:0.1:10,y0); %0到3分鐘，每分鐘一個(gè)計(jì)算點(diǎn)u=y1(:,1);v=y1(:,2);plot(x1,u,r-)xlabel(時(shí)間，M)ylabel(微生物濃度)hold onplot(x1,v,k:)disp(Results by using ode45()：)disp( x u v )disp(x1 y1)% -function dy = f(x,y) % 定義降溫速率的微分方程f1=0.09*y(1)*(1-y(1)/20)

7、-0.45*y(1)*y(2);f2=0.06*y(2)*(1-y(2)/15)-0.001*y(1)*y(2);dy=f1;f2;1.2 初值問(wèn)題求解 例3：當(dāng)X較大時(shí)，兩種方法計(jì)算結(jié)果有較大不同，為什么？單個(gè)微分方程有零點(diǎn)問(wèn)題？function L43ODEs% 在7.0 版本上調(diào)試通過(guò)% 由華南理工大學(xué)方利國(guó)編寫(xiě)，年2月29日% 歡送讀者調(diào)用，如有問(wèn)題請(qǐng)告知 clear allclcy0 = 1;x1,y1 = ode45(f,0:0.05:10,y0); %0到10，每間隔一個(gè)計(jì)算點(diǎn)x2,y2 = ode23(f,0:0.05:10,y0);%0到10，每間隔一個(gè)計(jì)算點(diǎn)plot(x1,

8、y1,r-)xlabel(x)ylabel(y)hold ondisp(Results by using ode45()：)disp( x y(1) )disp(x1 y1)disp(Results by using ode23()：)disp( x y(2) )disp(x2 y2)plot(x2,y2,b-)% -function dy = f(x,y) % 定義微分方程%dy=y2*cos(x); dy=x2+y*cos(x)計(jì)算值x y(1) Columns 1 through 13 Columns 14 through 26 Columns 27 through 29 高階微分方程求

9、解求解思路：將高階微分方程通過(guò)變量轉(zhuǎn)換，轉(zhuǎn)變成一級(jí)微分方程組進(jìn)展求解。例4：高階微分方程求解程序及解Columns 1 through 13 Columns 14 through 21 剛性方程求解有些微分組的系數(shù)變化很大，這時(shí)用ODE45就很難收斂求解，這時(shí)可用專門(mén)解決此類微分方程的ODE23S來(lái)求解，需要注意的是在解的圖像繪制時(shí)，也需要考慮數(shù)值的波動(dòng)幅度很大，需要引入對(duì)數(shù)坐標(biāo)。例5： dy1/dx=-0.03*y1+1e4*y2*y3 dy2/dx=0.03*y1-2e4*y2*y3-5e7*y22 dy3/dx=5e7*y22 y1(0)=1, y2(0)=0, y3(0)=0funct

10、ion gangxinDEs% 剛性問(wèn)題計(jì)算,在7.0 版本上調(diào)試通過(guò)% 由華南理工大學(xué)方利國(guó)編寫(xiě)，年3月13日% 歡送讀者調(diào)用，如有問(wèn)題請(qǐng)告知 % dy1=-0.03*y1+1e4*y2*y3% dy2=0.03*y1-2e4*y2*y3-5e7*y22% dy3=5e7*y22clear allclcfigurexspan=0 6*logspace(-6,6);y0 = 1 0 0;x1,y1 = ode15s(f,xspan,y0); %用 ode45計(jì)算剛性方程，可能有問(wèn)題u=y1(:,1);v=1e4*y1(:,2);w=y1(:,3);semilogx(x1,u,r-,linewi

11、dth,2)xlabel(x)ylabel(1e4*v)hold onsemilogx(x1,v,k:,linewidth,2)hold onsemilogx(x1,w,g-,linewidth,2)gridaxis(10(-10) 1010 -0.2 1.2)legend(u,v,w)disp(Results by using ode45()：)disp( x u v w )format longdisp(x1 y1)% -function dy = f(x,y) % 定義微分方程f1=-0.03*y(1)+1e4*y(2)*y(3);f2=0.03*y(1)-2e4*y(2)*y(3)-5

12、e7*y(2)2;f3=5e7*y(2)2;dy=f1;f2;f3;1.3 邊值問(wèn)題求解邊值問(wèn)題相對(duì)于初值問(wèn)題而言，多了一個(gè)端點(diǎn)的約束，假設(shè)在高階或微分方程組中端點(diǎn)約束過(guò)多，微分方程組可能無(wú)解，端點(diǎn)約束有一定限制。可以通過(guò)建立離散的方程組，再利用ODE45進(jìn)展求解，但可以利用MATLAB的專用工具求解最好。下面介紹ODE-BVPs的求解器，主要有bvpinit,bvp4c,deval,solinit等。通過(guò)實(shí)際例子介紹這些內(nèi)部函數(shù)的功能。1.3 邊值問(wèn)題求解solinit=bvpinitx,yinit): 產(chǎn)生在初始網(wǎng)格上的初始解，以便bvp4c調(diào)用，其中x為自變量網(wǎng)格，yinit為對(duì)應(yīng)函數(shù)的

13、初值。sol=bvp4c(odefun,BCfun, solinit,)Bcfun:為定義邊界條件方程， Bcfun(ya,yb),其中ya、yb分別表示左右邊界。其他符號(hào)意義同上。deval(sol,xint):計(jì)算任意點(diǎn)處的函數(shù)值。1.3 邊值問(wèn)題求解 例6：程序function BVP4c1%求解兩點(diǎn)邊值問(wèn)題的例如在7.0 版本上調(diào)試通過(guò)% 由華南理工大學(xué)方利國(guó)編寫(xiě)，年3月13日% 歡送讀者調(diào)用，如有問(wèn)題請(qǐng)告知 clear allclca=0;b=10;solinit = bvpinit(linspace(a,b,101),0 0);sol = bvp4c(ODEfun,BCfun,so

14、linit);format shortx = 0:0.1:10;y = deval(sol,x);y1=y(1,:)y2=y(2,:)plot(x,y1,r-)xlabel(x)ylabel(y)hold ongridplot(x,y2,k:)legend(y,dy)disp(Results by using bvp4c：)disp( x y dy )disp(x1 y1)% -function dy = ODEfun(x,y)f1 =y(2);f2=0.05*(1+x2)*y(1)+2;dy=f1;f2;% -function bc = BCfun(ya,yb)bc = ya(1)-40;y

15、b(1)-80;Results by using bvp4c： x y dy Columns 1 through 13 計(jì)算結(jié)果1.3 邊值問(wèn)題求解 例7：主要程序function dy = ODEfun(x,y)f1 =y(2);f2=4*y(1);dy=f1;f2;% -function bc = BCfun(ya,yb)bc = ya(1)-1;yb(1)-exp(4);function BVP4c2%求解兩點(diǎn)邊值問(wèn)題的例如在7.0 版本上調(diào)試通過(guò)% 由華南理工大學(xué)方利國(guó)編寫(xiě)，年3月13日% 歡送讀者調(diào)用，如有問(wèn)題請(qǐng)告知 clear allclca=0;b=2solinit = bvpi

16、nit(linspace(a,b,21),0 0);sol = bvp4c(ODEfun,BCfun,solinit);format shortx = 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)hold ongridplot(x,y2,k:,linewidth,2)hold onplot(x,y3,b:,linewidth,2)legend(y,dy,real-y)disp(Results by using bvp4c：)disp( x

17、y dy )disp(x; y1;y2)% -function dy = ODEfun(x,y)f1 =y(2);f2=4*y(1);dy=f1;f2;% -function bc = BCfun(ya,yb)bc = ya(1)-1;yb(1)-exp(4);2、偏微分方程組求解 2.1 問(wèn)題描繪 2 .2一維動(dòng)態(tài)PDE方程求解2 .3 二維穩(wěn)態(tài)PDE方程求解(自學(xué)問(wèn)題描繪當(dāng)A，B，C為常數(shù)時(shí)，稱為擬線性偏微分方程，當(dāng)A，B，C滿足不同條件時(shí)，分為三種不同的類型：不同類型的方程，MATLAB求解時(shí)，采用不同或一樣的內(nèi)置函數(shù)或工具箱PDE 求解數(shù)學(xué)模型通式 m=0，表示平板，m=1表示圓柱，m

18、=2表示球形，f項(xiàng)表示通量項(xiàng)，,s項(xiàng)表示源項(xiàng)，c項(xiàng)為對(duì)角陣，元素必須大于等于0才可以求解。調(diào)用方法sol=pdepe(m,pdefun, ,iCfun BCfun, xspan,tspan,)Bcfun:為定義邊界條件方程iCfun:為定義初始條件方程。其他符號(hào)意義同上，注意邊界條件必須寫(xiě)成以上形式：應(yīng)用策略Pdepe 內(nèi)部函數(shù)詳細(xì)應(yīng)用時(shí)，需將實(shí)際的偏微分方程對(duì)照標(biāo)準(zhǔn)模型，確定c、f、s函數(shù)的詳細(xì)形式及邊界條件p、q的詳細(xì)形式及m 值。蒸汽入口流體入口，u，t0冷凝液出口流體出口，u，tL套管動(dòng)態(tài)傳熱問(wèn)題原方程：代入詳細(xì)數(shù)據(jù)，并對(duì)長(zhǎng)度歸一化無(wú)量綱處理后，得到以下方程：對(duì)照標(biāo)準(zhǔn)模型，T就是標(biāo)準(zhǔn)模

19、型中的函數(shù)u，m=0, a=0,b=1,t0=0, t f = 1標(biāo)準(zhǔn)模型：初始條件：零時(shí)刻所有位置溫度為30，即u0=30 邊界條件：在零位置處任意時(shí)間溫度為30，在1位置處，偏導(dǎo)為0 pa=u-30,qa=0;pb=0,qb=1 程序清單function l51clcclear allglobal ua ubalpha = 1.0; % cm/sua = 30;ub = 30;m = 0;a = 0;b = 1;t0 = 0;tf =1x = linspace(a,b,11);t = linspace(t0,tf,101);sol = pdepe(m,PDEfun,ICfun,BCfun,x,t);u = sol(:,:,1)% surface plot of the solutionfigure;surf(x,t,u); title