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

Matlab求解化工常微分方程和偏微分方程方利国Matlab求解化工常微分方程和偏微分方程方利国1Matlab求解化工常微分方程和偏微分方程1、常微分方程（组）求解1.1问题描述及Matlab调用命令1.2初值问题求解1.3边值问题求解1.4加权问题求解（自学内容）2、偏微分方程（组）求解2.1问题描述及一维动态PDE方程求解2

.2二维求解Matlab求解化工常微分方程和偏微分方程1、常微分方程（21、常微分方程（组）求解1.1问题描述及Matlab调用命令

常微分方程：（初值问题）常微分方程：（两点边值问题）

1、常微分方程（组）求解1.1问题描述及Matlab调用命31、常微分方程（组）求解1.1问题描述及Matlab调用命令

微分方程组：

1、常微分方程（组）求解1.1问题描述及Matlab调用命41、常微分方程（组）求解1.1问题描述及Matlab调用命令

高级微分方程：

1、常微分方程（组）求解1.1问题描述及Matlab调用命51、常微分方程（组）求解1.1问题描述及Matlab调用命令Matlab调用命令：ODE45:4-5阶龙格库塔法（非刚性）ODE23:2-3阶龙格库塔法（非刚性）ODE113:可变D-B-M法（非刚性）ODE15S:基于数值差分的可变阶方法法（刚性）ODE23S、ODE23t、ODE23tb（刚性）

1、常微分方程（组）求解1.1问题描述及Matlab调用命61、常微分方程（组）求解通用调用格式：[x,y]=ode***(@odefun,xspan,y0,<option,p1,p2..>)X:自变量向量，在实际调用时取名不一定要用x，也可以用其他名称，只要前后一致即可。Y:应变量向量，在实际调用时取名不一定要用Y，也可以用其他名称，只要前后一致即可。***:根据不同的问题调用不同格式，如45，23s@odefun:自定义函数的函数名，该函数为Xspan：自变量的积分限，[xa,xb],也可以是离散点，[x0,x1,x2,…xf]y0：应变量向量的初值<>:可以没有该选项，如有，具体应用见下面的实际例子

1、常微分方程（组）求解通用调用格式：71.2初值问题求解例1：已知某高温物体其温降过程符合以下规律，其中温度T的单位为K，时间

的单位为分钟，零时刻高温物体的温度为2000K，以1分钟作为时间步长，请计算零时刻以后每隔1分钟至170分钟的温度。单个微分方程1.2初值问题求解例1：已知某高温物体其温降过程符合以下8functionxODEs%铁球从2000K降温曲线,在7.0版本上调试通过%由华南理工大学方利国编写，2012年2月29日%欢迎读者调用，如有问题请告知lgfang@clearall；clcy0=[2000];[x1,y1]=ode45(@f,[0:1:170],y0);%0到170分钟，每分钟一个计算点[x2,y2]=ode23(@f,[0:1:170],y0);plot(x1,y1,'r-')xlabel('时间，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初值问题求解该问题相当与一个自变量，两个应变量问题，已知初值及微分表达式，可以利用ODE45求解。微分方程组求解1.2初值问题求解该问题相当与一个自变量，两个应变量问题12程序代码functionuvDEs%微生物消亡问题计算,在7.0版本上调试通过%由华南理工大学方利国编写，2012年3月12日%欢迎读者调用，如有问题请告知lgfang@clearall;Clcy0=[1.61.2];[x1,y1]=ode45(@f,[0:0.1:10],y0);%0到3分钟，每0.1分钟一个计算点u=y1(:,1);v=y1(:,2);plot(x1,u,'r-')xlabel('时间，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初值问题求解

例3：当X较大时，两种方法计算结果有较大不同，为什么？单个微分方程有零点问题？1.2初值问题求解例3：当X较大时，两种方法计算结果有14functionL43ODEs%在7.0版本上调试通过%由华南理工大学方利国编写，2012年2月29日%欢迎读者调用，如有问题请告知lgfang@clearallclcy0=[1];[x1,y1]=ode45(@f,[0:0.05:10],y0);%0到10，每0.05间隔一个计算点[x2,y2]=ode23(@f,[0:0.05:10],y0);%0到10，每0.05间隔一个计算点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计算值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计算值xy(1)16高阶微分方程求解求解思路：将高阶微分方程通过变量转换，转变成一级微分方程组进行求解。例4：高阶微分方程求解求解思路：将高阶微分方程通过变量转换，转变成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刚性方程求解有些微分组的系数变化很大，这时用ODE45就很难收敛求解，这时可用专门解决此类微分方程的ODE23S来求解，需要注意的是在解的图像绘制时，也需要考虑数值的波动幅度很大，需要引入对数坐标。例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刚性方程求解有些微分组的系数变化很大，这时用ODE45就很难19functiongangxinDEs%刚性问题计算,在7.0版本上调试通过%由华南理工大学方利国编写，2012年3月13日%欢迎读者调用，如有问题请告知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计算刚性方程，可能有问题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边值问题求解边值问题相对于初值问题而言，多了一个端点的约束，如果在高阶或微分方程组中端点约束过多，微分方程组可能无解，端点约束有一定限制。可以通过建立离散的方程组，再利用ODE45进行求解，但可以利用MATLAB的专用工具求解最好。下面介绍ODE-BVPs的求解器，主要有bvpinit,bvp4c,deval,solinit等。通过实际例子介绍这些内部函数的功能。1.3边值问题求解边值问题相对于初值问题而言，多了一个端点211.3边值问题求解solinit=bvpinit（x,yinit):产生在初始网格上的初始解，以便bvp4c调用，其中x为自变量网格，yinit为对应函数的初值。sol=bvp4c(@odefun,@BCfun,solinit,<option,p1,p2..>)@Bcfun:为定义边界条件方程，Bcfun(ya,yb),其中ya、yb分别表示左右边界。其他符号意义同上。deval(sol,xint):计算任意点处的函数值。1.3边值问题求解solinit=bvpinit（x,yi221.3边值问题求解

例6：1.3边值问题求解例6：23程序functionBVP4c1%求解两点边值问题的示例在7.0版本上调试通过%由华南理工大学方利国编写，2012年3月13日%欢迎读者调用，如有问题请告知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计算结果计算结果251.3边值问题求解

例7：1.3边值问题求解例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%求解两点边值问题的示例在7.0版本上调试通过%由华南理工大学方利国编写，2012年3月13日%欢迎读者调用，如有问题请告知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问题描述2

.2一维动态PDE方程求解

2

.3二维稳态PDE方程求解(自学）

2、偏微分方程（组）求解2.1问题描述29问题描述当A，B，C为常数时，称为拟线性偏微分方程，当A，B，C满足不同条件时，分为三种不同的类型：不同类型的方程，MATLAB求解时，采用不同或相同的内置函数或工具箱问题描述当A，B，C为常数时，称为拟线性偏微分方程，当A，B30PDE求解数学模型通式m=0，表示平板，m=1表示圆柱，m=2表示球形，f项表示通量项，,s项表示源项，c项为对角阵，元素必须大于等于0才可以求解。PDE求解数学模型通式m=0，表示平板，m=1表示圆柱，31调用方法sol=pdepe(m,@pdefun,,@iCfun@BCfun,xspan,tspan,<option,p1,p2..>)@Bcfun:为定义边界条件方程@iCfun:为定义初始条件方程。其他符号意义同上，注意边界条件必须写成以上形式：调用方法sol=pdepe(m,@pdefun,,@iCf32应用策略Pdepe内部函数具体应用时，需将实际的偏微分方程对照标准模型，确定c、f、s函数的具体形式及边界条件p、q的具体形式及m值。蒸汽入口流体入口，u，t0冷凝液出口流体出口，u，tL应用策略Pdepe内部函数具体应用时，需将实际的偏微分方程33套管动态传热问题原方程：代入具体数据，并对长度归一化无量纲处理后，得到以下方程：套管动态传热问题原方程：代入具体数据，并对长度归一化无量纲处34对照标准模型，T就是标准模型中的函数u，m=0,

a=0,b=1,t0=0,tf=1标准模型：初始条件：零时刻所有位置温度为30，即u0=30

边界条件：已知在零位置处任意时间温度为30，在1位置处，偏导为0pa=u-30,qa=0;pb=0,qb=1

对照标准模型，T就是标准模型中的函数u，m=0,a=0,b35程序清单functionl51clcclearallglobaluaubalpha=1.0;%cm/sua=30;ub=30;m=0;a=0;