第4章常微分方程数值解

1、第第4章章 常微分方程数值解常微分方程数值解 4.1微分方程在化工中的应用微分方程在化工中的应用 4. 2 欧拉（欧拉（euler）公式）公式 4.3 龙格龙格-库塔方法库塔方法 4.4 常微分方程组的数值解法常微分方程组的数值解法 4.5 程序示例程序示例 目录目录4.1微分方程在化工中的应用微分方程在化工中的应用 微分方程在化工中应用的简单而又典型的例子是套管式换热器的稳态温度分布。首先作以下假设： 1、套管内侧为液体，其温度只随套管的长度改变而改变，忽略温度的径向变化；套管环隙为蒸汽，其温度在任何位置均为恒定值，可认为是饱和蒸汽的温度。 2、忽略套管内侧流体的纵向热传导。 3、在整个套管

2、长度方向上，总传热系数k不变。 据以上假设，可以得到图中所示微元的能量平衡方程（示意图见图4- 1）： 流入的热量+传入的热量-流出的热量=0 （4-1） 即： 蒸汽入口流体入口，u，t0冷凝液出口流体出口，u，tl0)()(222dldtdltucrttrdlktucrpwp4.14.44.34. 2总目录总目录4.54.1微分方程在化工中的应用微分方程在化工中的应用 化简上式得其温度的微分方程： （4-3） 微分方程中各变量的含义如下： 通过求解微分方程（4-3），就可以得到管内流体的温度随管子长度而改变的曲线，为化工模拟和设计提供依据。如果方程（4-3）中传热系数、物流性质不随温度或位置

3、的变化，那么，方程（4-3）是可以解析求解的（数值求解当然可以），得到我们常见的换热器传热方程；如果上述性质可能随温度或位置而改变，那么方程（4-3）就只能用数值的方法求解了。事实上传热系数也好，物流性质也好都会随着温度的改变而改变，故在深入研究换热器各点温度分布时，拟采用微分方程的数值求解为好。 )( 2ttrcukdldtwpmrmkgkkgjcsmuktlktpw内套管半径，套管内流体的密度套管内流体的比热套管内流体的速度套管的管壁温度置流体在套管内所处的位套管内某一点的温度,/, /,/,m, ,3,4.14.44.34. 2总目录总目录4.54.1微分方程在化工中的应用微分方程在化工

4、中的应用 另一个在化工中常见的微分方程是物料冷却过程的数学模型，其模型可用下式表示： 它含有自变量t（时间）、未知函数 t（随时间变化的物料温度）、t0（环境温度）、k（降温速率）以及温度的一阶导数 ,是一个常微分方程。在微分方程中我们称自变量函数只有一个的微分方程为常微分方程，自变量函数个数为两个或两个以上的微分方程为偏微分方程。给定微分方程及其初始条件，称为初值问题；给定微分方程及其边界条件，称为边值问题。 在化工模拟中主要碰到的是常微分方程的初值问题： 或记为 只有一些特殊形式的 , 才能找到它的解析解；对于大多数常微分方程的初值问题, 只能计算它的数值解。常微分方程初值问题的数值解就是

5、求y(x)在求解区间a,b上各个分点序列xn，n =1,2,m的数值解yn。在计算中约定y(xn)表示常微分方程准确解的值，yn表示y(xn)的近似值。 下面我们将向大家介绍几种常用的常微分方程数值求解法。 )(0ttkdtdtdtdt)( ,)(),()( 0bxayayyxfxy)( ,)(),(0bxayayyxfdxdy),(yxf4.14.44.34. 2总目录总目录4.54. 2 欧拉（欧拉（euler）公式）公式 4.2.1向前欧拉公式向前欧拉公式 4.2.2向后欧拉公式向后欧拉公式 4.2.3中心欧拉公式中心欧拉公式 4.2.4梯形公式 4.14.44.34. 2总目录总目录4

6、.54.2.1向前欧拉公式向前欧拉公式 对于常微分方程初值问题式（4-5）,在求解区间a,b上作等距分割，步长 ，记xn = xn-1+h , n =1,2,m。用差商近似导数计算常微分方程。 做 的在x = x0处的一阶向前差商得： 又 ，于是得到： 故y(x1)的近似值y1可按 求得。类似地，由 得到计算近似值的向前欧拉公式： mabh)(xyhxyxyxy)()()( 010)(,()( 000 xyxfxy)(,()()(0001xyxfhxyxy),(),(00010001yxhfyyyxfhyy或)(,()( )()()( 1nnnnnnxyxfxyhxyxyxy以及),(1nnn

7、nyxhfyy 由yn直接算出yn+1值的计算格式称为显式格式，向前欧拉公式是显式格式。而欧拉方法的几何意义是以y1作为斜率，通过点(x0, y0)做一条直线，它与直线x=x1的交点就是y1。依此类推，是以yn+1作为斜率，经过点的直线与直线x=xn+1的交点。故欧拉法也称为欧拉折线法，如右图所示。 图4.1 欧拉折线法4.14.44.34. 2总目录总目录4.54.2.1向前欧拉公式向前欧拉公式 例例4.1：假定某物体的温度w因自热而产生的热量可以使物体在每秒钟内以4%的速度增长，同时该物体由于散热可使其温度在每秒种内下降100k，则物体温度随时间变化的微分方程： （t以秒为单位) 分别以初

8、始温度x(0)=1500k,y(0)=2500k,z(0)=3500k，用欧拉公式预测24秒后的物体温度趋势。 解： w0分别以x0=1500,y0=2500,z0=3500代入，计算结果见表4-1(下页)。 从表4-1可以看到当自热引起物体温度升高的速度小于散热引起温度下降的速度，物体的温度随时间而逐渐减少：当自热引起物体温度升高的速度与散热引起温度下降的速度平衡时，物体的温度保持不变；当自热引起物体温度升高的速度大于散热引起温度下降的速度，物体的温度随时间而增长。在图4-3（下页）中l1,l2,l3l1,l2,l3分别表示初始值3500,2500和1500的三条温度变化趋势曲线。 1000

9、4.0wdtdw1,10004. 1)10004. 0(1hwwhwwnnnn4.14.44.34. 2总目录总目录4.54.2.1向前欧拉公式向前欧拉公式n xn yn zn n xn yn zn 1 2 3 4 5 6 7 8 9 10 11 12 1460 1418.4 1375.14 1330.14 1283.35 1234.68 1184.07 1131.43 1076.69 1019.76 960.546 898.968 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 3540 3581.6 3624.86 3

10、669.86 3716.65 3765.32 3815.93 3868.57 3923.31 3980.24 4039.45 4101.03 13 14 15 16 17 18 19 20 21 22 23 24 834.926 768.324 699.056 627.019 552.1 474.183 393.151 308.877 221.232 130.081 35.2845 -63.3042 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 2500 4165.07 4231.68 4300.94 4372.98 4447.9

11、 4525.82 4606.85 4691.12 4778.77 4869.92 4964.72 5063.3 05101520250500100015002000250030003500400045005000wt图图4-3 三种初始值的温度变化曲线三种初始值的温度变化曲线 表表4-1 4.14.44.34. 2总目录总目录4.54.2.2向后欧拉公式向后欧拉公式 做出的在x1处的一阶向后差商式： 而， 得到的近似值y1的计算公式： 类似地，可得到计算y(xn+1)近似值yn+1的计算公式： 公式(4-7)称为向后欧拉公式。通常 为非线性函数，因此式(4-7)是关于yn+1的非线性方程，称为

12、隐式欧拉公式，需要通过迭代法求得yn+1。其中初始值 可由向前欧拉公式提供。 最简单的迭代公式为： 可以证明，h充分小时，以上迭代收敛。事实上，记 ，则 h充分小时，可以保证 ，其中l为李普希兹条件。 hxyxyxy)()()( 011)(,()( 111xyxfxy),(1101yxhfyy),(111nnnnyxhfyy),(yxf)0(1ny给定精度直到)(1)1(1)(11)1(1)0(1, 2 , 1 , 0,),(),(knknknnnknnnnnyykyxhfyyyxhfyy),()(1yxhfyynn),()( 1yxhfyny1),(1hlyxhfny4.14.44.34.

13、2总目录总目录4.54.2.3中心欧拉公式中心欧拉公式 做出y(x)的在x=x1处的中心差商式： 又 ，可得到y(x2)的近似值y2的计算公式： 类似地，可得到计算y(xn+1)近似值yn+1的计算公式： (4-8) 公式(4-8)称为中心格式。按公式(4-8)，需要知道yn-1, yn的值才能求得yn+1的值。因此，要先用其它公式计算出y1，再用中心格式算出y2, y3,。y1可用向前欧拉公式计算，为提高精度，也可用向后欧拉公式计算。 hxyxyxy2)()()( 021)(,()( 111xyxfxy),(21102yxhfyy),(211nnnnyxhfyy4.14.44.34. 2总目

14、录总目录4.54.2.4 梯形公式 在 两点之间进行梯形近似计算有： 则得梯形公式： 梯形公式也是隐式格式，计算中为了保证一定的精确度，又避免用迭代过程不菲的计算量，可先用显式公式算出初始值，再用隐式公式进行一次修正。称为预估-校正过程。例如，下面是用显式的欧拉公式和隐式的梯形公式给出的一次预估-校正公式： 上式也称为改进的欧拉公式，它可合并成 如果想要获得较高的计算精度，可进行多次迭代计算，也就是进行多次校正计算。下面的例子对每一个点进行了4次迭代计算。 1,nnxx)(,()(,(2 )()()(21)(11111nnnnnnnnxxxyxfxyxfhxyxyxxdxxynn),(),(2

15、111nnnnnnyxfyxfhyy),(),(2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy),(,(),(211nnnnnnnnyxhfyxfyxfhyy4.14.44.34. 2总目录总目录4.54.2.4 梯形公式 例例4.2：请用预估-校正公式（改进的欧拉公式）解下面初值问题： 解： 。 用下面的迭代公式，对每个点迭代4次，k=1,2,3,4。 该方程的精确解是 ，计算结果如表4-2所示。 0.40.0 , 1)0(2xyydxdy1 . 0, 10hy)(22)(12)1(12)0(1knnnknnnnyyhyyhyyyxy11n nx ny )(nxy )(n

16、nxyy 1 2 3 4 0.1 0.2 0.3 0.4 1.1118 1.2520 1.4311 1.6763 1.1111 1.2500 1.4326 1.6667 0.0007 0.0020 0.0095 0.0004 4.14.44.34. 2总目录总目录4.5vb程序程序比较x 向向前前y 向向后后y 中中心心y 1 0 0 0 1.1 1.1 1.1 1.10040916980883 1.2 1.20040950743689 1.20033346557798 1.20081833881153 1.3 1.30174836142175 1.30141370547213 1.30308

17、460881526 1.4 1.40419578896213 1.40336054223952 1.40569901443133 1.5 1.50769631079679 1.50609021418897 1.51006128256494 1.6 1.61203615561277 1.60938247891755 1.61432676578844 1.7 1.71690172242177 1.71293283220273 1.71972904125822 1.8 1.82192572264155 1.81639394221938 1.82428092462869 1.9 1.926724060

18、53033 1.91940879367697 1.92923501526952 2.0 2.03092532750262 2.02163718096732 2.03257402533708 4.14.44.34. 2总目录总目录4.54.3 龙格龙格-库塔方法库塔方法 龙格-库塔法是求解常微分方程较常用的一种方法，它通过巧妙的线性组合，在显式格式的情况下获得理想的计算精度，大大提高了计算速度。该方法的推导过程不是本课程要研究的对象，本课程主要研究其实际应用,，故直接给出各类龙格-库塔公式。 1、 二阶龙格-库塔 其中c1=1/2, c2=1/2, a=1, b=1或 其中c1=0,c2=1,a

19、=1/2,b=1/2。 二阶龙格-库塔公式的局部截断误差为o(h3), 是二阶精度的计算公式。也可建立高阶的龙格-库塔公式，如三阶龙格-库塔公式、四阶龙格-库塔公式。较常用的是四阶龙格-库塔公式，四阶龙格-库塔公式的局部截断误差界为o(h5), 是四阶精度的计算公式。 ),(),()(2121211hkyhxfkyxfkkkhyynnnnnn)2,2(),(12121khyhxfkyxfkhkyynnnnnn4.14.44.34. 2总目录总目录4.54.3 龙格龙格-库塔方法库塔方法2、三阶龙格、三阶龙格-库塔公式库塔公式 )2,()2,2(),()4(6) 1 (2131213211hkh

20、kyhxfkkhyhxfkyxfkkkkhyynnnnnnnn2312131132,3231,31,34 )2(hkyhxfkhkyhxfkyxfkkkhyynnnnnnnn23121321143,4321,21,4329 ) 3(hkyhxfkhkyhxfkyxfkkkkhyynnnnnnnn4.14.44.34. 2总目录总目录4.54.3 龙格龙格-库塔方法库塔方法3、四阶龙格、四阶龙格库塔公式库塔公式 342312143211,21,2121,21,226 ) 1 (hkyhxfkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnn321421312143211,31

21、,3231,31,338 )2(hkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnn4.14.44.34. 2总目录总目录4.54.3 龙格龙格-库塔方法库塔方法例例4.3：用四阶龙格库塔公式(4-16)求解下面初值问题：解：取步长h=0.2，计算公式为： 计算结果列表4-3中。 8 . 01 . 0 , 1) 0 (cos2xyxydxdy2 . 0cos2 . 01 . 0cos1 . 01 . 0cos1 . 0cos2262 . 02342232122143211nnnnnnnnnnxkykxkykxkykxykkkkkyyn nx ny n

22、xy nnxyy 1 2 3 4 0.2 1.24789 1.24792 0.00003 0.4 1.63762 1.63778 0.00016 0.6 2.29618 2.29696 0.00078 0.8 3.53389 3.53802 0.00413 4.14.44.34. 2总目录总目录4.54.3 龙格龙格-库塔方法库塔方法步长的选择步长的选择 怎样选取合适的步长，这在实际计算中是很重要的。由于步长愈小，每步计算的截断误差就愈小；但在一定的求解范围内，需要完成的步数就愈多，这不但引起计算量的增加，而且计算步数的增加往往造成舍入误差的严重积累，反而降低了计算精度。上面介绍的龙格-库塔方

23、法是对定步长（即步长h为常数）而言的，但为了保证精确度，一种有效的措施是在计算过程中自动进行步长调整，即变步长技巧。 下面以四阶龙格-库塔方法为例，说明如何自动选择步长，使计算结果满足给定精度的要求。 设从节点xn出发，先以h为步长，利用四阶龙格-库塔公式方法经过一步计算得y(xn+1)的近似值，记为 ，由于公式的局部截断误差是y(h5)，故有 当h不大时，c可近似地看作常数。然后将步长h对折，即取h/2为步长，从出发经过两步计算求y(xn+1)的近似值，记为 ，每一步计算的局部截断误差为c(h/2)5，于是就有 )(1hny5)(11)(chyxyhnn)2/(1hny5)2/(11)2/(

24、)(hcyxyhnn4.14.44.34. 2总目录总目录4.54.3 龙格龙格-库塔方法库塔方法步长的选择步长的选择 把它与（4-18）式相比，可得 经整理可得 这表明以 作为y(xn+1) 的近似值，其误差可用先后两次计算结果之差来表示，因而，只需考察 是否成立。若成立，则可将 作为y(xn+1)的近似值；若不成立，则将步长再次对折进行计算，直到不等式成立为止，并取最后的 作为计算结果。以上方法就是计算过程中自动选择步长的方法，也称为变步长方法。从表面上看，为了选择适当的步长，每一步的计算量增加了，但从总体考虑，工作量往往还是减少的。 龙格-库塔方法的主要优点是计算精确度较高，能满足通常的

25、计算要求；容易编制程序；每次计算y(xn+1) ，只用到前一步的计算结果yn，因此，在已知初始值y0的条件下，就可自动地进行计算，是单步法，而且计算过程中可随时改变步长。它的缺点是每前进一步需要计算多次f(x,y)的值，因此，计算工作量较大，且其截断误差难以估计。 在实际应用上，一般当要求更高精确度时，采用的办法是缩小步长，而不是采用更高阶的公式，因为高阶公式的计算太复杂，一般选用标准四阶龙格-库塔方法即可。 161)()()(11)2/(11hnnhnnyxyyxy)()()(11151)2/(11hnnhnnyxyyxy)2/(1hny)2/(11)(hnnyxy)2/(1hny)2/(1

26、hny4.14.44.34. 2总目录总目录4.54.3 龙格龙格-库塔方法库塔方法步长的选择步长的选择下面用一个例子说明： 由vb选用标准四阶龙格-库塔方法计算得 由积分法算得y(2)=2.23607,当h=0.5时相差0.00013，而h=0.25误差小于0.000001，但当 h=0.125时虽然足够精确但计算次数却比h=0.25多了一倍。因此合理选择步长既能保证精度又能减少计算量。 1)0( yyxdxdy步长 y(2) 计算次数 h=1 2.23940 2 h=0.5 2.23620 4 h=0.25 2.23607 8 h=0.125 2.23607 16 4.14.44.34.

27、2总目录总目录4.54.4 常微分方程组的数值解法常微分方程组的数值解法 4.4.1 一阶常微分方程组的数值解法一阶常微分方程组的数值解法 4.4.2 高阶常微分方程数值方法高阶常微分方程数值方法 4.14.44.34. 2总目录总目录4.54.4.1 一阶常微分方程组的数值解法一阶常微分方程组的数值解法将由m个一阶方程组成的常微分方程初值问题： 写成向量形式： 其中： btaayayayyyytfdtdyyyytfdtdyyyytfdtdymmmmmmm)()()(),(),(),(22112121222111)(),(ayytfdtdymmmmmmyytfyytfyytfytftytyty

28、ty211121121,),(),(),(),(,)()()()(4.14.44.34. 2总目录总目录4.54.4.1 一阶常微分方程组的数值解法一阶常微分方程组的数值解法 前面对常微分方程所用的各种方法，都可以平行地应用到常微分方程组的数值解中。下面以两个方程组为例，给出相应的计算公式。 常微分方程组： 欧拉公式： 预估校正公式： btazazyayzytfdtdzzytfdtdy )()(),(),(00),(),(11nnnnnnnnnnzythgzzzythfyy),(),(),(),(2),(),(1111111111nnnnnnnnnnnnnnnnnnnnnnnnnnzytgzy

29、tfzytgzytfhzyzyzytgzytfhzyzy4.14.44.34. 2总目录总目录4.54.4.1 一阶常微分方程组的数值解法一阶常微分方程组的数值解法四阶龙格库塔公式： )2(3)1(3)2(3)1(3)2(4)1(44)2(2)1(2)2(2)1(2)2(3)1(33)2(1)1(1)2(1)1(1)2(2)1(22)2(1)1(11)2(4)1(4)2(3)1(3)2(2)1(2)2(1)1(11143211,2,2,22,2,22,2,22,2,2),(),(226226hkzhkyhtghkzhkyhtfkkkkhzkhyhtgkhzkhyhtfkkkkhzkhyhtgk

30、hzkhyhtfkkkzytgzytfkkkkkkkkkkkhzyzykkkkhyynnnnnnnnnnnnnnnnnnnnnnnnnnnnnn4.14.44.34. 2总目录总目录4.54.4.1 一阶常微分方程组的数值解法一阶常微分方程组的数值解法 例例4.4：两种微生物，其数量分别是u=u(t)，v=v(t)，t的单位为分，其中一种微生物以吃另一种微生为生，两种微生物的增长函数如下列常微分方程组所示，预测3分钟后这一对微生物的数量。 解：记 在本题中f(t,u,v)=f(u,v)，g(t,u,v)=g(u,v)。 用欧拉预估校正公式(4-22)： 取，计算结果见表4-4。 2.1)0(6.1)0(001.015106.045.020109.0vuuvvvdtdvuvuudtduuvvvvuguvuuvuf001. 015