化工计算方法6微分方程

6常微分方程数值解法6.1基本概念及求解思路

微分方程：表示未知函数与未知函数的导数以及自变量之间关系的方程。（含有未知函数的导数的方程)只含一个自变量的导数的方程叫做常微分方程。方程中出现的导数的最高阶数叫微分方程的阶一阶常微分方程二阶常微分方程根据所给的(已知)定解条件，将常微分方程分为初值问题和边值问题两类给出自变量两端的函数值或导数值，称边值条件#初值问题给出自变量为0或通过变换可将自变量变为0（自变量一端）时的函数值或导数值，称初值条件；边值问题例：在一个封闭系统中，三个组分浓度分别为C1、C2、C3。系统在特定频率光照下发生反应，求每种物质浓度随时间的变化规律。已知自变量t为0时的变量值（初始条件），为常微分方程组的初值问题。#

解：此问题的数学方程为已知：C1C2C3k1C1k3C22k2C2C3例：一维均匀介质稳态导热问题，设其一端绝热，另一端恒温为T1，此问题的数学模型为:绝热T1L常微分方程初值问题的一般形式表示为已知自变量两端的条件，属常微分方程的边值问题

常微分方程的初值问题和边值问题求解方法不同,本课程只涉及初值问题的求解。一阶常微分方程初值问题的一般形式解的形式为y=g(x），是在xy平面上通过初始点（x0，y0）的一条曲线，是连续函数。x1x2…xixn…y=g(x）数值解并不是求这个连续函数的形式y=g(x），而是寻找函数在一系列离散点xi

上的函数的近似值yi用有限个离散点上的近似值的集合代替连续解。#微分方程数值求解的基本原理和步骤将求解区间分成若干离散点x1，x2，…，xi，…；在离散点上用差商近似地代替导数，将微分方程化为代数方程（差分方程），差分方程中的未知量就是未知函数在离散点xi上的近似值yi；常微分方程初值问题数值求解基本思路是在离散点上用差商代替导数，基本方法则是递推，从初始值开始，沿着节点排列的顺序向前推进。3.求解代数方程得到未知函数在离散点上的近似值y1，y2，…，yi

，…；x1x2…xixn…将点

x1,x2,…处的导数用差商表示，微分方程变成代数方程6.2欧拉（Euler）法2.在各离散点xi应用差商近似代替一阶导数一阶常微分方程的一般形式步长代入微分方程式，得点xi处的差分方程求

x=x0～xs

间y的变化情况（y随x变化的函数）

将求解区间分成若干离散点：x0，x1，x2，。。。，xi，。。。

xs欧拉（Euler）格式上式中，令欧拉格式：沿着节点排列的顺序向前推进，可逐步求得各离散点上的近似数值解差商代替微分方程式中的导数，得点xi处的差分方程y=g(x)Euler格式只有左端有未知数的格式，称为显式格式

是一条直线，即点(xi,yi)处的切线，故Euler格式的几何意义:用每一步起点的切线代替原方程的积分曲线，用连续折线代替连续曲线。＃从出发递推，求出代数方程组的解未知数xixi+1xn第一步第二步第i步递推求解:隐式格式无法直接计算未知量，故与Euler格式结合求解。方程两端均有未知数，称为隐式格式Euler格式2为改善精度，将Euler格式中f(xi,yi)用xi和xi＋1处的平均值代替，称为改进Euler格式。在Euler格式中，f（xi，yi）是xi

处的函数值，也是点（xi，yi）处切线的斜率，精度较低，误差积累大。xixi+1xn改进Euler格式几何意义:以点（xi,yi）和（xi+1,yi+1）处的斜率的算术平均值代表区间（xi,

xi＋1）上的斜率。Euler格式预报梯形公式校正Euler格式3采用预报－校正的方式求解改进Euler格式计算误差通常用局部截断误差来估计，用步长h的方次表示。误差的大小与方法的精确程度有关，如果一种方法的局部截断误差为，则称这种方法是p阶方法具有p阶精度。分析可知，欧拉法的精度为一阶，而改进欧拉法精度为二阶。步长h＝0.1，求x=0~1间各节点上的数值解。

解

先用Euler格式计算例6-1用改进欧拉法解初值问题Euler格式预报梯形公式校正再用改进欧拉法求解比较Euler格式结果:y1=1.1比较Euler格式结果:y2=1.1918182#No.X

Y

Y*()001110.11.0959091.09544520.21.1840971.18321630.31.2662011.26491140.41.3433601.34164150.51.4164021.41421460.61.4859551.48324070.71.5525141.54919380.81.6164741.61245290.91.6781661.673320101.01.7378671.732051比较数值解与解析解，改进Euler法至少有三位有效数字运行结果本题解析解为，与数值解一并列出以比较计算精度。计算程序参考教材p.66运行结果解析解结果例6-2

已知在管式反应器中进行液相反应A→R+S，反应为吸热反应，反应管外油浴温度为340℃，假定已知管内温度与转化率的关系为Rg=1.987，tc

=反应器外壁温度。如反应器入口温度为340℃，反应器出口转化率为90％，试用改进欧拉格式求反应器出口温度。速度常数为出口：转化率xA＝90％，t=？℃

#入口：转化率xA＝0，t=340℃例6-2求解xA0

0.05

0.10

0.15

0.20

0.250.30t340336.85

333.93

331.34

329.15327.42326.18xA0.350.400.450.500.550.60

0.65t

325.40325.02324.97325.20325.62326.21326.94xA0.70

0.750.800.850.90

t327.67

328.60329.87331.32

330.56

使用改进欧拉法计算，设自变量xA的步长△xA＝0.05，经过18步，可求出xA＝0.9时的出口温度程序分为两部分：(a)定义方程右端函数；(b)用改进欧拉法计算。参考程序见教材p.68。解：由题意知为常微分方程的初值问题，假定tc＝340℃（管外油浴温度），方程形式为6.3龙格-库塔（Runge-Kutta）法

龙格-库塔（Runge-Kutta）法使用更多的斜率来构造计算格式，所用的斜率个数称为阶数；根据所用斜率个数不同，有三阶、四阶、六阶龙格-库塔格式，一般来说，阶数越高，精度越高一个斜率在初值问题的数值解法中最常用的是经典的四阶龙格-库塔格式，其计算公式为Euler格式改进Euler格式如令二个斜率，精度高于用一个斜率四阶龙格-库塔格式用4个点处的斜率值，求其加权平均值作为小区间内的斜率值四阶龙格-库塔格式是4阶精度O(h4)四阶龙格-库塔-吉尔格式龙格-库塔-吉尔格式是四阶龙格-库塔格式的一种改进对有些问题的计算比经典的龙格-库塔格式稳定。MATLAB有龙格-库塔函数可直接调用仍用4个斜率值，但取的权重与龙格-库塔格式不同，递推格式也有差异常微分方程组及及高阶方程求解自变量x为一个，因变量(y1，y2，…，yn）和方程数（f1，f2，…，fn）均大于一个时构成常微分方程组。

2个变量2个方程构成方程组的一般形式为含n个方程的方程组的一般形式（向量形式）为求解时，需将龙格-库塔法解单个方程的计算格式推广到方程组的情形。例：2个方程构成的方程组求常微分方程组的四阶龙格-库塔格式

yji为第j个变量yj(x)在节点xi

处的近似解；

n为因变量和方程的个数高阶微分方程求解对于一阶以上的高阶微分方程，通常可将其化为一阶方程组求解。对以下形式的二阶常微分方程的初值问题一个二阶常微分方程转换为由两个一阶方程构成的常微分方程组。一个n阶的常微分方程可通过变换化为n个一阶方程构成的的常微分方程组，然后可按照一阶常微分方程组求解。初始条件改为例6-3：对例6-1改用四阶龙格-库塔法计算,比较改进欧拉法与四阶龙格-库塔法两种方法计算的精度。计算结果显示为（y*为精确值）

i x y

y* 0 0.00 1.000000

1.000000 1 0.10 1.095445

1.095445 2 0.20 1.1832161.183216 3 0.30 1.2649111.264911 4 0.40 1.3416411.341641 5 0.50 1.4142141.414214 6 0.60 1.4832401.483240 7 0.70 1.5491931.549193 8 0.80 1.6124521.612452 9 0.90 1.6733201.673320 10 1.00 1.7320511.732051解为编程计算方便，直接调用MATLAB函数库中龙格-库塔函数。结果中列出了解析解y*以进行精度比较。比较可知，步长h=0.1时，四阶龙格-库塔法有效数字可达到7位例6-4

求解以下微分方程组，区间为[0,0.2]，步长取h=0.01，要求每隔5步输出一次计算结果。解：本题为3个方程构成的常微分方程组，注意程序中右端函数的定义，直接调用函数计算。参考程序见教材p.71。计算结果显示

x

y1

y2y3 0.00 -1.000000 0.000000 1.000000 0.05 -0.998750 0.049979 0.951229 0.10 -0.995004 0.099833 0.904837 0.15 -0.988771 0.149438 0.860708 0.20 -0.980067 0.198669 0.818731例题6-7甲烷水蒸气催化转化反应制氢，计算：

1、从反应管入口到转化率xA大于0.644时，反应管的轴向温度(T)及转化率(xA)分布；

2、所需要的反应管长度和催化剂用量。反应管内径102mm(10根)假设反应管径向温度和组成分布均匀，反应管内温度和转化率只沿轴向（管长）变化甲烷水蒸气入口温度

t0＝650℃催化剂烟道气加热轴向径向已知物料衡算和热量衡算式如下步长取h＝△L=△x=0.5m初始条件为：解(1)本例为两个方程构成的常微分方程组自变量：管长L，令L=x；

因变量：轴向温度T，转化率xA

令

xA＝y(1)，t=y(2)物料衡算xA0＝0t0＝650℃甲烷水蒸气热量衡算xA(L)=？t(L)=？例题6-7程序(a)本题求温度T和转化率xA的分布，计算转化率xA≥0.644时的反应管长度对所给的微分方程组，从反应管入口L=0开始，到xA≥0.644时为止，每隔0.5m(步长)进行一次计算，即可得到温度和转化率的分布情况，计算到xA≥0.644时的反应管长，即为所要求的管长；(b)计算10根反应管的催化剂用量积分求得反应管长度L后，即可计算催化剂用量

m3调用函数ode计算，程序简短，计算快捷。参考程序见教材p.74。

运行后输出结果L=7.00mVr=0.5720m3L(m)xaT(K)0.000.0000923.150.500.0618916.761.000.1172915.721.500.1695917.082.000.2198919.622.500.2688922.813.000.3166926.403.500.3633930.284.000.4089934.414.500.4534938.775.000.4968943.375.500.5390948.216.000.5798953.326.500.6194958.707.000.6575964.39如选择小步长(如h=0.01m)，计算结果可更精确。进料温度＝620℃水蒸气与乙苯的比例为10比1质量流速G＝2500kg/m2·h反应混合物平均热容Cp＝2.18kJ/kg·K反应热H＝140000kJ/kmol催化剂堆密度＝1440m3/kg绝热反应器例题6-8用一维拟均相模型计算以下两种情况的乙苯脱氢制苯反应器。

（1）采用绝热反应器乙苯水蒸气传热系数

U＝0反应器直径D＝1.211m（2）采用管式换热反应器分别计算当乙苯转化率达到50％时两种反应器的(1)温度分布(2)转化率分布(3)反应器（管）长度产量、反应管