化工常微分方程和偏微分方程Matlab求解_精品文档.ppt
- 文档编号:2562781
- 上传时间:2022-11-01
- 格式:PPT
- 页数:40
- 大小:556.50KB
化工常微分方程和偏微分方程Matlab求解_精品文档.ppt
《化工常微分方程和偏微分方程Matlab求解_精品文档.ppt》由会员分享,可在线阅读,更多相关《化工常微分方程和偏微分方程Matlab求解_精品文档.ppt(40页珍藏版)》请在冰豆网上搜索。
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调用命令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:
自变量的积分限,xa,xb,也可以是离散点,x0,x1,x2,xfy0:
应变量向量的初值:
可以没有该选项,如有,具体应用见下面的实际例子1.2初值问题求解例1:
已知某高温物体其温降过程符合以下规律,其中温度T的单位为K,时间的单位为分钟,零时刻高温物体的温度为2000K,以1分钟作为时间步长,请计算零时刻以后每隔1分钟至170分钟的温度。
单个微分方程单个微分方程functionxODEs%铁球从2000K降温曲线,在7.0版本上调试通过%由华南理工大学方利国编写,2012年2月29日%欢迎读者调用,如有问题请告知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);Resultsbyusingode45():
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.30051.2初值问题求解该问题相当与一个自变量,两个应变量问题,已知初值及微分表达式,可以利用ODE45求解。
微分方程组求解微分方程组求解程序代码functionuvDEs%微生物消亡问题计算,在7.0版本上调试通过%由华南理工大学方利国编写,2012年3月12日%欢迎读者调用,如有问题请告知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;1.2初值问题求解例3:
当X较大时,两种方法计算结果有较大不同,为什么?
单个微分方程有单个微分方程有零点问题?
零点问题?
functionL43ODEs%在7.0版本上调试通过%由华南理工大学方利国编写,2012年2月29日%欢迎读者调用,如有问题请告知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=y2*cos(x);dy=x2+y*cos(x)计算值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高阶微分方程求解求解思路:
将高阶微分方程通过变量转换,转变成一级微分方程组进行求解。
例4:
高阶微分方程求解程序及解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刚性方程求解有些微分组的系数变化很大,这时用ODE45就很难收敛求解,这时可用专门解决此类微分方程的ODE23S来求解,需要注意的是在解的图像绘制时,也需要考虑数值的波动幅度很大,需要引入对数坐标。
例5:
dy1/dx=-0.03*y1+1e4*y2*y3dy2/dx=0.03*y1-2e4*y2*y3-5e7*y22dy3/dx=5e7*y22y1(0)=1,y2(0)=0,y3(0)=0functiongangxinDEs%刚性问题计算,在7.0版本上调试通过%由华南理工大学方利国编写,2012年3月13日%欢迎读者调用,如有问题请告知%dy1=-0.03*y1+1e4*y2*y3%dy2=0.03*y1-2e4*y2*y3-5e7*y22%dy3=5e7*y22clearallclcfigurexspan=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)1010-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;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分别表示左右边界。
其他符号意义同
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 化工 微分方程 Matlab 求解 精品 文档