* 第十三章 数值计算初步第一节 误差与方程求根第二节 拉格朗日插值公式第三节 曲线拟和的最小二乘法第四节 数值积分第五节 常微分方程的数值解法一,误差二,方程求根第一节 误差与方程求根第一节 误差与方程求根
1.绝对误差与相对误差设 *x 为准确数 x 的近似值，我们称 xxxe **)( 为近似数 *x
的绝对误差，简称误差,由于一般无法得到准确值 x,因此绝对误差
)( *xe 也无法直接算出，如果能估计其绝对值的范围
.**)( xxxe则称? 为近似数 *x 的绝对误差限，简称误差限,
有了绝对误差限就可以知道准确值 x 的范围,
.xxxxx
一、误差即 x 落在闭区间 ],[
**
xx 上，在应用上，常采用如下写法来刻画
*
x 的精度
.
*
xx
绝对误差限并不能完全表示近似值的好坏程度，例如
110x
51 0 0 0y
虽然 x 的绝对误差限比 y 的绝对误差限小，但 1 0 0 0 作为 y
的近似值要比 10 作为 x 的近似值要好，为了清楚的描述这一现象，
我们引进
,)()( *** ** x xxxxexr
并称 )( *xr 为近似数 *x 的相对误差,
显然,)( *xr 的准确值也无法直接得到，如果我们知道
,)()( *
*
*
x
xexr
则称 r? 为近似数 *x 的相对误差限,
若? 为?x 的绝对误差限，则
*xr
前面提到的 110x 的近似值 10*?x 的相对误差限为
1 0 %,而 y = 1 0 0 0 ± 5 的近似值 1000*?y 的相对误差限是 0,5 %,
由此可见，相对误差限愈小，近似程度愈高,
2．有效数字定义 （有效数字） 若近似值 *x 的误差限是某一数位的半个单位，且该位到 *x 的左边第一位非零数字共有 n 位,则称
*x 有 n 位有效数字,
例 1 按四舍五入的原则写出下列各数具有 5 位有效数字的近似值,1 8 7,9 3 2 5,0,0 3 7 8 5 5 1,8,0 0 0 0 3 3,2,7 1 8 2 8 解 按定义，上述各数具有 5 位有效数字的近似值分别为,
187.93,0.037 855,8.000,2.718 3,按有效数字的定义，近似数的最后一个有效数位确实能反映其绝对误差的大小，并且可以证明：有效数字位数越多，相对误差限越小,
3,数值运算的误差估计设函数,,( 21 xxfu ),,nx 如果自变量,,21 xx? nx,
的某组值,,21 xx? nx,的近似值分别为,,*2*1 xx? *,nx,则函数
u 的近似值 fu?*,,( *2*1 xx? ),*nx,且 *u 的误差
),,,(),,,()( 21**2*1** nn xxxfxxxfuuue
).(d *
1
*
*
* kkxx
n
k k
xx
xxx xxx
fu
kk
kk
kkk
),()( *
1
*
* kxx
n
k k
xexfue
kk
即从而,*u 的误差限为 ).()( *
1
*
*
k
xx
n
k k
xxfu
kk
*u 的相对误差限为
.)()(
*1
*
*
*
*
kk xxk
n
k
k
r x
f
u
x
u
u
例 2 已测得某场地长 l 的近似值为 *l =110 m,宽
d 的近似值为 *d =80 m,已知
ll?* ≤ 0.2 m,dd?* ≤ 0.1 m,
试求面积 ldA? 的绝对误差与相对误差限,
解 ldA?
ldAdlA,
)()()( *** dlldA = 8 0 ×0,2 + 1 1 0 ×0,1
=27 )( 2m,031.0
80100
27)(
*
*
AAr
4,数值运算应注意的若干原则一个实际的计算问题，其数值运算次数以千万次计，我们不可能对其每步计算都进行误差分析，因此，在设计算法时要特别注意误差的积累,具体地说要注意如下几点,
⑴要尽量减少运算次数，以减少舍入误差的积累；
⑵要尽量避免两相近数相减，以防有效数字的严重损失而影响精度； ⑶要尽量避免使用“小分母”，以防运算结果过大而造成溢出；
⑷要防止大数“吃掉”小数现象,例如，在五位计算机上计算
.9.052492
1000
1
k
A
把运算写成规格化形式?
1000
1
55 10009000.01052492.0
k
A
555 100 0 0 0 0 9.0100 0 0 0 0 9.0105 2 4 9 2.0
51052492.0
该结果显然不正确，这是由于在计算机内计算时要对阶，每个 0,9 都被看成了零,要避免这一现象，就应先把数量级相同的数相加，然后在加上大数 52 4 92,也就避免了大数吃掉小数,
方程的根：方程 0)(?xf 的根也叫做函数 )( xf 的根，或叫做 )( xf 的零点,方程的根的存在性：若
)( xf 在闭区间ba,上连续，且
0)()( bfaf,则 )( xf 在开区间 ),( ba 内至少有一实根,
1,方程求根的二分法设 )( xf 为闭区间ba,上的连续函数，且在开区间
),( ba 内有惟一的单实根 ξ,下面，我们给出求单实根 ξ 的近似值方法,
二、方程求根取区间ba,的中点
2
0
ba
x
,考察区间 ],[
0
xa 及 ],[
0
bx
中哪一个为有根区间，即检查 )(
0
xf 与 )( af 是否同号,如果确系同号，说明所求的根 ξ 在区间 ],[ 0 bx 中，即新的有根区间为 ],[ 0 bx,这时令 bbxa 101,；否则 ξ 必在 ],[
0
xa 中，这时令 aa?1,01 xb?,即不管出现哪种情形，新的有根区间
],[
11
ba 的长度仅为原有根区间的一半,
如此反复二分下去，即可得出一系列有根区间 ],[ ba,
],[ 11 ba,],[ 22 ba,,? ],[ kk ba,
其中每一个区间都是前一个区间的一半，因此，二分 k
次后的有根区间 ],[ kk ba 的长度为对于新的有根区间 ],[ 11 ba 又可施行同样的作法，即用中点
2
11
1
ba
x
将区间 ],[ 11 ba 再分为两半，然后判定所求的根 ξ 在 1x
的哪一侧？从而确定出新的有根区间 ],[ 22 ba,其长度为区间
],[ 11 ba 的一半,
.2 kkk abab
由此可见,如果二分过程无限的继续下去，这些有根区间必收缩于点 ξ,该点就是所求的根,每次二分后，若取有根区间的中点
2 kkk
bax作为根的近似值，则 在二分过程中，可获得一个近似 根序列
}{ kx,该序列必以根 ξ 为极限由于,22 1?
k
kkk ababx?
所以，只要二分足够多次（即 k 充分大），便有
. kx
这里? 为事先给定的精度，再注意到 kx ≤
2
kk ab?,
所以，在实际计算时，只需某个有根区间的长度小于 ε,我们就可以停止计算，并取该有根区间的中点作为根 ξ 的近似值,
上述求根的方法被称为方程求根的二分法，简称为二分法,
下面的图框（如下页图）简单明了地描述了该方法,
各框的具体含义为,
框①,输入有根区间ba,及近似根的精度? ；
框②,开始二分；
框③,判断二分后新的有跟区间；
框④,检查是否满足精度；
框⑤,输出近似根 x,
开始输入
,,a b c
) (
0 a f
y
0
0
y y
N
x
b
b a
输出 x
N
①
②
③
④
⑤
结束
) (,
2
x f y x?
b a?
x
a
Y
Y
例 3 求 052)( 3 xxxf 在区间 ]3,2[ 之间的根,
解 016)1()3()2( ff,所以（ 2,3 ） 是 )( xf 的有根区间,
具体过程如下
x )( xf 有根区间
2 - 1
3 6
2.5 5.625 （ 2,2.5 ）
2.25 1.890625 （ 2.00,2.25 ）
2.125 0.345703 （ 2,2.125 ）
2.0625 - 0.351318 （ 2.0625,2.125 ）
2.09375 - 0.008942
2,1 0 9 3 7 5 0,1 6 6 8 3 6 （ 2,0 9 3 7 5,2,1 0 9 3 7 5 ）
（ 2.09375,2.125）
所以,1015625.2?x 即为所求方程的近似解，其误差限为
0078125.02 23 7
2,牛顿迭代法 已知非曲线方程 0)(?xf 的根的近似值为
nx,我们用曲线
)( xfy? 在点 nx 的切线（如下图）
)()()()( nnn xfxxxfxY
y
x
n x
1? n x
O
x n + 2 作为曲线 )( xf 的近似表达式,求方程
0)(?xY 的根，将其记为 1?nx,得
)(
)(
1
n
n
nn xf
xfxx
（ *）
定理 若 )( xf 满足 （ 1 ） 在ba,上存在二阶导数,且 )( xf 在该区间上不变号 ；
（ 2 ） )( xf? 在区间ba,上不等于零 ；
（ 3 ） 0)()( bfaf ； （ 4 ）
],[0 bax?,且有 0)()( 00 xfxf,则由递推公式（ * ）而得的序列nx 收敛于 0)(?xf 在区间ba,内的 惟 一根 ξ,
用式（ * ）中
1?n
x 作为方程 0)(?xf 的根的近似值,只要给定方程 0)(?xf 的根的初始近似值
0
x,由公式（ * ）反复迭代就可得到一个近似根序列 }{
n
x,可以证明：只要该近似根序列收敛就可断言其极限值必是 0)(?xf 的根,但是，
由于有时初始值选择不当，由式（ * ）所得的迭代序列不收敛，为此,我们给出下面的定理,
用上述公式（ * ）求 0)(?xf 的近似根的方法被称为牛顿迭代法，也叫切线法,下图是牛顿迭代法求根的框图,
1 0
x x?
,1 k k
开始
N
结束
N x,,
0
输入
k=1
0 ) ( f?
x
0 1 x x
N k?
输出迭代失败标志
Y
输出
1
x
Y
) ( f ) f( 0 0 0 1 x x x x
N
N
输出奇异标志
Y
在上面框图中，用 || 01 xx 来控制精度是合理的,因为，可以证明：若用后一次迭代值 1x 作为根? 的近似值，
则有 ||
1x ≤ || 01 xxC?,其中 C 是只与 )( xf 有关的确定常数，
0x 是前一次迭代值，所以当 || 01 xx? 足够小时，就可保证 || 1x 足够小,例 4 用牛顿发法求 1 1 5 的近似值，要求绝对误差限为 610?,
解 为求 115?x,需解方程,0115
2x
令 115)( 2 xxf,由牛顿迭代公式
，)( )(1
n
n
nn xf
xfxx
取初值 100?x,按式（ * * ）迭代 3 次便得精度为 610? 的结果,见下表
k kx
0 10
1 10.75
2 10.723 837
3 10.723 805
4 10.723 805
,2 115
2
1
n
n
nn x
xxx
整理得,
2
115
1
n
n
n
x
x
x
得
（ **）
思考题
2,方程求根的程序中，有几种精度控制方式？
1,只要 )( xf 在 ],[ ba 上具有二阶导数，由方程 0)(?xf 的牛顿迭代公式所得近似根序列 }{ kx 均收敛吗？
一,多项式的存在惟一性二,插值多项式的构造三,插值多项式的余项第二节 拉格朗日插值公式第二节 拉格朗日插值公式 当函数的解析表达式不知道时，往往通过实验测得它在一系列点（这些点称为节点）上的函数值,只知道表格函数
ix nxxxx?210
)( ii xfy? nyyyy?210
而对于
i
x 以外的点（即非节点）处的函数 )( xf 的值一无所知，为了获得函数在非节点处的函数值，就需要我们设法构造一个尽可能简单的函数（如多项式函数） )( xg 来近似代替函数 )( xf,为了得到近似函数 )( xg,我们自然要求它在节点 ix 处的函数值与 )( xf 在该节点处的函数值 iy 相同，更确切地说，我们要研究下列问题,
已知函数 )( xfy? 在 n +1 个相异节点 ),,1,0( nix i 处的函数值 ),,1,0()( niyxf ii,求一个次数不高于 n 的多项式
)( xP n,并使其满足插值条件
),,2,1,0()()( nixfxP iin （ 1 ）
点 ix 被称为插值节点,iy 称为插值节点处的函数值,)( xfy? 称为被插值函数,)( xP n 称为插值多项式,这样的问题称为插值问题,
)( xf 的满足插值条件（ 1 ）的多项式 )( xP
n
存在吗？若存在是否惟一呢？可以证明满足插值条件（ 1 ）的 )( xf 的插值多项式就是惟一的,事实上，若 )( xp,)( xq 都是 )( xf 的满足插值条件（ 1 ）的次数不超过 n 的多项式，且二者不相等,
则多项式 )()( xqxp? 也必是次数不超过 n 的多项式，注意到一,多项式的存在惟一性
0)()( iiii yyxqxp ).,,2,1,0( ni
便可知 )()( xqxp? 具有 n +1 个互异的零点，这与它是次数不超过 n 的多项式相矛盾，所以二者必相等,因此，惟一性得证,具体到插值多项式的存在性，我们只要能够设法将其构造出来就足以证明,
为了构造 )( xfy? 的满足插值条件（ 1 ）的次数不超过 n 的插值多项式 )( xP n,我们设
).()()()( 1100 xlyxlyxlyxP nnn
为使 ),,2,1,0()( niyxP iin,我们要求
.,0
,1)(
ji
jixl
ji （ 2）
二、插值多项式的构造并称 )( xl i 为插值基函数,它们均是次数不超过 n 的多项式,
显然，插值基函数 )( xl i 有 n 个零点：,,10 xx,..,
11, ii xx,,..,nx,因此,)( xl i 可写成如下形式,
).())(())(()( 1110 niii xxxxxxxxxxAxl
其中 A 为待定常数，又因为 1)(?ii xl 所以
).())(())((1 1110 niiiiiii xxxxxxxxxxA
于是
0 1 1 1
1
.
( ) ( ) ( ) ( ) ( )i i i i i i i n
A
x x x x x x x x x x
因此
0 1 1 1
0 1 1 1
( ) ( ) ( ) ( ) ( )
( ),
( ) ( ) ( ) ( ) ( )
i i n
i
i i i i i i i n
x x x x x x x x x x
lx
x x x x x x x x x x
所以
00 0
()
( ) ( )
()
nnn
j
n i i i
ii i ij
ij
xx
P x y l x y
xx
（ 3）
式（ 3 ）也叫做 n 次拉格朗日插值公式,其计算框图见下下页,
特别地，当 1?n 时，拉格朗日插值公式为一次多项式，称其为线性插值公式，即
,)(
01
0
1
10
1
01 xx
xxy
xx
xxyxP
当 n =2 时，有
))((
))((
))((
))((
)(
2101
20
1
2010
21
02 xxxx
xxxx
y
xxxx
xxxx
yxP
))((
))((
1202
10
2
xxxx
xxxx
y
称其为二次插值或抛物插值公式,
例 1 取节点 1,0 10 xx 和,00?x 211?x,12?x,对函数 xy e 分别建立线性和抛物插值公式,
解 （ 1 ）将 1,0 10 xx 及,1e 00y 11 ey 代入线性插值公式得 xy e 的线性插值公式为
1 ( ) 1 0,6 3 2 1 2 0 6,P x x （ 2 ）将 1,2
1,0
210 xxx 及,10?y
2
1
1 e
y,12 ey 代入抛物插值公式得,3 0 9 6 3 6 2.09 4 1 7 5 6 8.01)( 22 xxxP
利用 x?e 的线性插值和抛物插值公式分别计算 x?e 在点 2.0?x
处的近似值，得
,8 7 3 5 8.0)2.0(1?P 8 2 4 0 3.0)2.0(2?P,
显见,)2.0(
2
P 要比 )2.0(
1
P 更接近于
x?
e 在 0.2 处的值
0,8 1 8 73,一般地说，二次插值要比一次插值来得更精确些，但是要特别注意：并不是插值多项式的次数愈高愈精确，反而随着插值多项式次数的增高，计算量的增大，舍入误差的影响就会增大,
实践证明，在大范围内使用高次插值，逼 近 效果往往是 不 理想的，因此，实践中常用的是分段线性插值或分段抛物插值,
开始输入
1) n,1,2,(i,, i i xy x
y 1 k 0,
1 j 1,T
k? j?
) ( ) ( * j
k j
x x x x T T
1? n j 1
j j
y * T y y
1? n k
输出 y
Y
Y
N
N
Y
1 k k
构造插值多项式 )( xP n 的目的在于 为 )( xf 提供一个近似函数,
因而我们自然要关心用 )( xP n 代替 )( xf 所产生的误差是多少? 对此,
我们有下面的定理,
定理 设区间 ],[ ba 含有节点 nxxx,,,10?,而 )( xf 在开区间
),( ba 内有直到 1?n 阶导数,且
.),,2,1,0()( niyxf ii
n
k
k
n
n xxn
fxPxf
0
1
)()!1( )()()(?
已给定,则当 ],[ bax? 时,对于由式 (1) 给定的 )( xP n,有
(4),式中? 与节点 ix 及 x 有关,且 ),( ba,并称式 ( 4 ) 为插值多项式 )( xP n 的余项公式，记为 )( xR n,
三．插值多项式的余项证 ⑴当 x 是插值节点时，公式（ 4 ）显然成立,
⑵当 x 不是插值节点时，由于
)()()( xPxfxR nn
含有 1?n 个零点 nxxx,,,10?,于是
)())()(()( 10 nn xxxxxxxkxR （ 5 ）
其中 )( xk 为待定 函数，现在把 x 看成 ],[ ba 上一个固定点，作函数
))()(()()()( 10 xtxtxktPtft n )( nxt? 则 )( t? 有 2?n 个零点,x,nxxx,,,10?,即
0)()()()( 10 nxxxx
由拉格朗日中值定理知，若函数 () t? 在闭区间 ],[ ba 上连续，在开区间 ),( ba 内可导，且在区间端点处函数值相等，则至少有一点 ),( ba,使 0)( 成立,该结论也称为罗尔
（ R olle ）中值定理,
据上述罗尔中值定理,)( t 在 )( t? 的任意两个相邻零点间至少有一个零点，所以，在 )( t 内至少有 1?n 个零点，依次类推，
( 1 )
()
n
t?
在 ),( ba 内至少有一个零点，记为 ),( ba,使得
.0)()!1()()( )1()1( xknf nn
于是
) !,1(
)(
)(
)1(
n
f
xk
n?
),( ba,且依赖于 x,将 )( xk 代入式（ 5 ）便得要证明的公式（ 4 ）,
思考题在本节定理的证明过程中，下面两式为何成立,
1,0))(( )1(nn tP
2,)!1())())((( )1(10 nxtxtxt nn?
一,直线拟和二,将非多项式曲线拟合转化为线形拟合第三节 曲线拟和的最小二乘法第三节 曲线拟和的最小二乘法在许多实际问题中，自变量 x 和因变量 y 之间的函数关系往往不知道,但我们能通过实验测得它们在若干个点上的对应值 (,) ( 1,2,,)
ii
x y i n?,由这些数值，可以按拉格朗日插值方法构造出一个次数不超过 1n? 的插值多项式
1
()
n
Px
来逼近 ()y f x?,但是，由于测量工具及实验设备限制，所测得的函 数值 ( 1,2,
i
yi,) n 一般总会有误差,这就是说，所给的数据点 ),,2,1)(,( niyx
ii
未必都在曲线
)( xfy?
上，但由于大多数点基本上反应了 y 与 x 之间的对应规律，因此，所求的近似曲线
)( xy
未必要求
ii
yx?)(?,但要求所给数据点
),,2,1)(,( niyx
ii
尽可能在其附近,这样，我们 就有理由认为这条曲线是
)( xfy?
的逼近曲线，也称其拟和曲线,
为求数据点 ),,2,1)(,( niyx
ii
的拟和曲线，一般先在平面直角作标系 中描出所给的 数据点
),,2,1)(,( niyx
ii
图形，称这一图形为散点图,其次,观察散点图上的大多数点近似地在哪一条曲线附近,并凭经验确定出这条曲线的一般形式,最后，根据“让所有点均离其最近”这一原则具体的求出这条曲线,
若所给得数据点 ),,2,1)(,( niyx ii 近似在某直线附近，
可设拟合曲线的一般形式为 Y a bx,其中 ba,为待定系数,
为使所 求 直线 方程 满 足,使 所 给 数 据 点
,2,1)(,(?iyx ii ),n 尽量靠近一直线,自然要求 ba,使总误差一、直线拟和于是得
n
i
iii
n
i
ii
xbxay
bxay
1
1
,0)())((2
,0)1())((2
可以证明方程组（ 2 ）的解 a,b 就是使式（ 1 ）达到最小的解，称式（ 2 ）为正规方程组,上述方法通常称为最小二乘法,
n
i
ibxa
1
2))(
i(yQ
达到最小，因此必有
.0,0 ba QQ
(1)
即
11
2
1 1 1
( ),
( ) ( ),
nn
ii
ii
n n n
i i i i
i i i
n a x b y
x a x b x y
(2）
例 经实验测得某物理量 x,y 的 5 对数据如下，用最小二乘法求其拟合曲线,
解之得 9392.60a,5183.1?b,
5 702 758,
702 998 64 108 396,
ab
ab
所以 xy 5 1 8 3.19 3 9.60 为所求的拟合曲线,
解 首先描出散点图，发现其近似地在一直线上，
假设拟合曲线为 bxay,
将所给数据点代入正规方程组（ 2 ）得
141
148
123
125
150
172
123
126
165
187iy
ix
⑴ 若某组数据点 ),,2,1)(,( niyx ii 的散点图近似于一条双曲线，我们就可以设拟合曲线的一般形式为
.1 xbay
如果我们令 z
y
1,t
x
1
,则有 z a b t,此时 ba,可按线 性 最小二乘问题求出；
⑵ 若拟合曲线的一般形式为 ebxya?
将其两边取对数得
(a,b 为待定系数 ).
,lnln bxay 令 Aazy ln,ln,得 bxAz
则 A,b 可按线 性 最小二乘问题求出,
二、将非多项式曲线拟合转化为线性拟合思考题 ⒈根据“让所有点均离其最近”这一原则，数学上是如何刻画的,
⒉用线性最小二乘法求指数拟合 e bxya? 曲线的一般步骤是什么?
一,求积公式的建立二,求积公式的误差估计三,复化求积公式四,变步长的求积公式第四节 数 值 积 分第四节 数 值 积 分用牛顿 - 莱布尼茨公式 )()(d)( aFbFxxf
b
a
来计算定积分
b
a
xxfI d)( 的值的前提是,)( xf 的原函数 )( xF 能够求出,当 )( xf 由表格形式给出，或其原函数不易求出时，用牛顿 -- 莱布尼茨公式就解决不了问题,在此，我们希望用一个易于求原函数的函数（如多项式函数）来近似代替被积函数，然后计算这个多项式函数的积分，就会得到定积分的近似计算公式,
在［ a,b ］区间上选取 n+ １ 个等距节点,khax
k
,
k =0,1,2,?,n,其中
n
ab
h
为步长，并记
kk
yxf?)(,
下面，我们通过数据点组 ),(
kk
yx （ k =0,1,2,? n ）
来构造被积函数 )( xf 的 n 次插值多项式 )( xP
n
，即 一、求积公式的建立
.)()(
0
n
k
kkn xlyxP
其中 )( xl k 为插值基函数，所以
,d)(d)()(
0
n
k
b
a kk
b
a n
b
a xxlyxxPdxxf
称式（ 1 ）为牛顿 - 科茨（ Newton - Cotes ）求积公式,)( nkC
为牛顿 - 科茨系数，为计算 )( nkC,作变量替换
x a t h,于是
ba knk xxlabC d)(1)(
若记?
b
a k
n
k xxlabC d)(
1)(,则有求积公式
n
k
k
n
k
b
a
yCabxxf
0
)()(d)(
(1)
b
a
n
kj
j jk
j x
xx
xx
ab
d1
0
th
jk
jt
nh
n n
kj
j
d1
0
0
即 ()
0
0
( 1 )
[ ( ) ] d
! ( ) !
nk nn
n
k
j
jk
C t j t
n k n k
（ 2 ）
按式（ 2 ）计算科茨系数只涉及多项式函数的积分，
故不会有困难,将科茨系数列于下面,
288
19
288
75
288
50
288
50
288
75
288
19
5
90
7
90
32
90
12
90
32
90
7
4
8
1
8
3
8
3
8
1
3
6
1
6
4
6
1
2
2
1
2
1
1
)( n
k
Cn
当 n =1 时，牛顿－科茨求积公式被称为梯形公式，即
1 ( ) ( ) ( ),
2
T b a f a f b
当 n =2 时，牛顿－科茨公式称为抛物公式（或辛普森公式），即
)(
2
4)()(
6
1
bf
ba
fafabS,当 n =4 时，牛顿－科茨公式简称为科茨公式，即
)(7)(32)(12)(32)(7
90 43210
xfxfxfxfxf
ab
C
,
其中 4,3,2,1,0,4
1)( kabkax
k,
例 1 试分别用梯形公式、抛物公式、科茨公式计算定积分
1
0,5
dxx?,
解 ⑴梯形公式,
1
0,5
1 0,5d ( 0,5 1,0 ) 0,4 2 6 7 7 6 7
2
xx,
⑵抛物公式,
)0.175.045.0(
6
5.01
d
1
5.0
xx
4 3 0 9 3 4 0 3.0?,
(3) 科茨公式，此时步长,1 2 5.04 5.01h
.43 09 64 07.0
)0.17875.032
75.012625.0325.07(
90
5.01
d
1
5.0
xx
按牛顿 — 莱 布 尼茨公式计算，易知 1
5.0
2
31
5.0 3
2d xxx
.4 3 0 9 6 4 1 4.0?
若 )( xP n 为函数 )( xf 的 n 次插值多形项式，则
),()()( xPxfxR nn两边积分得
.d)(d)(d)( ba nbaba n xxPxxfxxR
所以，有,d)( b
a nn xxRIIfR
其中 (,)ab,于是，可得如下余项公式,
(1) 梯形公式的余项
).,(,)(12 )( 3 baabfTIR T
二、求积公式的误差估计
(2) 辛普森公式的余项
),(),()
2
(
180
)4(4 bafababSIR
S?
,
(3) 科茨公式的余项
),(),()
4(9 4 5
)(2 )6(6 bafababCIR
C?
,
从拉格朗日插值多项式的余项公式可以看到,代替被积函数的插值多项式次数越高，对被积函数光滑要求越高,另外，积分区间越小，求积公式的截断误差越小,因此，我们经常先把积分区间分成若干小区间，在每个区间上采用低阶求积公式，如梯形公式或抛物公式，然后再把它加起来得到整个区间上的求积公式，这是复化求积公式的基本思想,
三、复化求积公式
1,复化梯形公式将 区 间 nba ],[ 等 分,步 长
n
ab
h
,分点
khax k ),,2,1,0( nk,在每个小区间 ],[ 1?kk xx 上应用梯形公式，即得复化梯形公式
b
a
n
k
x
x
k
k
xxfxxf
1
0
1 d)(d)(
)]()([2 1
1
0
kk
n
k
xfxfh
1
1
[ ( ) 2 ( ) ( ) ],2 n k
k
h f a f x f b?
于是，有复化梯形求积公式
1
1
[ ( ) 2 ( ) ( ) ]
2
n
nk
k
h
T f a f x f b
（ 3 ）,
2,复化辛普森公式若记子段 ],[
1?kk
xx 的中点为
22
1
h
khax
k
,则复化辛普森公式为
)]()(4)([
6
1
1
0
2
1
kkk
n
k
n
xfxfxf
h
S （ 4 ）
类似可得复化科茨公式，这里从略,
例 2 分别用复化梯形公式，复化辛普森公式计算积分
1
20
4
d,
1
Ix
x
并将计算结果与准确值比较,
解 （ 1 ）准确值为
1 102
0
4 d 4 a r c t a n 3,1 4 1 5 9
1
I x x
x
（ 2 ） 利 用 复 化 梯 形 公 式 计 算,取
1 2 5.0
8
1
8
01
h,节点
88
1
0
k
kx
k
（ k = 0,1,2,? 8 ）,被积函数
2
1
4
)(
x
xf
在节点处的函数值分别为
9 3 8 4 6.3)
8
1
(,4)0( ff
5 0 6 8 5.3)
8
3
(,7 6 4 7 1.3)
8
2
( ff
8 7 6 4 0.2)
8
5
(,2 0 0 0 0.3)
8
4
( ff
2 6 5 4 9.2)
8
7
(,5 6 0 0 0.2)
8
6
( ff
2)1(?f,
8
0,1 2 5 1 2 3 4 5 6 7[ ( 0) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) ( 1 ) ]2 8 8 8 8 8 8 8T f f f f f f f f f
.13899.3?
（ 3 ） 为了利用上面求出的被积函数在相应节点处的函数值，取 4?n,此时步长
4
01?
h
1
0,2 5
4
,利用复化求积公式（ 4 ）得
8684822878583814)1()0(625.04 fffffffffS
.14159.3?
注意，在计算 4S 时，基本节点 )4,3,2,1,0(
4
0 k
k
x
k
,
中 间 节 点 1
2
1 1 1
4 2 4 4 8
k
kk
x
( 0,1,2,3 )k?,如
8
3
,
8
1
2
1
2
1 10 xx 等,
比较 8T,4S 与 I 的准确值会发现，尽管 8T 与 4S 的计算量基本相同，但 4S 的精度要比 8T 的精度高，可见，
复化辛普森公式是一种精度较高的求积公式,
应用复化求积公式计算定积分，关键是步 长选择,
步长太小，计算量加大，舍入误差也随之增大；步长太大，截断误差增大，精度也难以保证,所以,利用复化求积公式计算时，通常采用变步长的求积方法,即 把步长逐次减半，在步长逐次减半的过程中，反复利用复化求积公式进行计算，直到达到所 需 的精度,
若记 nI 为步长为 h 时的积分值,nI 2 为将 h 分半后的积分值,? 为给定的精度要求，可以证明四 步长的求积公式当 2nnII 时，便有
2
( ) d
b
n
a
f x x I,这样，用二分前后的两个积分值 nI 与 2 nI 之差的绝对值 2nnII? 是否小于预先给定的精度?,来决定二分过程是否停止,
由于用数值积分公式计算定积分，主要工作量在函数值的计算上，因此在应用变步长的复化求积公式时，应考虑步长二分前后计算公式间的关系，以便充分利用已求出的信息（函数值）,
1.变步长梯形法则的递推公式对积分 xxf
b
a
d)(? 应用变步长的复合梯形求积公式,把区间 [ a,b ] 逐次二等分，设现已把 [ a,b ] n 等分，步长为
n
ab
h
,分点为
,,,,,,,,1210 bxxxxxxa nkk
则由复化梯形公式（ 3 ），知
.)()(2)(2
1
1
n
k
kn bfxfaf
hT
然后，将步长
n
ab
h
二分，即在每个子区间
],[ 1 kk xx? 中增加一个二分点 ),,2,1(
2
1 nkx k,便将
[ a,b ] 区间 2 n 等分，此时步长为
n
abh
22
,分点为
,,,,,,,,,,212121 122,110 bxxxxxxxxxa nkkk
则由复化梯形公式，知
1
1
1
1
2
)()(2)(2)(
4 2
1
n
k
n
k
kkn
bfxfxfaf
h
T
n
k
k
n
k
k xf
h
bfxfaf
h
1
1
1
)(
2
)()(2)(
4 2
1
n
k
kn
xf
h
T
1
)(
22
1
2
1,
即即
n
k
knn xf
hTT
1
2 )(22
1
2
1 （ 5 ）
为变步长的梯形递推公式,使用该公式时，要特别注意，
n
ab
h
为二分前的步长,?
n
k
k
xf
1
)(
2
1
为被积函数
)( xf 在二分后所有新增节点处函数值之和,下页图给出了变步长的梯形求积计算框图,
同样的道理，可讨论变步长的辛普森求积法则,
开始输入 a,b,c
a b h
2 )) ( ) ( ( 1 h b f a f T
s? 0,
x? a + h / 2
b? x
ε
T
T
1 2
输出
2 T
Y
N
N
Y
s? s + f ( x ),x = x + h
T 2
h s / 2
h? h/2,
T
1
例 3 用变步长方法计算 xx xI ds i n10,解 先对整个区间 [0,1] 用梯形公式,对于
x
xxf s i n)(?,1s i nl i m)0(
0
x
xf
x
,8 4 1 4 7 1 0.0)1(?f,故
1
1 [ ( 0 ) ( 1 ) ] 0,9 2 0 7 3 5 5,
2T f f
然后将区间 [0,1] 二分，由于 9 5 8 8 5 1 1.0)
2
1
(?f,利用递推公式（ 5 ）得
21
1 1 1
( ) 0,9 3 9 7 9 3 3,
2 2 2
T T f
再将每一个小区间二分一次，并计算新增分点上的函数值,
9 8 9 6 1 5 8.0)
4
1
(?f,9088517.0)
4
3
(?f,再利用公式 （ 5 ）,得
9 4 4 5 1 3 5.0)43()41(2 5.021 24 ffTT 这样不断二分下去，计算结果见 下 表 （表中 k 代表二分次数，区间等分数 kn 2? ）,积分准确值为
0,9 4 6 0 8 3 1,用变步长方法二分 10 次得到了这个结果,
0.94605965
0.9460831100.94598504
0.946083090.94569093
0.946082780.94451352
0.946081570.93979331
0.946076960.9207355 0
k 2kT k 2kT
思考题
1,为什么说抛物公式比梯形公式精确度高？
2,复化梯形公式与变步长的梯形递推公式在算法上有何区别与联系？
一,欧拉方法二,改进的欧拉方法三,龙格 — 库塔方法四,误差的分析第五节 常微分方程的数值解法应用解析方法求解常微分方程初值问题,?
00 )(
),(
d
d
yxy
yxf
x
y
( 1 )
只能求出一些特殊类型的方程 的 解,在大多数情况下，对初值问题（ 1 ），只 需 给出其数值解就能满足工程需要,初值问题 ( 1 ) 的数值解就是对未知函数 )( xfy? 给出在一系列节点 nxxx,,,10? 处的函数值 nyyyy,,,,210?,
第五节 常微分方程的数值解法将方程 ))(,( xyxfy 两边对 x 从 ix 到 1?ix 积分，得 1
1( ) ( ) (,( ) ) d
i
i
x
ii xy x y x f x y x x
( 2 )
用公式（ 2 ）计算函数值 )(
1?i
xy 的难点在于积分
1
d))(,(
i
i
x
x
xxyxf 的计算,下面我们就用不同的数值积分公式代替此积分，得出相应的微分方程数值解公式,
由于?
1
))(,(d))(,(
i
i
x
x
ii xyxhfxxyxf
其中 ii xxh 1 为步长，于是有一、欧拉方法用 iy 替代上式右端的 )( ixy,并将由式（ 3 ）右端算出的值 1?iy 作为 )( 1?ixy 的近似值，这样建立的计算格式
),,2,1,0(),(1 niyxhfyy iiii （ 4 ）
称为欧拉公式,用此公式求解初值问题（ 1 ）的数值解的方法叫做欧拉方法,欧拉方法的特点是算法简单，但精度不高,
))(,()()( 1 iiii xyxhfxyxy （ 3）
为了提高精度，改用梯形公式计算求积项
) ) ](,())(,([
2
d))(,(
11
1
iiii
x
x
xyxfxyxf
h
xxyxf
i
i
,
将上式代入（ 2 ）式得
))](,())(,([
2
)()(
111
iiiiii
xyxfxyxf
h
xyxy,( 5 )
将其中 )(),(
1 ii
xyxy
分别用
ii
yy,
1?
代替，则得到
)].,(),([2 111 iiiiii yxfyxfhyy
称其为梯形公式,该公式的右端含有未知的 1?iy,因此，
称其为隐式格式,为了应用公式（ 5 ）计算 1?iy,我们先用欧拉公式（ 4 ）计算 1?iy 的近似值 1?iy,称其为预二、改进的欧拉方法报值，然后用预报值 1?iy 代替公式（ 5 ）中右端的 1?iy,再计算出校正值 1?iy,这样建立起来的预报 — 校正系统称为改进的欧拉公式，记为
）（校正预报
niyxfyxf
h
yy
yxhfyy
iiiiii
iiii
,,2,1,0)()],(),([
2
)(),(
111
1
（ 6）
为了减少计算量，可将式（ 6 ）写为
)(
2
1
),,2,1,0(),(
),(
1
1
cpi
piic
iiip
yyy
niyxhfyy
yxhfyy
（ 7 ）
请画出改进的欧拉方法的计算框图,
例 1 取步长 h =0.1,求初值问题
d2
,
d
( 0 ) 1
yx
y
xy
y
在 x =0 到 x =1 各节点上的数值解,
解 所给的初值问题的改进的欧拉公式为
2
0.1 ( )
i
p i i
i
x
y y y
y
,
1
2
0,1 ( )
i
c i p
p
x
y y y
y
,
1
1
()
2
i p c
y y y
,
表 2
计算结果见下 表,
i
x 0 0,1 0,2 0,3 0,4 0,5
i
y 1 1,09 5 9 1,18 4 1 1,26 6 2 1,34 3 4 1,41 6 4
()
i
yx 1 1,09 5 4 4 5 1,18 3 2 1 6 1,26 4 9 1 1 1341641 1,41 4 2 1 4
i
x 0 0,6 0,7 0,8 0,9 1,0
i
y
1 1,4 860 1,55 2 5 1,61 6 5 1,67 8 2 1,73 7 9
()
i
yx 1 1,48 3 2 4 0 1,54 9 1 9 3 1,61 2 4 5 2 1,67 3,32 0 1,73 2 0 5 1
该初始值问题的解析解为 xy 21,按这个式子计算出 y 在相应节点处的值（下表第三行是相应的精确值），与准确解比较，会发现按改进的欧拉公式计算出的函数大致有三位有效数字,下面给出的龙格 — 库塔（ Runge — Kutta ）方法具有更高的精度,
对于初值问题（ 1 ），我们不加证明给出如下计算格式,
)22(
6
),(
)
2
,
2
(
)
2
,
2
(
),(
43211
34
23
12
1
kkkk
h
yy
hkyhxfk
k
h
y
h
xfk
k
h
y
h
xfk
yxfk
ii
ii
ii
ii
ii
公式（ 8 ）称为四阶龙格 — 库塔公式，其计算框图 见下页,
（ 8）
三、龙格 — 库塔方法开始结束
Y
输出输入
n h,,y,x
0 0
) 2 hk y,2 h f( x k
2
0 0
3
) 2 hk y
2,h f( x k 1
0 0
2
) hk y,f(x k 3 0
1
4
) y,f ( x k
0 0 1
定义
1 i y ),f ( x,?
x 1 = x
0 + h
6 ) k 2k 2k h(k y y 4 3 2 1 0 1
1 1 y,x
n? i?
N
,
1 0 x x
1 0
y y?
,1 i i
例 2 四阶龙格 — 库塔方法求解例 1中的初值问题,
解 初值问题
2,
( 0) 1
xyy
y
y
(0 ≤ x ≤ 1)
的四阶龙格 — 库塔公式为（步长 2.0?h ）
i
i
i y
xyk 2
1
21
1
0.2
2( )0.2
2,
0.22
2
i
i
i
x
k y k
yk
计算结果见下表直接验证可知,由 四阶龙格 — 库塔方法所得到的结果具有四位有效数字,
ix
iy
0 0.2 0.4 0.6 0.8 1.0
1 1.183229 1.341667 1.483281 1.612514 1.732142
,
2
2.0
)
2
2.0
(2
2
2.0
2
23
ky
x
kyk
i
i
i
3
34 2.0
)2.0(22.0
ky
xkyk
i
i
i?
)22(6 2.0 43211 kkkkyy ii
,
.
若取步长 h =0.05,用四阶龙格 — 库塔方法计算例
1 中 的 初 值 问 题 在 节 点 2.01?x,?2x
0.4,8.0,6.0 43 xx 处函数 )( xy 的值，有
,3 4 1 6 4 1.1,1 8 3 2 1 6.1 21 yy
.6 1 2 4 5 2.1,4 8 3 2 4 0.1 43 yy再将它们与准确值比较，会发现比步长取 0.2 时的解还要精确,
在实际计算中，并不是步长取得愈小愈好，因为，
步长变小，计算量增大，舍入误差也相应地增大，一般地，我们根据事先给定的精度要求，在计算过程中，
逐步缩小步长,具体作法如下,
四、误差的分析设给定的精度要求为 ε,从节点
i
x 出发，先以
h 为步长，利用四阶龙格 — 库塔方法得到 )(
1?i
xy 的近似值为
)(
1
h
i
y
，然后将步长 h 分半，即以
2
h
为步长，从
i
x 出发，经两步按四阶龙格 — 库塔方法计算出 )( 1?ixy 的近似值
)(
1
2
h
i
y
,可以证明,若不等式
)()( 12 hii yy h （ 9）
成立，则将
)(
1
2
h
iy? 作为 )( 1?ixy 的近似值即满足精度要求,否则，进一步将“步长折半”进行计算，直到不等式（ 9 ）成立后，再取最后一个
)(
1
2
h
iy? 作为
)( 1?ixy 的近似值即可,
以上计算过程中，步长自动选取，故称其为变步长的龙格 - 库塔方法,用不等式（ 9 ）也可以实现变步长的改进欧拉方法等,
下页图给出了变步长的四阶龙格 - 库塔方法计算框图，仿该计算框图可画出 变 步长的改进欧拉方法的计算框图,
1 1 1
t t + h?
x=x
0
+h
0
结束开始
b 输入
n,ε,y ),f(x,,y,x 0 0
输出
2
y x,
b? x?
Y
x,
x 0?
2 0
y y?
1 2
y
y
N
N
1
y
2
y?
x? t
1?
2 0
y g?
N
h
0
=(b － x
0
)/n
计算
y
1
1
1 0 1
h
h,g y
2
10
tx?
计算
y
2
Y
0 1 h h
<
Y
计算 y1
计算 y2
k1=f(x0,y0),
k2=f(x0+h0/2,y0+ h0 k1/2),
k3=f(x0+h0/2,y0+ h0 k2/2),
k4=f(x,y0+ h0k3),
y1= y0+( k1+2k2+2k3+ k4) h0/6.
k1=f(t1,g0),
k2=f(t1+h1/2,g0+ h1 k1/2),
k3=f(t1+h1/2,g0+ h1 k2/2),
k4=f(t1 +h1,g0+ h1k3),
y2= g0+( k1+2k2+2k3+ k4) h1/6.
思考题
⒈ 是否每个微分方程都能求其数值解,
⒉ 用欧拉方法能够求出初值问题的足够精确的解吗？
1.绝对误差与相对误差设 *x 为准确数 x 的近似值，我们称 xxxe **)( 为近似数 *x
的绝对误差，简称误差,由于一般无法得到准确值 x,因此绝对误差
)( *xe 也无法直接算出，如果能估计其绝对值的范围
.**)( xxxe则称? 为近似数 *x 的绝对误差限，简称误差限,
有了绝对误差限就可以知道准确值 x 的范围,
.xxxxx
一、误差即 x 落在闭区间 ],[
**
xx 上，在应用上，常采用如下写法来刻画
*
x 的精度
.
*
xx
绝对误差限并不能完全表示近似值的好坏程度，例如
110x
51 0 0 0y
虽然 x 的绝对误差限比 y 的绝对误差限小，但 1 0 0 0 作为 y
的近似值要比 10 作为 x 的近似值要好，为了清楚的描述这一现象，
我们引进
,)()( *** ** x xxxxexr
并称 )( *xr 为近似数 *x 的相对误差,
显然,)( *xr 的准确值也无法直接得到，如果我们知道
,)()( *
*
*
x
xexr
则称 r? 为近似数 *x 的相对误差限,
若? 为?x 的绝对误差限，则
*xr
前面提到的 110x 的近似值 10*?x 的相对误差限为
1 0 %,而 y = 1 0 0 0 ± 5 的近似值 1000*?y 的相对误差限是 0,5 %,
由此可见，相对误差限愈小，近似程度愈高,
2．有效数字定义 （有效数字） 若近似值 *x 的误差限是某一数位的半个单位，且该位到 *x 的左边第一位非零数字共有 n 位,则称
*x 有 n 位有效数字,
例 1 按四舍五入的原则写出下列各数具有 5 位有效数字的近似值,1 8 7,9 3 2 5,0,0 3 7 8 5 5 1,8,0 0 0 0 3 3,2,7 1 8 2 8 解 按定义，上述各数具有 5 位有效数字的近似值分别为,
187.93,0.037 855,8.000,2.718 3,按有效数字的定义，近似数的最后一个有效数位确实能反映其绝对误差的大小，并且可以证明：有效数字位数越多，相对误差限越小,
3,数值运算的误差估计设函数,,( 21 xxfu ),,nx 如果自变量,,21 xx? nx,
的某组值,,21 xx? nx,的近似值分别为,,*2*1 xx? *,nx,则函数
u 的近似值 fu?*,,( *2*1 xx? ),*nx,且 *u 的误差
),,,(),,,()( 21**2*1** nn xxxfxxxfuuue
).(d *
1
*
*
* kkxx
n
k k
xx
xxx xxx
fu
kk
kk
kkk
),()( *
1
*
* kxx
n
k k
xexfue
kk
即从而,*u 的误差限为 ).()( *
1
*
*
k
xx
n
k k
xxfu
kk
*u 的相对误差限为
.)()(
*1
*
*
*
*
kk xxk
n
k
k
r x
f
u
x
u
u
例 2 已测得某场地长 l 的近似值为 *l =110 m,宽
d 的近似值为 *d =80 m,已知
ll?* ≤ 0.2 m,dd?* ≤ 0.1 m,
试求面积 ldA? 的绝对误差与相对误差限,
解 ldA?
ldAdlA,
)()()( *** dlldA = 8 0 ×0,2 + 1 1 0 ×0,1
=27 )( 2m,031.0
80100
27)(
*
*
AAr
4,数值运算应注意的若干原则一个实际的计算问题，其数值运算次数以千万次计，我们不可能对其每步计算都进行误差分析，因此，在设计算法时要特别注意误差的积累,具体地说要注意如下几点,
⑴要尽量减少运算次数，以减少舍入误差的积累；
⑵要尽量避免两相近数相减，以防有效数字的严重损失而影响精度； ⑶要尽量避免使用“小分母”，以防运算结果过大而造成溢出；
⑷要防止大数“吃掉”小数现象,例如，在五位计算机上计算
.9.052492
1000
1
k
A
把运算写成规格化形式?
1000
1
55 10009000.01052492.0
k
A
555 100 0 0 0 0 9.0100 0 0 0 0 9.0105 2 4 9 2.0
51052492.0
该结果显然不正确，这是由于在计算机内计算时要对阶，每个 0,9 都被看成了零,要避免这一现象，就应先把数量级相同的数相加，然后在加上大数 52 4 92,也就避免了大数吃掉小数,
方程的根：方程 0)(?xf 的根也叫做函数 )( xf 的根，或叫做 )( xf 的零点,方程的根的存在性：若
)( xf 在闭区间ba,上连续，且
0)()( bfaf,则 )( xf 在开区间 ),( ba 内至少有一实根,
1,方程求根的二分法设 )( xf 为闭区间ba,上的连续函数，且在开区间
),( ba 内有惟一的单实根 ξ,下面，我们给出求单实根 ξ 的近似值方法,
二、方程求根取区间ba,的中点
2
0
ba
x
,考察区间 ],[
0
xa 及 ],[
0
bx
中哪一个为有根区间，即检查 )(
0
xf 与 )( af 是否同号,如果确系同号，说明所求的根 ξ 在区间 ],[ 0 bx 中，即新的有根区间为 ],[ 0 bx,这时令 bbxa 101,；否则 ξ 必在 ],[
0
xa 中，这时令 aa?1,01 xb?,即不管出现哪种情形，新的有根区间
],[
11
ba 的长度仅为原有根区间的一半,
如此反复二分下去，即可得出一系列有根区间 ],[ ba,
],[ 11 ba,],[ 22 ba,,? ],[ kk ba,
其中每一个区间都是前一个区间的一半，因此，二分 k
次后的有根区间 ],[ kk ba 的长度为对于新的有根区间 ],[ 11 ba 又可施行同样的作法，即用中点
2
11
1
ba
x
将区间 ],[ 11 ba 再分为两半，然后判定所求的根 ξ 在 1x
的哪一侧？从而确定出新的有根区间 ],[ 22 ba,其长度为区间
],[ 11 ba 的一半,
.2 kkk abab
由此可见,如果二分过程无限的继续下去，这些有根区间必收缩于点 ξ,该点就是所求的根,每次二分后，若取有根区间的中点
2 kkk
bax作为根的近似值，则 在二分过程中，可获得一个近似 根序列
}{ kx,该序列必以根 ξ 为极限由于,22 1?
k
kkk ababx?
所以，只要二分足够多次（即 k 充分大），便有
. kx
这里? 为事先给定的精度，再注意到 kx ≤
2
kk ab?,
所以，在实际计算时，只需某个有根区间的长度小于 ε,我们就可以停止计算，并取该有根区间的中点作为根 ξ 的近似值,
上述求根的方法被称为方程求根的二分法，简称为二分法,
下面的图框（如下页图）简单明了地描述了该方法,
各框的具体含义为,
框①,输入有根区间ba,及近似根的精度? ；
框②,开始二分；
框③,判断二分后新的有跟区间；
框④,检查是否满足精度；
框⑤,输出近似根 x,
开始输入
,,a b c
) (
0 a f
y
0
0
y y
N
x
b
b a
输出 x
N
①
②
③
④
⑤
结束
) (,
2
x f y x?
b a?
x
a
Y
Y
例 3 求 052)( 3 xxxf 在区间 ]3,2[ 之间的根,
解 016)1()3()2( ff,所以（ 2,3 ） 是 )( xf 的有根区间,
具体过程如下
x )( xf 有根区间
2 - 1
3 6
2.5 5.625 （ 2,2.5 ）
2.25 1.890625 （ 2.00,2.25 ）
2.125 0.345703 （ 2,2.125 ）
2.0625 - 0.351318 （ 2.0625,2.125 ）
2.09375 - 0.008942
2,1 0 9 3 7 5 0,1 6 6 8 3 6 （ 2,0 9 3 7 5,2,1 0 9 3 7 5 ）
（ 2.09375,2.125）
所以,1015625.2?x 即为所求方程的近似解，其误差限为
0078125.02 23 7
2,牛顿迭代法 已知非曲线方程 0)(?xf 的根的近似值为
nx,我们用曲线
)( xfy? 在点 nx 的切线（如下图）
)()()()( nnn xfxxxfxY
y
x
n x
1? n x
O
x n + 2 作为曲线 )( xf 的近似表达式,求方程
0)(?xY 的根，将其记为 1?nx,得
)(
)(
1
n
n
nn xf
xfxx
（ *）
定理 若 )( xf 满足 （ 1 ） 在ba,上存在二阶导数,且 )( xf 在该区间上不变号 ；
（ 2 ） )( xf? 在区间ba,上不等于零 ；
（ 3 ） 0)()( bfaf ； （ 4 ）
],[0 bax?,且有 0)()( 00 xfxf,则由递推公式（ * ）而得的序列nx 收敛于 0)(?xf 在区间ba,内的 惟 一根 ξ,
用式（ * ）中
1?n
x 作为方程 0)(?xf 的根的近似值,只要给定方程 0)(?xf 的根的初始近似值
0
x,由公式（ * ）反复迭代就可得到一个近似根序列 }{
n
x,可以证明：只要该近似根序列收敛就可断言其极限值必是 0)(?xf 的根,但是，
由于有时初始值选择不当，由式（ * ）所得的迭代序列不收敛，为此,我们给出下面的定理,
用上述公式（ * ）求 0)(?xf 的近似根的方法被称为牛顿迭代法，也叫切线法,下图是牛顿迭代法求根的框图,
1 0
x x?
,1 k k
开始
N
结束
N x,,
0
输入
k=1
0 ) ( f?
x
0 1 x x
N k?
输出迭代失败标志
Y
输出
1
x
Y
) ( f ) f( 0 0 0 1 x x x x
N
N
输出奇异标志
Y
在上面框图中，用 || 01 xx 来控制精度是合理的,因为，可以证明：若用后一次迭代值 1x 作为根? 的近似值，
则有 ||
1x ≤ || 01 xxC?,其中 C 是只与 )( xf 有关的确定常数，
0x 是前一次迭代值，所以当 || 01 xx? 足够小时，就可保证 || 1x 足够小,例 4 用牛顿发法求 1 1 5 的近似值，要求绝对误差限为 610?,
解 为求 115?x,需解方程,0115
2x
令 115)( 2 xxf,由牛顿迭代公式
，)( )(1
n
n
nn xf
xfxx
取初值 100?x,按式（ * * ）迭代 3 次便得精度为 610? 的结果,见下表
k kx
0 10
1 10.75
2 10.723 837
3 10.723 805
4 10.723 805
,2 115
2
1
n
n
nn x
xxx
整理得,
2
115
1
n
n
n
x
x
x
得
（ **）
思考题
2,方程求根的程序中，有几种精度控制方式？
1,只要 )( xf 在 ],[ ba 上具有二阶导数，由方程 0)(?xf 的牛顿迭代公式所得近似根序列 }{ kx 均收敛吗？
一,多项式的存在惟一性二,插值多项式的构造三,插值多项式的余项第二节 拉格朗日插值公式第二节 拉格朗日插值公式 当函数的解析表达式不知道时，往往通过实验测得它在一系列点（这些点称为节点）上的函数值,只知道表格函数
ix nxxxx?210
)( ii xfy? nyyyy?210
而对于
i
x 以外的点（即非节点）处的函数 )( xf 的值一无所知，为了获得函数在非节点处的函数值，就需要我们设法构造一个尽可能简单的函数（如多项式函数） )( xg 来近似代替函数 )( xf,为了得到近似函数 )( xg,我们自然要求它在节点 ix 处的函数值与 )( xf 在该节点处的函数值 iy 相同，更确切地说，我们要研究下列问题,
已知函数 )( xfy? 在 n +1 个相异节点 ),,1,0( nix i 处的函数值 ),,1,0()( niyxf ii,求一个次数不高于 n 的多项式
)( xP n,并使其满足插值条件
),,2,1,0()()( nixfxP iin （ 1 ）
点 ix 被称为插值节点,iy 称为插值节点处的函数值,)( xfy? 称为被插值函数,)( xP n 称为插值多项式,这样的问题称为插值问题,
)( xf 的满足插值条件（ 1 ）的多项式 )( xP
n
存在吗？若存在是否惟一呢？可以证明满足插值条件（ 1 ）的 )( xf 的插值多项式就是惟一的,事实上，若 )( xp,)( xq 都是 )( xf 的满足插值条件（ 1 ）的次数不超过 n 的多项式，且二者不相等,
则多项式 )()( xqxp? 也必是次数不超过 n 的多项式，注意到一,多项式的存在惟一性
0)()( iiii yyxqxp ).,,2,1,0( ni
便可知 )()( xqxp? 具有 n +1 个互异的零点，这与它是次数不超过 n 的多项式相矛盾，所以二者必相等,因此，惟一性得证,具体到插值多项式的存在性，我们只要能够设法将其构造出来就足以证明,
为了构造 )( xfy? 的满足插值条件（ 1 ）的次数不超过 n 的插值多项式 )( xP n,我们设
).()()()( 1100 xlyxlyxlyxP nnn
为使 ),,2,1,0()( niyxP iin,我们要求
.,0
,1)(
ji
jixl
ji （ 2）
二、插值多项式的构造并称 )( xl i 为插值基函数,它们均是次数不超过 n 的多项式,
显然，插值基函数 )( xl i 有 n 个零点：,,10 xx,..,
11, ii xx,,..,nx,因此,)( xl i 可写成如下形式,
).())(())(()( 1110 niii xxxxxxxxxxAxl
其中 A 为待定常数，又因为 1)(?ii xl 所以
).())(())((1 1110 niiiiiii xxxxxxxxxxA
于是
0 1 1 1
1
.
( ) ( ) ( ) ( ) ( )i i i i i i i n
A
x x x x x x x x x x
因此
0 1 1 1
0 1 1 1
( ) ( ) ( ) ( ) ( )
( ),
( ) ( ) ( ) ( ) ( )
i i n
i
i i i i i i i n
x x x x x x x x x x
lx
x x x x x x x x x x
所以
00 0
()
( ) ( )
()
nnn
j
n i i i
ii i ij
ij
xx
P x y l x y
xx
（ 3）
式（ 3 ）也叫做 n 次拉格朗日插值公式,其计算框图见下下页,
特别地，当 1?n 时，拉格朗日插值公式为一次多项式，称其为线性插值公式，即
,)(
01
0
1
10
1
01 xx
xxy
xx
xxyxP
当 n =2 时，有
))((
))((
))((
))((
)(
2101
20
1
2010
21
02 xxxx
xxxx
y
xxxx
xxxx
yxP
))((
))((
1202
10
2
xxxx
xxxx
y
称其为二次插值或抛物插值公式,
例 1 取节点 1,0 10 xx 和,00?x 211?x,12?x,对函数 xy e 分别建立线性和抛物插值公式,
解 （ 1 ）将 1,0 10 xx 及,1e 00y 11 ey 代入线性插值公式得 xy e 的线性插值公式为
1 ( ) 1 0,6 3 2 1 2 0 6,P x x （ 2 ）将 1,2
1,0
210 xxx 及,10?y
2
1
1 e
y,12 ey 代入抛物插值公式得,3 0 9 6 3 6 2.09 4 1 7 5 6 8.01)( 22 xxxP
利用 x?e 的线性插值和抛物插值公式分别计算 x?e 在点 2.0?x
处的近似值，得
,8 7 3 5 8.0)2.0(1?P 8 2 4 0 3.0)2.0(2?P,
显见,)2.0(
2
P 要比 )2.0(
1
P 更接近于
x?
e 在 0.2 处的值
0,8 1 8 73,一般地说，二次插值要比一次插值来得更精确些，但是要特别注意：并不是插值多项式的次数愈高愈精确，反而随着插值多项式次数的增高，计算量的增大，舍入误差的影响就会增大,
实践证明，在大范围内使用高次插值，逼 近 效果往往是 不 理想的，因此，实践中常用的是分段线性插值或分段抛物插值,
开始输入
1) n,1,2,(i,, i i xy x
y 1 k 0,
1 j 1,T
k? j?
) ( ) ( * j
k j
x x x x T T
1? n j 1
j j
y * T y y
1? n k
输出 y
Y
Y
N
N
Y
1 k k
构造插值多项式 )( xP n 的目的在于 为 )( xf 提供一个近似函数,
因而我们自然要关心用 )( xP n 代替 )( xf 所产生的误差是多少? 对此,
我们有下面的定理,
定理 设区间 ],[ ba 含有节点 nxxx,,,10?,而 )( xf 在开区间
),( ba 内有直到 1?n 阶导数,且
.),,2,1,0()( niyxf ii
n
k
k
n
n xxn
fxPxf
0
1
)()!1( )()()(?
已给定,则当 ],[ bax? 时,对于由式 (1) 给定的 )( xP n,有
(4),式中? 与节点 ix 及 x 有关,且 ),( ba,并称式 ( 4 ) 为插值多项式 )( xP n 的余项公式，记为 )( xR n,
三．插值多项式的余项证 ⑴当 x 是插值节点时，公式（ 4 ）显然成立,
⑵当 x 不是插值节点时，由于
)()()( xPxfxR nn
含有 1?n 个零点 nxxx,,,10?,于是
)())()(()( 10 nn xxxxxxxkxR （ 5 ）
其中 )( xk 为待定 函数，现在把 x 看成 ],[ ba 上一个固定点，作函数
))()(()()()( 10 xtxtxktPtft n )( nxt? 则 )( t? 有 2?n 个零点,x,nxxx,,,10?,即
0)()()()( 10 nxxxx
由拉格朗日中值定理知，若函数 () t? 在闭区间 ],[ ba 上连续，在开区间 ),( ba 内可导，且在区间端点处函数值相等，则至少有一点 ),( ba,使 0)( 成立,该结论也称为罗尔
（ R olle ）中值定理,
据上述罗尔中值定理,)( t 在 )( t? 的任意两个相邻零点间至少有一个零点，所以，在 )( t 内至少有 1?n 个零点，依次类推，
( 1 )
()
n
t?
在 ),( ba 内至少有一个零点，记为 ),( ba,使得
.0)()!1()()( )1()1( xknf nn
于是
) !,1(
)(
)(
)1(
n
f
xk
n?
),( ba,且依赖于 x,将 )( xk 代入式（ 5 ）便得要证明的公式（ 4 ）,
思考题在本节定理的证明过程中，下面两式为何成立,
1,0))(( )1(nn tP
2,)!1())())((( )1(10 nxtxtxt nn?
一,直线拟和二,将非多项式曲线拟合转化为线形拟合第三节 曲线拟和的最小二乘法第三节 曲线拟和的最小二乘法在许多实际问题中，自变量 x 和因变量 y 之间的函数关系往往不知道,但我们能通过实验测得它们在若干个点上的对应值 (,) ( 1,2,,)
ii
x y i n?,由这些数值，可以按拉格朗日插值方法构造出一个次数不超过 1n? 的插值多项式
1
()
n
Px
来逼近 ()y f x?,但是，由于测量工具及实验设备限制，所测得的函 数值 ( 1,2,
i
yi,) n 一般总会有误差,这就是说，所给的数据点 ),,2,1)(,( niyx
ii
未必都在曲线
)( xfy?
上，但由于大多数点基本上反应了 y 与 x 之间的对应规律，因此，所求的近似曲线
)( xy
未必要求
ii
yx?)(?,但要求所给数据点
),,2,1)(,( niyx
ii
尽可能在其附近,这样，我们 就有理由认为这条曲线是
)( xfy?
的逼近曲线，也称其拟和曲线,
为求数据点 ),,2,1)(,( niyx
ii
的拟和曲线，一般先在平面直角作标系 中描出所给的 数据点
),,2,1)(,( niyx
ii
图形，称这一图形为散点图,其次,观察散点图上的大多数点近似地在哪一条曲线附近,并凭经验确定出这条曲线的一般形式,最后，根据“让所有点均离其最近”这一原则具体的求出这条曲线,
若所给得数据点 ),,2,1)(,( niyx ii 近似在某直线附近，
可设拟合曲线的一般形式为 Y a bx,其中 ba,为待定系数,
为使所 求 直线 方程 满 足,使 所 给 数 据 点
,2,1)(,(?iyx ii ),n 尽量靠近一直线,自然要求 ba,使总误差一、直线拟和于是得
n
i
iii
n
i
ii
xbxay
bxay
1
1
,0)())((2
,0)1())((2
可以证明方程组（ 2 ）的解 a,b 就是使式（ 1 ）达到最小的解，称式（ 2 ）为正规方程组,上述方法通常称为最小二乘法,
n
i
ibxa
1
2))(
i(yQ
达到最小，因此必有
.0,0 ba QQ
(1)
即
11
2
1 1 1
( ),
( ) ( ),
nn
ii
ii
n n n
i i i i
i i i
n a x b y
x a x b x y
(2）
例 经实验测得某物理量 x,y 的 5 对数据如下，用最小二乘法求其拟合曲线,
解之得 9392.60a,5183.1?b,
5 702 758,
702 998 64 108 396,
ab
ab
所以 xy 5 1 8 3.19 3 9.60 为所求的拟合曲线,
解 首先描出散点图，发现其近似地在一直线上，
假设拟合曲线为 bxay,
将所给数据点代入正规方程组（ 2 ）得
141
148
123
125
150
172
123
126
165
187iy
ix
⑴ 若某组数据点 ),,2,1)(,( niyx ii 的散点图近似于一条双曲线，我们就可以设拟合曲线的一般形式为
.1 xbay
如果我们令 z
y
1,t
x
1
,则有 z a b t,此时 ba,可按线 性 最小二乘问题求出；
⑵ 若拟合曲线的一般形式为 ebxya?
将其两边取对数得
(a,b 为待定系数 ).
,lnln bxay 令 Aazy ln,ln,得 bxAz
则 A,b 可按线 性 最小二乘问题求出,
二、将非多项式曲线拟合转化为线性拟合思考题 ⒈根据“让所有点均离其最近”这一原则，数学上是如何刻画的,
⒉用线性最小二乘法求指数拟合 e bxya? 曲线的一般步骤是什么?
一,求积公式的建立二,求积公式的误差估计三,复化求积公式四,变步长的求积公式第四节 数 值 积 分第四节 数 值 积 分用牛顿 - 莱布尼茨公式 )()(d)( aFbFxxf
b
a
来计算定积分
b
a
xxfI d)( 的值的前提是,)( xf 的原函数 )( xF 能够求出,当 )( xf 由表格形式给出，或其原函数不易求出时，用牛顿 -- 莱布尼茨公式就解决不了问题,在此，我们希望用一个易于求原函数的函数（如多项式函数）来近似代替被积函数，然后计算这个多项式函数的积分，就会得到定积分的近似计算公式,
在［ a,b ］区间上选取 n+ １ 个等距节点,khax
k
,
k =0,1,2,?,n,其中
n
ab
h
为步长，并记
kk
yxf?)(,
下面，我们通过数据点组 ),(
kk
yx （ k =0,1,2,? n ）
来构造被积函数 )( xf 的 n 次插值多项式 )( xP
n
，即 一、求积公式的建立
.)()(
0
n
k
kkn xlyxP
其中 )( xl k 为插值基函数，所以
,d)(d)()(
0
n
k
b
a kk
b
a n
b
a xxlyxxPdxxf
称式（ 1 ）为牛顿 - 科茨（ Newton - Cotes ）求积公式,)( nkC
为牛顿 - 科茨系数，为计算 )( nkC,作变量替换
x a t h,于是
ba knk xxlabC d)(1)(
若记?
b
a k
n
k xxlabC d)(
1)(,则有求积公式
n
k
k
n
k
b
a
yCabxxf
0
)()(d)(
(1)
b
a
n
kj
j jk
j x
xx
xx
ab
d1
0
th
jk
jt
nh
n n
kj
j
d1
0
0
即 ()
0
0
( 1 )
[ ( ) ] d
! ( ) !
nk nn
n
k
j
jk
C t j t
n k n k
（ 2 ）
按式（ 2 ）计算科茨系数只涉及多项式函数的积分，
故不会有困难,将科茨系数列于下面,
288
19
288
75
288
50
288
50
288
75
288
19
5
90
7
90
32
90
12
90
32
90
7
4
8
1
8
3
8
3
8
1
3
6
1
6
4
6
1
2
2
1
2
1
1
)( n
k
Cn
当 n =1 时，牛顿－科茨求积公式被称为梯形公式，即
1 ( ) ( ) ( ),
2
T b a f a f b
当 n =2 时，牛顿－科茨公式称为抛物公式（或辛普森公式），即
)(
2
4)()(
6
1
bf
ba
fafabS,当 n =4 时，牛顿－科茨公式简称为科茨公式，即
)(7)(32)(12)(32)(7
90 43210
xfxfxfxfxf
ab
C
,
其中 4,3,2,1,0,4
1)( kabkax
k,
例 1 试分别用梯形公式、抛物公式、科茨公式计算定积分
1
0,5
dxx?,
解 ⑴梯形公式,
1
0,5
1 0,5d ( 0,5 1,0 ) 0,4 2 6 7 7 6 7
2
xx,
⑵抛物公式,
)0.175.045.0(
6
5.01
d
1
5.0
xx
4 3 0 9 3 4 0 3.0?,
(3) 科茨公式，此时步长,1 2 5.04 5.01h
.43 09 64 07.0
)0.17875.032
75.012625.0325.07(
90
5.01
d
1
5.0
xx
按牛顿 — 莱 布 尼茨公式计算，易知 1
5.0
2
31
5.0 3
2d xxx
.4 3 0 9 6 4 1 4.0?
若 )( xP n 为函数 )( xf 的 n 次插值多形项式，则
),()()( xPxfxR nn两边积分得
.d)(d)(d)( ba nbaba n xxPxxfxxR
所以，有,d)( b
a nn xxRIIfR
其中 (,)ab,于是，可得如下余项公式,
(1) 梯形公式的余项
).,(,)(12 )( 3 baabfTIR T
二、求积公式的误差估计
(2) 辛普森公式的余项
),(),()
2
(
180
)4(4 bafababSIR
S?
,
(3) 科茨公式的余项
),(),()
4(9 4 5
)(2 )6(6 bafababCIR
C?
,
从拉格朗日插值多项式的余项公式可以看到,代替被积函数的插值多项式次数越高，对被积函数光滑要求越高,另外，积分区间越小，求积公式的截断误差越小,因此，我们经常先把积分区间分成若干小区间，在每个区间上采用低阶求积公式，如梯形公式或抛物公式，然后再把它加起来得到整个区间上的求积公式，这是复化求积公式的基本思想,
三、复化求积公式
1,复化梯形公式将 区 间 nba ],[ 等 分,步 长
n
ab
h
,分点
khax k ),,2,1,0( nk,在每个小区间 ],[ 1?kk xx 上应用梯形公式，即得复化梯形公式
b
a
n
k
x
x
k
k
xxfxxf
1
0
1 d)(d)(
)]()([2 1
1
0
kk
n
k
xfxfh
1
1
[ ( ) 2 ( ) ( ) ],2 n k
k
h f a f x f b?
于是，有复化梯形求积公式
1
1
[ ( ) 2 ( ) ( ) ]
2
n
nk
k
h
T f a f x f b
（ 3 ）,
2,复化辛普森公式若记子段 ],[
1?kk
xx 的中点为
22
1
h
khax
k
,则复化辛普森公式为
)]()(4)([
6
1
1
0
2
1
kkk
n
k
n
xfxfxf
h
S （ 4 ）
类似可得复化科茨公式，这里从略,
例 2 分别用复化梯形公式，复化辛普森公式计算积分
1
20
4
d,
1
Ix
x
并将计算结果与准确值比较,
解 （ 1 ）准确值为
1 102
0
4 d 4 a r c t a n 3,1 4 1 5 9
1
I x x
x
（ 2 ） 利 用 复 化 梯 形 公 式 计 算,取
1 2 5.0
8
1
8
01
h,节点
88
1
0
k
kx
k
（ k = 0,1,2,? 8 ）,被积函数
2
1
4
)(
x
xf
在节点处的函数值分别为
9 3 8 4 6.3)
8
1
(,4)0( ff
5 0 6 8 5.3)
8
3
(,7 6 4 7 1.3)
8
2
( ff
8 7 6 4 0.2)
8
5
(,2 0 0 0 0.3)
8
4
( ff
2 6 5 4 9.2)
8
7
(,5 6 0 0 0.2)
8
6
( ff
2)1(?f,
8
0,1 2 5 1 2 3 4 5 6 7[ ( 0) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) 2 ( ) ( 1 ) ]2 8 8 8 8 8 8 8T f f f f f f f f f
.13899.3?
（ 3 ） 为了利用上面求出的被积函数在相应节点处的函数值，取 4?n,此时步长
4
01?
h
1
0,2 5
4
,利用复化求积公式（ 4 ）得
8684822878583814)1()0(625.04 fffffffffS
.14159.3?
注意，在计算 4S 时，基本节点 )4,3,2,1,0(
4
0 k
k
x
k
,
中 间 节 点 1
2
1 1 1
4 2 4 4 8
k
kk
x
( 0,1,2,3 )k?,如
8
3
,
8
1
2
1
2
1 10 xx 等,
比较 8T,4S 与 I 的准确值会发现，尽管 8T 与 4S 的计算量基本相同，但 4S 的精度要比 8T 的精度高，可见，
复化辛普森公式是一种精度较高的求积公式,
应用复化求积公式计算定积分，关键是步 长选择,
步长太小，计算量加大，舍入误差也随之增大；步长太大，截断误差增大，精度也难以保证,所以,利用复化求积公式计算时，通常采用变步长的求积方法,即 把步长逐次减半，在步长逐次减半的过程中，反复利用复化求积公式进行计算，直到达到所 需 的精度,
若记 nI 为步长为 h 时的积分值,nI 2 为将 h 分半后的积分值,? 为给定的精度要求，可以证明四 步长的求积公式当 2nnII 时，便有
2
( ) d
b
n
a
f x x I,这样，用二分前后的两个积分值 nI 与 2 nI 之差的绝对值 2nnII? 是否小于预先给定的精度?,来决定二分过程是否停止,
由于用数值积分公式计算定积分，主要工作量在函数值的计算上，因此在应用变步长的复化求积公式时，应考虑步长二分前后计算公式间的关系，以便充分利用已求出的信息（函数值）,
1.变步长梯形法则的递推公式对积分 xxf
b
a
d)(? 应用变步长的复合梯形求积公式,把区间 [ a,b ] 逐次二等分，设现已把 [ a,b ] n 等分，步长为
n
ab
h
,分点为
,,,,,,,,1210 bxxxxxxa nkk
则由复化梯形公式（ 3 ），知
.)()(2)(2
1
1
n
k
kn bfxfaf
hT
然后，将步长
n
ab
h
二分，即在每个子区间
],[ 1 kk xx? 中增加一个二分点 ),,2,1(
2
1 nkx k,便将
[ a,b ] 区间 2 n 等分，此时步长为
n
abh
22
,分点为
,,,,,,,,,,212121 122,110 bxxxxxxxxxa nkkk
则由复化梯形公式，知
1
1
1
1
2
)()(2)(2)(
4 2
1
n
k
n
k
kkn
bfxfxfaf
h
T
n
k
k
n
k
k xf
h
bfxfaf
h
1
1
1
)(
2
)()(2)(
4 2
1
n
k
kn
xf
h
T
1
)(
22
1
2
1,
即即
n
k
knn xf
hTT
1
2 )(22
1
2
1 （ 5 ）
为变步长的梯形递推公式,使用该公式时，要特别注意，
n
ab
h
为二分前的步长,?
n
k
k
xf
1
)(
2
1
为被积函数
)( xf 在二分后所有新增节点处函数值之和,下页图给出了变步长的梯形求积计算框图,
同样的道理，可讨论变步长的辛普森求积法则,
开始输入 a,b,c
a b h
2 )) ( ) ( ( 1 h b f a f T
s? 0,
x? a + h / 2
b? x
ε
T
T
1 2
输出
2 T
Y
N
N
Y
s? s + f ( x ),x = x + h
T 2
h s / 2
h? h/2,
T
1
例 3 用变步长方法计算 xx xI ds i n10,解 先对整个区间 [0,1] 用梯形公式,对于
x
xxf s i n)(?,1s i nl i m)0(
0
x
xf
x
,8 4 1 4 7 1 0.0)1(?f,故
1
1 [ ( 0 ) ( 1 ) ] 0,9 2 0 7 3 5 5,
2T f f
然后将区间 [0,1] 二分，由于 9 5 8 8 5 1 1.0)
2
1
(?f,利用递推公式（ 5 ）得
21
1 1 1
( ) 0,9 3 9 7 9 3 3,
2 2 2
T T f
再将每一个小区间二分一次，并计算新增分点上的函数值,
9 8 9 6 1 5 8.0)
4
1
(?f,9088517.0)
4
3
(?f,再利用公式 （ 5 ）,得
9 4 4 5 1 3 5.0)43()41(2 5.021 24 ffTT 这样不断二分下去，计算结果见 下 表 （表中 k 代表二分次数，区间等分数 kn 2? ）,积分准确值为
0,9 4 6 0 8 3 1,用变步长方法二分 10 次得到了这个结果,
0.94605965
0.9460831100.94598504
0.946083090.94569093
0.946082780.94451352
0.946081570.93979331
0.946076960.9207355 0
k 2kT k 2kT
思考题
1,为什么说抛物公式比梯形公式精确度高？
2,复化梯形公式与变步长的梯形递推公式在算法上有何区别与联系？
一,欧拉方法二,改进的欧拉方法三,龙格 — 库塔方法四,误差的分析第五节 常微分方程的数值解法应用解析方法求解常微分方程初值问题,?
00 )(
),(
d
d
yxy
yxf
x
y
( 1 )
只能求出一些特殊类型的方程 的 解,在大多数情况下，对初值问题（ 1 ），只 需 给出其数值解就能满足工程需要,初值问题 ( 1 ) 的数值解就是对未知函数 )( xfy? 给出在一系列节点 nxxx,,,10? 处的函数值 nyyyy,,,,210?,
第五节 常微分方程的数值解法将方程 ))(,( xyxfy 两边对 x 从 ix 到 1?ix 积分，得 1
1( ) ( ) (,( ) ) d
i
i
x
ii xy x y x f x y x x
( 2 )
用公式（ 2 ）计算函数值 )(
1?i
xy 的难点在于积分
1
d))(,(
i
i
x
x
xxyxf 的计算,下面我们就用不同的数值积分公式代替此积分，得出相应的微分方程数值解公式,
由于?
1
))(,(d))(,(
i
i
x
x
ii xyxhfxxyxf
其中 ii xxh 1 为步长，于是有一、欧拉方法用 iy 替代上式右端的 )( ixy,并将由式（ 3 ）右端算出的值 1?iy 作为 )( 1?ixy 的近似值，这样建立的计算格式
),,2,1,0(),(1 niyxhfyy iiii （ 4 ）
称为欧拉公式,用此公式求解初值问题（ 1 ）的数值解的方法叫做欧拉方法,欧拉方法的特点是算法简单，但精度不高,
))(,()()( 1 iiii xyxhfxyxy （ 3）
为了提高精度，改用梯形公式计算求积项
) ) ](,())(,([
2
d))(,(
11
1
iiii
x
x
xyxfxyxf
h
xxyxf
i
i
,
将上式代入（ 2 ）式得
))](,())(,([
2
)()(
111
iiiiii
xyxfxyxf
h
xyxy,( 5 )
将其中 )(),(
1 ii
xyxy
分别用
ii
yy,
1?
代替，则得到
)].,(),([2 111 iiiiii yxfyxfhyy
称其为梯形公式,该公式的右端含有未知的 1?iy,因此，
称其为隐式格式,为了应用公式（ 5 ）计算 1?iy,我们先用欧拉公式（ 4 ）计算 1?iy 的近似值 1?iy,称其为预二、改进的欧拉方法报值，然后用预报值 1?iy 代替公式（ 5 ）中右端的 1?iy,再计算出校正值 1?iy,这样建立起来的预报 — 校正系统称为改进的欧拉公式，记为
）（校正预报
niyxfyxf
h
yy
yxhfyy
iiiiii
iiii
,,2,1,0)()],(),([
2
)(),(
111
1
（ 6）
为了减少计算量，可将式（ 6 ）写为
)(
2
1
),,2,1,0(),(
),(
1
1
cpi
piic
iiip
yyy
niyxhfyy
yxhfyy
（ 7 ）
请画出改进的欧拉方法的计算框图,
例 1 取步长 h =0.1,求初值问题
d2
,
d
( 0 ) 1
yx
y
xy
y
在 x =0 到 x =1 各节点上的数值解,
解 所给的初值问题的改进的欧拉公式为
2
0.1 ( )
i
p i i
i
x
y y y
y
,
1
2
0,1 ( )
i
c i p
p
x
y y y
y
,
1
1
()
2
i p c
y y y
,
表 2
计算结果见下 表,
i
x 0 0,1 0,2 0,3 0,4 0,5
i
y 1 1,09 5 9 1,18 4 1 1,26 6 2 1,34 3 4 1,41 6 4
()
i
yx 1 1,09 5 4 4 5 1,18 3 2 1 6 1,26 4 9 1 1 1341641 1,41 4 2 1 4
i
x 0 0,6 0,7 0,8 0,9 1,0
i
y
1 1,4 860 1,55 2 5 1,61 6 5 1,67 8 2 1,73 7 9
()
i
yx 1 1,48 3 2 4 0 1,54 9 1 9 3 1,61 2 4 5 2 1,67 3,32 0 1,73 2 0 5 1
该初始值问题的解析解为 xy 21,按这个式子计算出 y 在相应节点处的值（下表第三行是相应的精确值），与准确解比较，会发现按改进的欧拉公式计算出的函数大致有三位有效数字,下面给出的龙格 — 库塔（ Runge — Kutta ）方法具有更高的精度,
对于初值问题（ 1 ），我们不加证明给出如下计算格式,
)22(
6
),(
)
2
,
2
(
)
2
,
2
(
),(
43211
34
23
12
1
kkkk
h
yy
hkyhxfk
k
h
y
h
xfk
k
h
y
h
xfk
yxfk
ii
ii
ii
ii
ii
公式（ 8 ）称为四阶龙格 — 库塔公式，其计算框图 见下页,
（ 8）
三、龙格 — 库塔方法开始结束
Y
输出输入
n h,,y,x
0 0
) 2 hk y,2 h f( x k
2
0 0
3
) 2 hk y
2,h f( x k 1
0 0
2
) hk y,f(x k 3 0
1
4
) y,f ( x k
0 0 1
定义
1 i y ),f ( x,?
x 1 = x
0 + h
6 ) k 2k 2k h(k y y 4 3 2 1 0 1
1 1 y,x
n? i?
N
,
1 0 x x
1 0
y y?
,1 i i
例 2 四阶龙格 — 库塔方法求解例 1中的初值问题,
解 初值问题
2,
( 0) 1
xyy
y
y
(0 ≤ x ≤ 1)
的四阶龙格 — 库塔公式为（步长 2.0?h ）
i
i
i y
xyk 2
1
21
1
0.2
2( )0.2
2,
0.22
2
i
i
i
x
k y k
yk
计算结果见下表直接验证可知,由 四阶龙格 — 库塔方法所得到的结果具有四位有效数字,
ix
iy
0 0.2 0.4 0.6 0.8 1.0
1 1.183229 1.341667 1.483281 1.612514 1.732142
,
2
2.0
)
2
2.0
(2
2
2.0
2
23
ky
x
kyk
i
i
i
3
34 2.0
)2.0(22.0
ky
xkyk
i
i
i?
)22(6 2.0 43211 kkkkyy ii
,
.
若取步长 h =0.05,用四阶龙格 — 库塔方法计算例
1 中 的 初 值 问 题 在 节 点 2.01?x,?2x
0.4,8.0,6.0 43 xx 处函数 )( xy 的值，有
,3 4 1 6 4 1.1,1 8 3 2 1 6.1 21 yy
.6 1 2 4 5 2.1,4 8 3 2 4 0.1 43 yy再将它们与准确值比较，会发现比步长取 0.2 时的解还要精确,
在实际计算中，并不是步长取得愈小愈好，因为，
步长变小，计算量增大，舍入误差也相应地增大，一般地，我们根据事先给定的精度要求，在计算过程中，
逐步缩小步长,具体作法如下,
四、误差的分析设给定的精度要求为 ε,从节点
i
x 出发，先以
h 为步长，利用四阶龙格 — 库塔方法得到 )(
1?i
xy 的近似值为
)(
1
h
i
y
，然后将步长 h 分半，即以
2
h
为步长，从
i
x 出发，经两步按四阶龙格 — 库塔方法计算出 )( 1?ixy 的近似值
)(
1
2
h
i
y
,可以证明,若不等式
)()( 12 hii yy h （ 9）
成立，则将
)(
1
2
h
iy? 作为 )( 1?ixy 的近似值即满足精度要求,否则，进一步将“步长折半”进行计算，直到不等式（ 9 ）成立后，再取最后一个
)(
1
2
h
iy? 作为
)( 1?ixy 的近似值即可,
以上计算过程中，步长自动选取，故称其为变步长的龙格 - 库塔方法,用不等式（ 9 ）也可以实现变步长的改进欧拉方法等,
下页图给出了变步长的四阶龙格 - 库塔方法计算框图，仿该计算框图可画出 变 步长的改进欧拉方法的计算框图,
1 1 1
t t + h?
x=x
0
+h
0
结束开始
b 输入
n,ε,y ),f(x,,y,x 0 0
输出
2
y x,
b? x?
Y
x,
x 0?
2 0
y y?
1 2
y
y
N
N
1
y
2
y?
x? t
1?
2 0
y g?
N
h
0
=(b － x
0
)/n
计算
y
1
1
1 0 1
h
h,g y
2
10
tx?
计算
y
2
Y
0 1 h h
<
Y
计算 y1
计算 y2
k1=f(x0,y0),
k2=f(x0+h0/2,y0+ h0 k1/2),
k3=f(x0+h0/2,y0+ h0 k2/2),
k4=f(x,y0+ h0k3),
y1= y0+( k1+2k2+2k3+ k4) h0/6.
k1=f(t1,g0),
k2=f(t1+h1/2,g0+ h1 k1/2),
k3=f(t1+h1/2,g0+ h1 k2/2),
k4=f(t1 +h1,g0+ h1k3),
y2= g0+( k1+2k2+2k3+ k4) h1/6.
思考题
⒈ 是否每个微分方程都能求其数值解,
⒉ 用欧拉方法能够求出初值问题的足够精确的解吗？
