第二章数值积分法2011

第二章 连续系统仿真的数值积分法
2.1数值积分法 2.2欧拉法（Euler） 2.3龙格-库塔法（Runge-Kutta） 2.4亚当姆斯法（Adams） 2.5吉尔法（Gear） 2.6方法与步长的选择
2.1数值积分法
一、在仿真中研究数值积分法的意义
数值积分是数值分析的一个基本问题，用于求一个积 分的具体数值。 在数学上，用数值积分法求微分方程的初值问题。
y (tn 1 ) y (tn h) y (tn ) hy (tn ) h ! (tn ) ... h! y ( r ) (tn ) O(h r 1 ) y 2 r
2 r
问题：高阶导数计算不便 解决：受改进Euler法（将两点的一阶导数线性组合，从而使
精度提高）启发，将若干点的一阶导数值进行线性组
合，代替高阶导数，达到提高精度的目的。
一、龙格-库塔法的基本原理
对公式
y(tn1 ) y(tn )
t n1
tn
f (t , y)dt
(1)
若令
则有
yn y(tn )
Qn
t n 1
tn
f (t , y)dt
y(t n1 ) yn1 yn Q n
(2)
Qn的求解决定了数值积分法的计算公式。
yn 1 yn f ( , y ) d
tn
t n1
f(t,y) fn+1 fn
可认为 f（，y）= f（tn，yn）
t
用矩形面积 代替曲线下 的面积
0
tn
tn+1
若区间[ tn，tn+1 ]足够小
则有
yn+1 = yn + h f （tn，yn） h = tn +1 - t n
3
校正公式反复迭代的计算公
式
( yn0)1 yn hf (tn , yn ) ( ( yn1)1 yn h [ f (tn , yn ) f (tn 1 , yn0)1 )] 2 ( ( yns 1) yn h [ f (tn , yn ) f (tn 1 , yns )1 )] 1 2
截断误差正比于h3，称为二阶龙格-库塔法（简称RK-2）。 按照这一思想，可构造出截断误差正比于h5的四阶 龙格-库塔法（简称RK-4）。
二、固定步长的龙格库塔法
z (1 h) 0 z 1 h
为保证差分方程的解稳定，其特征根应位于单位圆内，即
z 1 h 1
解得：
2 h 0
可见，数值积分法的稳定性与计算步长h、原系统特征值有关。
2.3 龙格库塔法（Runge-Kutta）
提高数值积分法的精度===> 提高截断误差的阶数r ===> 在台劳展式中多截取几项===>

5、数值解的稳定性问题
用数值积分法仿真，所得仿真模型（即递推公式）应
保持原系统的稳定性。递推公式实际上是一个差分方
程，故我们可用差分方程的特征根判断其稳定性。
用Euler法求解检验方程：
y y, y(t0 ) y0
递推公式为： yn1 yn hf n yn yn (1 h) yn 这是一个差分方程，其特征方程和特征根分别为：
f(t,y) fn+1 fn
h为步长
欧拉公式
计算较简单。
用矩形面积代替曲线下 的面积，有较大误差。
一般不单独使用
t 0 tn tn+1
多与隐式公式配合使用，用来求初值。
例一 解:
. 已知 y’ + y² = 0，y(0) = 1，步长h=0.1 用欧拉法求数值解，并与解析解比较. 解析解为 y(t)= 1 / (1+t)， 数值解为 yn+1= yn+hf(tn,yn)，y’ = -y2 t1=0.1 t2=0.2 t3=0.3 t4=0.4 t5=0.5 t6=0.6 t7=0.7 t8=0.8 t9=0.9 t10=1 y(t1)=0.9091 y(t2)=0.8333 y(t3)=0.7692 y(t4)=0.7143 y(t5)=0.6667 y(t6)=0.6250 y(t7)=0.5882 y(t8)=0.5556 y(t9)=0.5263 y(t10)=0.5 y1=1+0.1(-1² )=0.9 y2=0.9+0.1(-0.9² )=0.819 y3=0.819+0.1(-0.819² )=0.7519 y4=0.6954 y5=0.6470 y6=0.6052 误差的数 y7=0.5685 量级在 y8=0.5362 10¯ , 即 h² ² y9=0.5075 y10=0.4817
改进欧拉法
例二. y´+ y²=0,

y(0)=0,
h=0.1,改进欧拉法
0.410-3 0.710¯ ³ 0.810¯ ³ 0.810¯ ³ 0.810¯ ³ 0.810¯ ³ 0.810¯ ³ 0.710¯ ³ 0.710¯ ³ 0.710¯ ³
在t0附近将y展开为台劳级数，只保留 h2 项，则有： 1 f dy f y1 y0 f (t0 , y0 )h ( ) t h2 (3) 2 y dt t
0
（3）式可以写成如下形式：
y1 y0 (a1 k1 a2 k 2 )h
其中
k1 f (t0 , y0 )

二.改进欧拉法
yn 1 yn f ( , y ) d
tn t n1
f(t,y) fn+1 fn
用梯形面积代替曲线 下的面积,可提高精度。
t 0 tn tn+1
用梯形面积 代替曲线下 的面积
yn+1= yn + h/2[ f(tn,yn)+ f(tn+1,yn+1)]
yn+1=yn+h/2 (3fn-fn-1)
3.截断误差
用台劳级数作为分析数值积分法精度的工具。 设yn是精确的，将y(tn+1)在tn点展开成台劳级数。
y(tn 1 ) y(tn h)
一阶精度 截断误差为O(h2)
h2 y(tn ) hy(tn ) 2!
h r y ( r ) (t ) O(h r 1 ) (tn ) ... r! y n
欧拉法
省略
每一步由这种方法引入的误差称为局部截断误差，简称截 断误差。 不同的数值积分法有不同的截断误差。一般，若截断误差 是O(hr+1)，则方法为r阶的。所以方法的阶数可以作为衡量 方法精度的一个重要指标。
改进欧拉法由于采用了平均斜率，相当于取台劳级数的前3 项，故为二阶精度。
4、舍入误差
计算机字长有限 ==> 对计算结果进行舍入处理==> 引起舍入误差 计算步数越多，舍入误差越大。
y f (t , y), y(t0 ) y0
（1）
对于一个连续系统，其模型可用微分方程（组）来表 示，而对其进行仿真即是从初始条件出发，求其过渡过程。 这在本质上是微分方程的初值问题。
二、数值积分法的基本原理
对微分方程（1）两边求积分可得下式：
y (t ) y (t0 ) f ( , y )d
0
(5)
将(5)式与(3)式进行比较，可得如下三个方程： a1+a2=1 a2b1=1/2 a2b2=1/2 令a1=a2，则可 解出这四个参 数。 a1=1/2 a2=1/2 b1=1 b2=1
计算公式为：
h (k1 k 2 ) 2 k1 f (t0 , y0 )，k 2 f (t0 h，y0 k1h) y1 y0
数值解 y1=0.9095 y2=0.8340 y3=0.7700 y4=0.7151 y5=0.6675 y6=0.6238 y7=0.5890 y8=0.5563 y9=0.5270 y10=0.5007
可见，其精度高于欧拉法。
若要得到更高的精度，可将校正公式反复迭代。
误差的数 量级在 3 10¯ , 即 h
三、几个基本概念
1.显式方法与隐式方法
公式右端无待求量的方法称为显式方法，如欧拉法。
公式右端有待求量的方法称为隐式方法，如改进欧拉法中的 梯形公式。隐式方法不能自启动，需要采用预估校正法。
2.单步法与多步法
求解tn+1时刻的值只用到tn时刻的值，这种方法称为单步法。
若求解tn+1时刻的值不仅用到tn时刻的值，而且需要tn-1、 tn-2… 时刻的值，则这种方法称为多步法。
k 2 f (t 0 b1h，y0 b2 k1h)
(4)
对k2式右端的函数展成台劳级数，保留h项，可得：
k 2 f (t0 , y0 ) (b1 f f b2 k1 ) t0 h t y
将k1、k2带入（4）可得：
y1 y0 a1hf (t 0 , y0 ) a2 h[ f (t 0 , y0 ) (b1 f f b2 k1 ) t h] t y
精度略 有提高
例三.
y´+ y² =0，y(0)=0，h=0.1， y1= 0.9086 取 s=1 ， 即 迭 代 y2= 0.8334 一次，得到计算 y3= 0.7695 结果列于右侧。 y4= 0.7148 y5= 0.6672 y6= 0.6256 实际应用时要综 y7= 0.5889 合考虑精度和计 y8= 0.5562 算量等问题。 y9= 0.5269 y10= 0.5006
e1= -4.5042e-004 e2= 3.2511e-005 e3= 3.1396e-004 e4= 4.7607e-004 e5= 5.6580e-004 e6= 6.1087e-004 e7= 6.2814e-004 e8= 6.2806e-004 e9= 6.1726e-004 e10= 5.9997e-004
2.2欧拉方法（Euler）
一、 欧拉法
欧拉法是最简单的数值积分法，常用来说明数值积分 法的几何意义和基本概念。
y f (t , y),
y(t0 ) y0
f n+1 fn
f(t,y)
yn 1 yn f ( , y ) d
tn
t n1
t 0 tn tn+1
曲线下的面积
箭头所指的部分结构形式相同，只是积分上限不同。
故有：
y (t n 1 ) y (t n ) f ( , y )d
tn
t n1
若用近似值yn+1、yn、Qnwk.baidu.com替上式中的各项，则可写出
yn+1 = yn+Qn
t=t0即 n=0时，y1=y0+Q0，其中， y0已知，若用数值积分法 求出Q0即可得到y1。依此类推，可求出t1 、t2、… 、tn时刻的 近似值 y1、y2、 … 、 yn 可见，这是一种步进式的近似解法，故存在精度和稳 定性问题。 用不同的方法求积分值 Qn就形成了不同的数值积分法， 每种方法的精度、稳定性和计算速度各不相同。 我们希望的数值积分法：精度高，速度快，稳定性好。
解析解 y=1/(1+t) t1=0.1 y(t1)=0.9091 t2=0.2 y(t2)=0.8333 t3=0.3 y(t3)=0.7692 t4=0.4 y(t4)=0.7143 t5=0.5 y(t5)=0.6667 t6=0.6 y(t6)=0.625 t7=0.7 y(t7)=0.5882 t8=0.8 y(t8)=0.5556 t9=0.9 y(t9)=0.5263 t10=1 y(t10)=0.5
yn+1= yn + h/2[ f(tn,yn)+ f(tn+1,yn+1)]
由于梯形公式(2.2)的右端有待求量，不能自启动， 故需要与欧拉法结合使用，这就是改进欧拉法，实 际上这是一种预估 —校正法。 预估 校正 yn+1 = yn + hf(tn,yn) yn+1 = yn + h/2[f(tn,yn)+f(tn+1,yn+1)]
t0 t
即：
y (t ) y (t0 ) f ( , y )d
t0
t n1 t0 tn t0 t n1 tn
t
当t=tn+1时，上式为： y (t n 1 ) y (t0 ) f ( , y )d 将积分项分为两项：
y (t n 1 ) y (t0 ) f ( , y )d f ( , y )d