常微分方程初值问题的高阶泰勒法与龙格-库塔法之应用对比

第39卷 第12期 高 师 理 科 学 刊 Vol. 39 No.12 2019年 12月 Journal of Science of Teachers′College and University Dec. 2019
文章编号：1007-9831（2019）12-0012-04
常微分方程初值问题的高阶泰勒法与
龙格-库塔法之应用对比
江山1，张岩1，孙美玲2
（1. 南通大学 理学院，江苏 南通 226019；2. 南通职业大学 公共教学部，江苏 南通 226007）
摘要：常微分方程的初值问题有着广泛应用，其数值解的精确化和高效化一直是人们追求的目标．基于严谨数学理论提升局部截断误差的阶数，从而达到提升整体截断误差阶数的目的，给出高阶泰勒法与龙格-库塔法的迭代公式．通过具体算例，比较各种高阶算法的计算精度和计算效率，给出对应的数值误差与图表说明，充分验证２类高阶方法各自的优势与缺陷，为求解实际应用问题提供参考依据．
关键词：常微分方程；初值问题；泰勒法；龙格-库塔法
中图分类号：O241.81文献标识码：A doi：10.3969/j.issn.1007-9831.2019.12.004
Application comparisons of higher order Taylor method and Runge-Kutta method for initial value problem of ordinary differentiable equations
JIANG Shan1，ZHANG Yan1，SUN Mei-ling2
（1. School of Science，Nantong University，Nantong 226019，China；
2. Department of Public Courses，Nantong Vocational University，Nantong 226007，China）
Abstract：The initial value problem of ordinary differentiable equation is universally used，its most accurate and most efficient numerical solution is always a goal．Based on mathematical theories to upgrade the order of local truncation error，consequently the order of global truncation error is upgraded too，the higher-order Taylor and Runge-Kutta iterative formula are presented．The computational accuracy and efficiency from different algorithms are compared through numerical experiments，and their corresponding errors are plotted and listed．In this way，sufficiently validates the advantages and disadvantages of these higher-order methods，and provides references for solving practical application problems．
Key words：ordinary differentiable equation；initial value problem；Taylor method；Runge-Kutta method
常微分方程是描述连续变化的数学语言，它是研究宏观和微观事物、自然现象、运动规律和变化过程的重要工具，在众多领域有着广泛的应用．不过并不是所有 解方程都具备解析解，对于某些复杂的常微分问题，其数值解更有实际意义[1-2]．而最大程度地追求高精度、高效率的数值结果，并将其应用至科学计算领域，一直是研究者孜孜追寻的目标[3-10]．
收稿日期：2019-07-18
基金项目：江苏省高校青蓝工程优秀青年骨干教师资助项目；南通大学教学改革研究课题（2018B06）
作者简介：江山（1980-），男，湖南湘潭人，教授，博士，从事微分方程数值解及其应用研究．E-mail：jiangshan@
第12期 江山，等：常微分方程初值问题的高阶泰勒法与龙格-库塔法之应用对比 13
1 初值问题与存在唯一性、初值敏感性
考虑常微分方程初值问题
[]()000(, ) , n y f x y x x x y x y '⎧=∈⎪⎨=⎪⎩
（1） 理1
[1]279
若(, )f x y 连续，且存在常数0L >，使得对于任意12, y y ∈R ，有()()12, , f x y f x y -≤
12L y y -，则称其满足李普希茨条件．此时，对于任意[]0, n x x x ∈，0y ∈R ，常微分初值问题（1）存在
唯一的连续可微解()y x ． 理2
[1]279
若(, )f x y 连续，且满足李普希茨条件，设常微分初值问题（1）的解与初值0y 有关，表
示为()0, y x y ，则有()()0
1212, , e
L x x y x y y x y y y --≤-．
理２表明 解对初值的依赖敏感性．当李普希茨常数L 较小时，解对初值和(, )f x y 相对不敏感，
可视为良态条件；当L 较大时，解对初值和(, )f x y 相对敏感，可视为病态条件．
2 高阶的泰勒法与龙格-库塔法
虽然求解常微分方程有解析方法，但往往只能求解一些特殊情形．实际生产生活中的常微分方程求解主要采用数值方法，即求解精确解()y x 在一系列等距离散节点01, , , n x x x L （1j j x x h +=+，h 为步长）对应的数值解{}n y ，作为精确解很好的近似．低阶方法有一阶的欧拉法和二阶的预测-校正法等，但其数值精度较低、误差较大．本文主要研究高阶方法，即高阶泰勒法和高阶龙格-库塔法，其数值精度更高，误差更小． 2.1 高阶泰勒法
运用高阶泰勒法，需要存在光滑的各阶导数(3)(1), , , , n f f f f -'''L ，理论推导上利用泰勒展开式，记()()()()1()
(1)
, , , , 2!!
n n n j j j j j j j j h h T x y f x y f x y f x y n --'=+++L ，则n 阶泰勒法公式为()
()1, n j j j j y y hT x y +=+．泰
勒法具有高阶的整体截断误差，但同时也要求(, )f x y 有较严格的光滑性，计算时各阶导数既要存在且连
续，还要容易求解，这会在应用上带来一 的局限性． 2.2 高阶龙格-库塔法
利用多元函数的泰勒展开，通过提升局部截断误差达1n +阶时，整体截断误差也提升为n 阶．再通过求解方程个数小于未知量个数的欠 方程组，可以得到满足限制条件约束的不同系数，从而构造出不同阶数的迭代公式．如二阶龙格-库塔法12j j y y hk +=+，其中：()
1, j j k f x y =； ()
210.5, 0.5j j k f x h y hk =++． 三阶龙格-库塔法()112346
j j h
y y k k k +=+
++，其中：()
1, j j k f x y =；()
210.5, 0.5j j k f x
h y hk =++；()312, 2j j k f x h y hk hk =+-+．四阶龙格-库塔法()11234226
j j h
y y k k k k +=++++，其中：()
1, j j k f x y =；
()
210.5, 0.5j j k f x h y hk =++；()320.5, 0.5j j k f x h y hk =++；()
43, j j k f x h y hk =++．
不同于泰勒法，龙格-库塔法不需要光滑的高阶导数并求其表达式，它仅需依次求解(, )f x y 在各点的函数值并对应记作i k ，再利用简单实用的线性组合来迭代计算即可．
可以看出，不论是高阶泰勒法，还是高阶龙格-库塔法，它们都是显式方法，即通过上一步第j 次迭代值即可算出下一步第1j +次迭代值，整体运算量和存储量都不大．
3 算例
分别用二阶、三阶、四阶泰勒法，二阶、三阶、四阶龙格-库塔法求初值问题21()(2)1y x y y ⎧'=+-⎪
⎨=⎪⎩
的数值
14 高 师 理 科 学 刊 第39卷
解，求解区间[]2, 3x ∈，取步长0.1h =，做10次迭代（该问题精确解为1()1y x x x
=+-）． 因2(, )1()f x y x y =+-，故()32()12()f x y y x y ''=--=--，则()
()()
(2), , , 2
j j j j j j h
T x y f x y f x y '=+=
()
()
2
3
1j j
j j x y h x y +---，故二阶泰勒法()
(2)1, j j j j y y hT x y +=+；同理类似，得三阶泰勒法的()
(3), j j T x y ，
四阶泰勒法的()
(4), j j T x y ，并代入高阶公式计算．
运用二阶、三阶、四阶龙格-库塔法分别求解1234, , , k k k k ，再作线性组合代入各自高阶公式计算．
结合计算机的快速运算能力，将推导转变为简便高效的Matlab 编程，求解结果见表1~2及图１~2．
表1 真解与各数值解
j x
真解 二阶泰勒法 三阶泰勒法 四阶泰勒法 二阶龙-库法 三阶龙-库法 四阶龙-库法 2.0 1.000 000 000 1.000 000 000 1.000 000 000 1.000 000 000 1.000 000 000 1.000 000 000 1.000 000 000 2.1 1.190 909 090 1.190 000 000 1.191 000 000 1.190 900 000 1.190 250 000 1.190 924 670 1.190 908 813 2.2 1.366 666 666 1.365 274 290 1.366 799 946 1.366 653 853 1.365 656 254 1.366 689 466 1.366 666 271 2.3 1.530 769 230 1.529 134 897 1.530 919 862 1.530 755 221 1.529 582 239 1.530 794 962 1.530 768 794 2.4 1.685 714 285 1.683 977 463 1.685 869 154 1.685 700 282 1.684 451 900 1.685 740 710 1.685 713 846 2.5 1.833 333 333 1.831 575 327 1.833 485 578 1.833 319 894 1.832 054 675 1.833 359 286 1.833 332 908 2.6 1.975 000 000 1.973 268 016 1.975 146 142 1.974 987 365 1.973 739 489 1.975 024 894 1.974 999 599 2.7 2.111 764 705 2.110 085 555 2.111 903 128 2.111 752 953 2.110 541 960 2.111 788 270 2.111 764 331 2.8 2.244 444 444 2.242 832 563 2.244 574 559 2.244 433 570 2.243 270 085 2.244 466 582 2.244 444 096 2.9 2.373 684 210 2.372 146 473 2.373 805 998 2.373 674 172 2.372 563 347 2.373 704 922 2.373 683 888 3.0
2.500 000 000
2.498 538 652
2.500 113 747
2.499 990 738
2.498 934 364
2.500 019 337
2.499 999 702
表2 各数值解的绝对误差
j x
二阶泰勒误差
（×10-2
） 三阶泰勒误差
（×10-3
） 四阶泰勒误差
（×10-4
） 二阶龙-库误差
（×10-2
） 三阶龙-库误差
（×10-4
） 四阶龙-库误差
（×10-6
） 2.0 0.000 000 000 0.000 000 000 0.000 000 000 0.000 000 000 0.000 000 000 0.000 000 000 2.1 0.090 909 090 0.090 909 090 0.090 909 090 0.065 909 090 0.155 799 242 0.277 241 310 2.2 0.139 237 666 0.133 279 334 0.128 135 330 0.101 041 228 0.227 998 421 0.395 509 739 2.3 0.163 433 305 0.150 632 045 0.140 096 745 0.118 699 156 0.257 316 214 0.436 522 516 2.4 0.173 682 227 0.154 868 530 0.140 037 415 0.126 238 515 0.264 252 047 0.439 607 047 2.5 0.175 800 618 0.152 244 973 0.134 385 059 0.127 865 745 0.259 535 653 0.424 399 196 2.6 0.173 198 343 0.146 142 733 0.126 343 158 0.126 051 029 0.248 943 949 0.400 949 173 2.7 0.167 915 064 0.138 422 327 0.117 526 820 0.122 274 552 0.235 643 920 0.374 463 192 2.8 0.161 188 066 0.130 114 864 0.108 744 127 0.117 435 929 0.221 384 110 0.347 627 658 2.9 0.153 773 735 0.121 788 301 0.100 384 207 0.112 086 293 0.207 123 534 0.321 790 574 3.0
0.146 134 707
0.113 747 501
0.092 616 195
0.106 563 581
0.193 374 085
0.297 580 231
第12期 江山，等：常微分方程初值问题的高阶泰勒法与龙格-库塔法之应用对比 15
由表1~2可以看出，同一类方法体系下，方法的阶数越高，数值解的精度越好，如泰勒法的绝对误差数量级由210-依次提升至310-，410-，龙格-库塔法的绝对误差由210-依次提升至410-，610-；而横向比较而言，阶数相同的龙格-库塔法又比泰勒法的精度更优，不需要(, )f x y 存在光滑的高阶导数且便于计算的限制，从而更具优势．图1~2再次验证，阶数更高的方法，精度越好，误差越小，而总体而言，同阶的泰勒法精度和误差又不如龙格-库塔法．
二阶泰勒法、三阶泰勒法、四阶泰勒法、二阶龙-库法、三阶龙-库法和四阶龙-库法的CPU 计算时间分别为0.002 893，0.008 727，0.009 436，0.011 298，0.011 324，0.012 209 s，观察到随着方法阶数的上升，高阶方法呈现出计算时间的线性增长（而龙格-库塔法呈现更为平缓的增长，不过其所需时间比同阶泰勒法稍多），这也再次证实 与理论知识相一致的数值结果．
本文从理论推导和数值实验２方面，给出并验证 高阶泰勒法与龙格-库塔法求解常微分方程初值问题数值解的精确性与稳 性，从计算精度、计算误差和计算效率等方面对比 ２类数值方法，得出 阶数更高的方法更有应用价值，阶数相同的方法更宜选用龙格-库塔公式的结论． 参考文献：
[1] 李庆扬，王能超，易大义．数值分析[M]．5版．北京：清华大学出版社，2008
[2] Butcher J C．Numerical Methods for Ordinary Differential Equations[M]．2nd ed．England：Wiley press，2008
[3]
Famelis I T，Tsitouras C．On modifications of Runge-Kutta-Nystrom methods for solving 4)()y f x y =(，[J]．Applied Mathematics and Computation，2016（273）：726-734
[4] Baeza A，Boscarino S，Mulet P，et al．Approximate Taylor methods for ODEs[J]．Computers and Fluids，2017（159）：156-166 [5] Higueras I，Roldan T．New third order low-storage SSP explicit Runge-Kutta methods[J]．Journal of Scientific Computing，2019，
（79）：1882-1906
[6] Vaquero J M，Kleefeld A．A fifth-order extrapolated stabilized explicit Runge-Kutta method[J]．Journal of Computational and
Applied Mathematics，2019（356）：22-36
[7] 孙美玲．常微分方程数值解法的Matlab 计算与可视比较[J]．高教学刊，2016（19）：60-61
[8] 马泽涛，文鹏，施琳达，等．基于Matlab 常微分方程数值解的分析与比较[J]．大学教育，2017（12）：50-52 [9] 周忠，周瑞芳．可视化方法在常微分方程数值解教学中的应用[J]．教育现代化，2018（1）：171-173 [10] 丛玉豪，胡洋，王艳沛．含分布时滞的时滞微分系统多步龙格-库塔方法的时滞相关稳 性[J]．计算数学，2019，41（1）：
104-112
x
x。