化工常微分方程和偏微分方程Matlab求解_精品文档.ppt

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、常微分方程（组）求

2、解1.1 问题描述及Matlab调用命令 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

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

4、0=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)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*

5、exp(0.001*(y(1)-300)*(-300+y(1);Results by using ode45()：x y(1)1.0e+003*Columns 1 through 13 0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 2.0000 0.8788 0.6133 0.4880 0.4181 0.3762 0.3498 0.3328 0.3218 0.3145 0.3097 0.3065 0.3043 Columns 14 through 18 0.1300 0

6、.1400 0.1500 0.1600 0.1700 0.3029 0.3019 0.3013 0.3009 0.3006Results by using ode23()：x y(2)1.0e+003*Columns 1 through 13 0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000 0.1100 0.1200 2.0000 0.8779 0.6124 0.4873 0.4176 0.3756 0.3493 0.3323 0.3213 0.3140 0.3092 0.3061 0.3041 C

7、olumns 14 through 18 0.1300 0.1400 0.1500 0.1600 0.1700 0.3027 0.3018 0.3012 0.3008 0.30051.2 初值问题求解该问题相当与一个自变量，两个应变量问题，已知初值及微分表达式，可以利用ODE45求解。微分方程组求解微分方程组求解程序代码function uvDEs%微生物消亡问题计算,在7.0 版本上调试通过%由华南理工大学方利国编写，2012年3月12日%欢迎读者调用，如有问题请告知 clear all;Clcy0=1.6 1.2;x1,y1=ode45(f,0:0.1:10,y0);%0到3分钟，每0.1

8、分钟一个计算点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)-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较大时，两种方法计算结果有较大不同，为什么？单

9、个微分方程有单个微分方程有零点问题？零点问题？function L43ODEs%在7.0 版本上调试通过%由华南理工大学方利国编写，2012年2月29日%欢迎读者调用，如有问题请告知 clear allclcy0=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)hold ondisp(Results by using ode45()：)disp(x y(1)disp(x1 y1)disp(Re

10、sults 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 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 1.0000 1.0526 1.1109 1.1757 1.2479 1.3287 1.4195 1.5218 1.6378 1.7698

11、 1.9210 2.0951 2.2970 Columns 14 through 26 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500 1.2000 1.2500 2.5329 2.8107 3.1411 3.5381 4.0206 4.6152 5.3598 6.3078 7.5435 9.1928 11.4614 14.7283 19.5991 Columns 27 through 29 1.3000 1.3500 1.4000 27.4283 41.2030 68.6630高阶微分方

12、程求解求解思路：将高阶微分方程通过变量转换，转变成一级微分方程组进行求解。例4：高阶微分方程求解程序及解Columns 1 through 13 0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000 0.9000 1.0000 1.1000 1.2000 0 0.1002 0.2015 0.3052 0.4129 0.5263 0.6476 0.7789 0.9230 1.0827 1.2614 1.4627 1.6908 1.0000 1.0053 1.0227 1.0544 1.1026 1.1698 1.2588 1.3723

13、 1.5136 1.6861 1.8935 2.1398 2.4296 Columns 14 through 21 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000 1.9502 2.2461 2.5841 2.9705 3.4123 3.9171 4.4936 5.1513 2.7677 3.1595 3.6110 4.1286 4.7196 5.3918 6.1541 7.0159刚性方程求解有些微分组的系数变化很大，这时用ODE45就很难收敛求解，这时可用专门解决此类微分方程的ODE23S来求解，需要注意的是在解的图像绘制时，

14、也需要考虑数值的波动幅度很大，需要引入对数坐标。例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)=0function gangxinDEs%刚性问题计算,在7.0 版本上调试通过%由华南理工大学方利国编写，2012年3月13日%欢迎读者调用，如有问题请告知 %dy1=-0.03*y1+1e4*y2*y3%dy2=0.03*y1-2e4*y2*y3-5e7*y22%dy3=5e7*y22clear allclcfigurexspan=0 6*logs

15、pace(-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-,linewidth,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(

16、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)-5e7*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=bvpinit（x,yinit):产生在初始网格上的初始解，以便bvp4c调用，其中x为自变量网格，yinit为对应函数的初值。sol=bvp4c(odefun,BCfun,solinit,)Bcfun:为定义边界条件方程，Bcfun(ya,yb),其中ya、yb分别表示左右边界。其他符号意义同