第4章常微分方程数值解

第4章

常微分方程数值解4.1微分方程在化工中的应用

4.2欧拉（Euler）公式

4.3龙格-库塔方法

4.4常微分方程组的数值解法

4.5程序示例及应用

总目录4.1微分方程在化工中的应用微分方程在化工中应用的简单而又典型的例子是套管式换热器的稳态温度分布。首先作以下假设：1、套管内侧为液体，其温度只随套管的长度改变而改变，忽略温度的径向变化；套管环隙为蒸汽，其温度在任何位置均为恒定值，可认为是饱和蒸汽的温度。2、忽略套管内侧流体的纵向热传导。3、在整个套管长度方向上，总传热系数K不变。总目录本章目录4.14.24.34.44.54.1微分方程在化工中的应用蒸汽入口流体入口，u，t0冷凝液出口流体出口，u，tL图4-1套管式换热器温度分布示意图流入的热量+传入的热量-流出的热量=0（4-1）（4-2）（4-3）总目录本章目录4.14.24.34.44.54.1微分方程在化工中的应用另一个在化工中常见的微分方程是物料冷却过程的数学模型，其模型可用下式表示：

（4-3）在微分方程中我们称自变量函数只有一个的微分方程为常微分方程，自变量函数个数为两个或两个以上的微分方程为偏微分方程。给定微分方程及其初始条件，称为初值问题；给定微分方程及其边界条件，称为边值问题。总目录本章目录4.14.24.34.44.54.1微分方程在化工中的应用在化工模拟中主要碰到的是常微分方程的初值问题：（4-5）或对于大多数常微分方程的初值问题,只能计算它的数值解。常微分方程初值问题的数值解就是求y(x)在求解区间[a,b]上各个分点序列xn，n=1,2,…,m的数值解yn。在计算中约定y(xn)表示常微分方程准确解的值，yn表示y(xn)的近似值。总目录本章目录4.14.24.34.44.54.2欧拉（Euler）公式4.2.1向前欧拉公式4.2.2向后欧拉公式

4.2.3中心欧拉公式4.2.4梯形公式

总目录本章目录4.14.24.34.44.54.2.1向前欧拉公式（4-6）下式为计算近似值的向前欧拉公式：图4-2欧拉折线法几何示意图总目录本章目录4.14.24.34.44.54.2.1向前欧拉公式实例例4.1：假定某物体的温度w因自热而产生的热量可以使物体在每秒钟内以4%的速度增长，同时该物体由于散热可使其温度在每秒种内下降100k，则物体温度随时间变化的微分方程：（t以秒为单位)

分别以初始温x(0)=1500k，y(0)=2500k，z(0)=3500k用欧拉公式预测24秒后的物体温度趋势。总目录本章目录4.14.24.34.44.54.2.1向前欧拉公式实例解：

w0分别以x0=1500，y0=2500，z0=3500代入。计算结果见表4-1。图4-3三种初始值的温度变化曲线

表4-1总目录本章目录4.14.24.34.44.54.2.1向前欧拉公式实例从表4-1可以看到当自热引起物体温度升高的速度小于散热引起温度下降的速度，物体的温度随时间而逐渐减少：当自热引起物体温度升高的速度与散热引起温度下降的速度平衡时，物体的温度保持不变；当自热引起物体温度升高的速度大于散热引起温度下降的速度，物体的温度随时间而增长。在图4-3中L1，L2，L3分别表示初始值3500，2500和1500的三条温度变化趋势曲线。总目录本章目录4.14.24.34.44.5h充分小时，以上迭代收敛。记，则

h充分小时，可保证，其中L为李普希兹条件。4.2.2向后欧拉公式向后欧拉公式：（4-6）式(4-7)是yn+1的非线性方程，即隐式欧拉公式，用迭代法求得yn+1。初始值由向前欧拉公式提供。最简单的迭代公式为：总目录本章目录4.14.24.34.44.54.2.3中心欧拉公式y(x)的在x=x1处的中心差商式：又，可得到y(x2)的近似值y2计算公式：类似地，可得到计算y(xn+1)近似值yn+1的计算公式：

公式(4-8)称为中心格式。按公式(4-8)，需要知道yn-1,yn的值才能求得yn+1的值。因此，要先用其它公式计算出y1，再用中心格式算出y2,y3,…。y1可用向前欧拉公式计算，为提高精度，也可用向后欧拉公式计算。(4-8)总目录本章目录4.14.24.34.44.54.2.4梯形公式梯形公式：(4-9)梯形公式也是隐式格式。用显式的欧拉公式和隐式的梯形公式给出的一次预估-校正公式：(4-10)上式也称为改进的欧拉公式，它可合并成：总目录本章目录4.14.24.34.44.54.2.4梯形公式实例例4.2：请用预估-校正公式（改进的欧拉公式）解右面初值问题：解：用下面的迭代公式，对每个点迭代4次，k=1,2,3,4。

总目录本章目录4.14.24.34.44.54.2.4梯形公式实例该方程的精确解是

计算结果如表4-2所示。表4-2计算结果总目录本章目录4.14.24.34.44.54.3龙格-库塔方法本课程主要研究其实际应用,，故直接给出各类龙格-库塔公式。1、二阶龙格-库塔其中c1=0,c2=1,a=1/2,b=1/2。其中c1=1/2,c2=1/2,a=1,b=1或(4-11)(4-12)总目录本章目录4.14.24.34.44.54.3龙格-库塔方法2、三阶龙格-库塔公式

(4-15)(4-14)(4-13)总目录本章目录4.14.24.34.44.54.3龙格-库塔方法3、四阶龙格—库塔公式

(4-17)(4-16)总目录本章目录4.14.24.34.44.54.3龙格-库塔方法实例例4.3：用四阶龙格—库塔公式(4-16)求解下面初值问题解：取步长h=0.2，计算公式为：表4-3计算结果总目录本章目录4.14.24.34.44.54.3龙格-库塔方法步长的选择下面以四阶龙格-库塔方法为例，说明如何自动选择步长，使计算结果满足给定精度的要求。设从节点xn出发，先以h为步长，利用四阶龙格-库塔公式方法经过一步计算得y(xn+1)的近似值，记为，由于公式的局部截断误差是y(h5)，故有当h不大时，c可近似地看作常数。然后将步长h对折，即取h/2为步长，从出发经过两步计算求y(xn+1)的近似值，记为，每一步计算的局部截断误差为c(h/2)5，于是就有(4-18)总目录本章目录4.14.24.34.44.54.3龙格-库塔方法步长的选择把它与（4-18）式相比，可得：经整理可得：总目录本章目录4.14.24.34.44.54.3龙格-库塔方法步长的选择这表明以作为y(xn+1)的近似值，其误差可用先后两次计算结果之差来表示，因而，只需考察是否成立。若成立，则可将作为y(xn+1)的近似值；若不成立，则将步长再次对折进行计算，直到不等式成立为止，并取最后的作为计算结果。以上方法就是计算过程中自动选择步长的方法，也称为变步长方法。总目录本章目录4.14.24.34.44.54.3龙格-库塔方法步长的选择用一个例子说明：由VB选用标准四阶龙格-库塔方法计算得

因此合理选择步长既能保证精度又能减少计算量。总目录本章目录4.14.24.34.44.54.4常微分方程组的数值解法4.4.1一阶常微分方程组的数值解法4.4.2高阶常微分方程数值方法

总目录本章目录4.14.24.34.44.54.4.1一阶常微分方程组的数值解法将由m个一阶方程组成的常微分方程初值问题：向量形式：其中：(4-19)(4-20)总目录本章目录4.14.24.34.44.54.4.1一阶常微分方程组的数值解法下面以两个方程组为例，给出相应的计算公式。常微分方程组：欧拉公式：预估—校正公式：

(4-21)(4-22)总目录本章目录4.14.24.34.44.54.4.1一阶常微分方程组的数值解法四阶龙格—库塔公式：(4-23)总目录本章目录4.14.24.34.44.54.4.1一阶常微分方程组的数值解法

实例例4.4：两种微生物，其数量分别是u=u(t)，v=v(t)，t的单位为分，其中一种微生物以吃另一种微生为生，两种微生物的增长函数如下列常微分方程组所示，预测3分钟后这一对微生物的数量。总目录本章目录4.14.24.34.44.54.4.1一阶常微分方程组的数值解法

实例解：记用欧拉预估—校正公式(4-22)表4-4计算结果总目录本章目录4.14.24.34.44.54.4.2高阶常微分方程数值方法以三阶常微分方程为例说明高阶常微分方程的数值计算步骤。

(4-24)令得到一阶方程组：总目录本章目录4.14.24.34.44.54.5程序示例及应用任务1用改进的欧拉公式求解常微分方程初值问题：

算法描述对给定的F(x,y)，用改进的欧拉公式求解常微分方程初值问题的解。

(计算实例VB程序见课本)总目录本章目录4.14.24.34.44.54.5程序示例及应用任务2：用四阶龙格—库塔方法求解常微分方程初值问题：

算法描述 对给定的F（x,y），用四阶龙格—库塔方法求解常微分方程初值问题。(计算实例VB程序见