化工常微分方程和偏微分方程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分别表示左右边界。
其他符号意义同
升级会员 