第三章 弹道数值算法

式反复迭代。通常在工程问题中，为简化计算，只迭代一次。这样可得改进的欧拉公式
34
ynp1 yn h f ( yn , tn ) h f ( yn , tn ) f p ( ynp1 , tn 1 ) yn 1 yn 2
(3.8)
上式中第一式称为预估公式，第二式称为校正公式。通常称这类方法为预估—校正方法，也 称为改进的欧拉法。欧拉法每计算一步只要对导数 f 调用一次，改进的欧拉法由于加了校正 过程，计算量增加了一倍，付出这种代价的目的为了提高精度。其算法框图如图 3-4 所示。
y (t1 ) y0 f ( y, t )dt
t0 t1
(3.3)
式(3.3)中的积分项是曲线 f 及 t t0 , t1 包围的面积（如图 3-1 所示） ，当步长 h tn t n 1 足 够小时，可以用矩形面积来近似，即 y (t1 ) y (t0 ) f ( y0 , t0 )(t1 t0 ) 令 y (t1 ) 的近似值为 y1 ，则有
h yn 1 yn (k1 2k2 2k3 k4 ) 6 k1 f ( yn , tn ) h h k2 f ( yn k1 , tn ) 2 2 h h k3 f ( yn k 2 , t n ) 2 2 k4 f ( yn hk3 , tn h)
37
此为 4 阶五级公式，还可推导出一个 3 阶四级公式
h yn 1 yn (3k1 9k3 12k4 ) 6
(3.13)
令误差 则
En
En yn 1 yn 1
h (2k1 9k3 8k4 k5 ) 6
(3.14)
式(3.12)、(3.13)、(3.14)简称为RKM3-4 法。根据该步的绝对误差 En ，既可按步长的控制 策略进行步长的控制。通常用对分策略。 设定一个最小误差限 min ，一个最大误差限 max ，每一步的局部误差取为
en En / ( yn 1)
(3.15)
其中 En 为变步长各式计算出的误差估计。由式(3.15)可知，当 yn 较大时， en 是相对误差，而 当 yn 的绝对值很小时， en 就成了绝对误差，这样可避免当 y 值很小时， en 变得过大。 其控制策略是：当 en 大于 max 时，将步长对分减半，并重新计算该步；当 en 在 max 与 min 之间时，步长不变；当 en 小于 min 时，将步长加倍。即 1 若 en emax ，则 hn hn ，重算此步； 2 若 emin en emax ，则 hn 1 hn ，继续计算； 若 en emin ，则 hn 1 2hn ，继续计算。 这种对分策略简便易行，每步附加计算量小，但不能达到每步最优。还有一种最优步长 ， 控制策略，其基本思想是在保证精度的前提下，每个积分步取得最大步长（或称最优步长） 这样可以减少计算量。具体做法是根据本步误差估计，近似确定下一步可能的最大步长。其 策略如下。 给定相对误差限 0 ，设本步步长为 hn ，本步相对误差估计值为 En en yn 1 对于 k 阶积分算法，认为
第三章 仿真数值算法
仿真算法是将系统数学模型转换成适于计算机仿真模型的一类算法。本章主要介绍了工 程领域最常见的连续系统的仿真算法，包括数值积分法、插值法。Equation Section 3
3.1 数值积分法
数值积分法就是对常微分方程组建立离散形式的数学模型——差分方程，并求出其数值 解。例如已知某系统的一阶向量微分方程为 f ( y, t ), y (t0 ) y0 y (3.1)
(3.9)
取泰勒展开式的前两项可以得到欧拉公式，其误差较大，如果要得到精度更高的近似解，必 须计算式中的高阶导数，这项工作往往相当困难。德国数学家 C.Runge 和 M.W.Kutta 两人先 后提出了间利用泰勒展开式的方法， 即用几个点上的函数 f 值的线性组合来确定其中的系数。 基于这一思想，得到龙格库塔(RK)法的一般形式为
方法名称 欧拉法 折线法 2 阶 RK 法 3 阶 RK 法
阶次 1 2 2 3
ij 21 1 21 1
21
1 2 31 1, 32 2 1 2 1 2
wi w1 1 w1 0, w2 1 w1 w2 1 2
1 2 w1 w3 , w2 2 3
en 1
(tn )hnk1
与式(3.8)相比，二者完全相同，所以预估—校正公式实际上是二阶RK公式。当阶次大于 4 阶 后，龙格库塔公式右端函数的计算次数要大于阶数，使积分工作量大大增加，所以通常只使 用 4 阶或 4 阶以下的方法。表 3-1 给出了常用的 1~4 阶龙格库塔方法的系统。 4 阶龙格库塔法 是使用较多的一种方法，其公式如下
y1 y0 hf ( y0 , t0 )
(3.4) (3.5)
把 t1 当作积分初始点， y1 作为初始值重复上述做法，可进一步得到 y (t2 ) 的近似公式，继续重 复可得到递推公式
yn 1 yn hf ( yn , tn )
(3.6)
33
式(3.6)称为欧拉公式，也称为矩形法。由式(3.6)可以看出，任何一个新的数值解 yn 1 都是基于 前一个数值解 yn 以及它的导数 f ( yn , tn ) 求得的。若已知初值 y0 ，利用式(3.6)进行迭代计算， 既可以求得式(3.1)在 t t1 , t2 , tn 处的近似解 y (t1 ), y (t2 ), , y (tn ) 。
yn 1 ai y n i h i f
i 0 i 1 m m n i
(3.2)
3.1.1 几种常见的积分法
（1）欧拉法 欧拉法是最简单的一种数值积分法。虽然它的计算精度较低，实际中很少采用，但其推 导简单，能说明构造数值解法一般计算公式的基本思想。 对式(3.1)两端由 t0 到 t1 进行积分，得到
En ( ) hnk
(3.16)
其中 ( ) 是 f ( y, t ) 在积分区间 (tn ~ tn h) 内一些偏导数的组合，通常可取 tn ，则
en
38
(tn ) hnk
yn 1
(3.17)
据此作判断： （1）若 en 0 ，则本步积分成功，先确定下一步的最大步长 hn 1 。 假定 hn 1 足够小，则 (tn hn 1 ) (tn ) ，下一步误差为
所谓数值解法，就是寻求式(3.1)中 y 在一系列离散点 t1 , t2 , , tn 的近似解 y1 , y2 , , yn ，相 邻两个点之间 h tn tn 1 ，称为计算步长或步距。根据已知的初始条件 y0 ，采用不同的递推算 法可逐步递推计算出各时刻的数值 yi 。常用的方法有欧拉法、梯形法、四阶龙格库塔法、亚 当姆斯法等。对式(3.1)，数值积分可写成统一公式
y
y0
y1
误差
y2
f (t )
f
t0
t1
t2
t
Hale Waihona Puke Baidu
0
t0
t1
t
图 3-2 欧拉折线
图 3-3 梯形近似及其误差
（2）梯形法 在上面推导中，若用图中的梯形面积来近似(3.3)中的积分项，则可得到梯形公式
yn 1 yn h f ( yn , tn ) f ( yn 1 , tn 1 ) 2
(3.11)
龙格库塔法属于单步法，只有给定方程的初值 y0 就可以一步步求出 y1 , y2 , yn 的值。故 单步法由下列优点：1）需要存储的数据量少，占用的存储空间少；2）只需要知道初值，既 可启动递推公式进行计算，可自启动；3）容易实现变步长运算。
36
表 3-1 常用的 1~4 阶龙格库塔法的系数
1 6 1 w2 w3 3 w1 w4
21
4 阶 RK 法 4
31 0, 32 43 1
3.1.2 变步长法
在实际使用时，对前述数值积分方法，仿真人员可根据实际情况在仿真的不同阶段选取 不同的步长， 也就是变步长积分法。 例如变步长龙格—库塔—默森 （Runge-Kutta-Merson） 法。 选取的原则是，在保证仿真过程满足一定精度的前提下，为使计算尽可能小，尽量选取适用 的较大的步长，这样仿真步长需不断改变。变步长应根据一定的条件，其前提是要有一个好 的局部误差根据公式，根据局部误差的大小来改变步长。对于龙格库塔算法的误差估计，通 常是设法找到另一个低阶（一般是低一阶）的龙格库塔公式，要求这两个公式中的 ki 相同， 则两个公式计算结果之差可以看作是误差。 假设微分方程为
开成 h 的幂级数，然后和台劳展开式的系统进行对比，以确定 ij ， wi 的值。当 r 1 时，得到 的数值解即为欧拉公式
yn 1 yn hf ( y n , tn )
当 r 2 时，取 21 1 ，则 w1 w2 0.5 ，继而可得到
1 yn 1 yn 2 h(k1 k2 ) k1 f ( yn , tn ) k f ( y k h, t h) n 1 n 2
f (t )
f
t0
t1
t
图 3-1 矩形近似及其误差
欧拉法的几何意义十分清楚。 图 3-2 通过 (t0 , y0 ) 点作积分曲线的切线， 其斜率为 f ( y0 , t0 ) ， 此切线与 t1 处平行于 y 轴直线的交点即为 y1 ，再过 (t1 , y1 ) 点作积分曲线的切线，它与过 t2 平行 于 y 轴直线的交点即为 y2 。这样过 (t0 , y0 ) 、 (t1 , y1 ) 、 (t2 , y2 ) ， ，得到一条折线，称为欧拉 折线。
35
yn 1 yn h wi ki
i 1
r
(3.10)
式中
k1 f ( yn , tn ) ki f ( yn h ij k j , tn i h) (i 2,3, r )
j 1 i 1
i ij
j 1
i 1
i ， ij ， wi 为待定系数。 r 为使用 k 值得个数（即阶数） 。在给定 r 值后，通过把式(3.10)展
(t ) f ( y, t ) y y (t0 ) y0
计算公式
h yn 1 yn (k1 4k4 k5 ) 6
(3.12)
其中
k1 k 2 k3 k4 k 5 f ( yn , t n ) h h f ( y n k1 , tn ) 3 3 h h f ( y n (k1 k2 ), tn ) 6 3 h h f ( y n (k1 3k3 ), tn ) 8 2 h f ( y n (k1 3k3 4k4 ), tn h) 2
设置初值： y0 t 0 h
计算导数 fn
计算预估值 ynp1
tn1 tn h
p 计算 fnp 1 f yn1 , tn1
计算校正值 yn1
输出 t , y
图 3-4 预估—校正法程序框图
（3）龙格库塔法 将式(3.1)在 tn 点展成泰勒级数
(tn ) y (tn h) y (tn ) hy h2 y (tn ) o(h3 ) 2
(3.7)
由上式可见，它是隐函数形式。公式右端隐含有待求量 yn 1 ，故梯形法不能自启动。通常可用 欧拉法启动求出初值，算出 y (tn 1 ) 的近似值 ynp1 ，然后将其代入原微分方程，计算 f n 1 的近似
p 值 f np 1 f ( yn 1 , t n 1 ) ，最后利用梯形公式求出修正后的 yn 1 。为了提高计算精度，可用梯形公