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

1、化工常微分方程和偏微分方程Matlab求解Matlab 求解化工常微分方程和偏微分方程1、常微分方程组求解1.1 问题描绘及Matlab调用命令1.2 初值问题求解1.3 边值问题求解1.4 加权问题求解自学内容2、偏微分方程组求解2.1 问题描绘及一维动态PDE方程求解2 .2 二维求解1、常微分方程组求解1.1 问题描绘及Matlab调用命令 常微分方程：初值问题 常微分方程：两点边值问题 1、常微分方程组求解1.1 问题描绘及Matlab调用命令 微分方程组： 1、常微分方程组求解1.1 问题描绘及Matlab调用命令 高级微分方程： 1、常微分方程组求解1.1 问题描绘及Matlab调

2、用命令 Matlab 调用命令：ODE45:4-5 阶龙格库塔法非刚性ODE23:2-3 阶龙格库塔法非刚性ODE113:可变D-B-M法非刚性ODE15S:基于数值差分的可变阶方法法刚性ODE23S、 ODE23t、 ODE23tb 刚性 1、常微分方程组求解通用调用格式： x,y=ode*(odefun,xspan,y0,)X:自变量向量，在实际调用时取名不一定要用x，也可以用其他名称，只要前后一致即可。Y:应变量向量，在实际调用时取名不一定要用Y，也可以用其他名称，只要前后一致即可。 *:根据不同的问题调用不同格式，如45，23sodefun: 自定义函数的函数名，该函数为Xspan：自

3、变量的积分限，xa,xb,也可以是离散点，x0,x1,x2,xfy0 ：应变量向量的初值: 可以没有该选项，如有，详细应用见下面的实际例子 1.2 初值问题求解例1：某高温物体其温降过程符合以下规律，其中温度T的单位为K，时间 的单位为分钟，零时刻高温物体的温度为2000K，以1分钟作为时间步长，请计算零时刻以后每隔1分钟至170分钟的温度。单个微分方程function xODEs% 铁球从2000K降温曲线,在7.0 版本上调试通过% 由华南理工大学方利国编写，年2月29日% 欢送读者调用，如有问题请告知 clear all；clcy0 = 2000;x1,y1 = ode45(f,0:1:

4、170,y0); %0到170分钟，每分钟一个计算点x2,y2 = ode23(f,0:1:170,y0);plot(x1,y1,r-)xlabel(时间，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 初值问题求解该问题相当与一个自变量，两个应变量问题，初值及微分表达式，可以利用ODE45求解。微分方程组求解程序代码function uvDEs% 微生物消亡问题计算,在7.0 版本上调试通过% 由

6、华南理工大学方利国编写，年3月12日% 欢送读者调用，如有问题请告知 clear all;Clcy0 = 1.6 1.2;x1,y1 = ode45(f,0:0.1:10,y0); %0到3分钟，每分钟一个计算点u=y1(:,1);v=y1(:,2);plot(x1,u,r-)xlabel(时间，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 初值问题求解 例3：当X较大时，两种方法计算结果有较大不同，为什么？单个微分方程有零点问题？function L43ODEs% 在7.0 版本上调试通过% 由华南理工大学方利国编写，年2月29日% 欢送读者调用，如有问题请告知 clear allclcy0 = 1;x1,y1 = ode45(f,0:0.05:10,y0); %0到10，每间隔一个计算点x2,y2 = ode23(f,0:0.05:10,y0);%0到10，每间隔一个计算点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)计算值x y(1) Columns 1 through 13 Columns 14 through 26 Columns 27 through 29 高阶微分方程求

9、解求解思路：将高阶微分方程通过变量转换，转变成一级微分方程组进展求解。例4：高阶微分方程求解程序及解Columns 1 through 13 Columns 14 through 21 刚性方程求解有些微分组的系数变化很大，这时用ODE45就很难收敛求解，这时可用专门解决此类微分方程的ODE23S来求解，需要注意的是在解的图像绘制时，也需要考虑数值的波动幅度很大，需要引入对数坐标。例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% 刚性问题计算,在7.0 版本上调试通过% 由华南理工大学方利国编写，年3月13日% 欢送读者调用，如有问题请告知 % 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计算刚性方程，可能有问题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 边值问题求解边值问题相对于初值问题而言，多了一个端点的约束，假设在高阶或微分方程组中端点约束过多，微分方程组可能无解，端点约束有一定限制。可以通过建立离散的方程组，再利用ODE45进展求解，但可以利用MATLAB的专用工具求解最好。下面介绍ODE-BVPs的求解器，主要有bvpinit,bvp4c,deval,solinit等。通过实际例子介绍这些内部函数的功能。1.3 边值问题求解solinit=bvpinitx,yinit): 产生在初始网格上的初始解，以便bvp4c调用，其中x为自变量网格，yinit为对应函数的

13、初值。sol=bvp4c(odefun,BCfun, solinit,)Bcfun:为定义边界条件方程， Bcfun(ya,yb),其中ya、yb分别表示左右边界。其他符号意义同上。deval(sol,xint):计算任意点处的函数值。1.3 边值问题求解 例6：程序function BVP4c1%求解两点边值问题的例如在7.0 版本上调试通过% 由华南理工大学方利国编写，年3月13日% 欢送读者调用，如有问题请告知 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 计算结果1.3 边值问题求解 例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%求解两点边值问题的例如在7.0 版本上调试通过% 由华南理工大学方利国编写，年3月13日% 欢送读者调用，如有问题请告知 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 问题描绘 2 .2一维动态PDE方程求解2 .3 二维稳态PDE方程求解(自学问题描绘当A，B，C为常数时，称为拟线性偏微分方程，当A，B，C满足不同条件时，分为三种不同的类型：不同类型的方程，MATLAB求解时，采用不同或一样的内置函数或工具箱PDE 求解数学模型通式 m=0，表示平板，m=1表示圆柱，m

18、=2表示球形，f项表示通量项，,s项表示源项，c项为对角阵，元素必须大于等于0才可以求解。调用方法sol=pdepe(m,pdefun, ,iCfun BCfun, xspan,tspan,)Bcfun:为定义边界条件方程iCfun:为定义初始条件方程。其他符号意义同上，注意边界条件必须写成以上形式：应用策略Pdepe 内部函数详细应用时，需将实际的偏微分方程对照标准模型，确定c、f、s函数的详细形式及边界条件p、q的详细形式及m 值。蒸汽入口流体入口，u，t0冷凝液出口流体出口，u，tL套管动态传热问题原方程：代入详细数据，并对长度归一化无量纲处理后，得到以下方程：对照标准模型，T就是标准模

19、型中的函数u，m=0, a=0,b=1,t0=0, t f = 1标准模型：初始条件：零时刻所有位置温度为30，即u0=30 边界条件：在零位置处任意时间温度为30，在1位置处，偏导为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