![数值分析--第4章数值积分与数值微分[1]详解](https://img.doczj.com/img1f/13pvsxt5m9v6zum1c7eb4j4a4yk1fq9m-f1.webp)
![数值分析--第4章数值积分与数值微分[1]详解](https://img.doczj.com/img1f/13pvsxt5m9v6zum1c7eb4j4a4yk1fq9m-d2.webp)
第4章 数值积分与数值微分
1 数值积分的基本概念
实际问题当中常常需要计算定积分。在微积分中，我们熟知，牛顿—莱布尼兹公式是计算定积分的一种有效工具，在理论和实际计算上有很大作用。对定积分()b
a
I f x dx =
?
，若()f x 在区间
[,]a b 上连续，且()f x 的原函数为()F x ，则可计算定积分
()()()b
a
f x dx F b F a =-?
似乎问题已经解决，其实不然。如
1）()f x 是由测量或数值计算给出数据表时，Newton-Leibnitz 公式无法应用。 2）许多形式上很简单的函数，例如
222sin 1
(),sin ,cos ,,ln x x f x x x e x x
-= 等等，它们的原函数不能用初等函数的有限形式表示。
3）即使有些被积函数的原函数能通过初等函数的有限形式表示，但应用牛顿—莱布尼兹公式计算，仍涉及大量的数值计算，还不如应用数值积分的方法来得方便，既节省工作量，又满足精度的要求。例如下列积分
41arc 1)arc 1)1dx tg tg C x ?=+++-+?+? 对于上述这些情况，都要求建立定积分的近似计算方法——数值积分法。
1.1 数值求积分的基本思想
根据以上所述，数值求积公式应该避免用原函数表示，而由被积函数的值决定。由积分中值定理：对()[,]f x C a b ∈，存在[,]a b ξ∈，有
()()()b
a
f x dx b a f ξ=-?
表明，定积分所表示的曲边梯形的面积等于底为b a -而高为()f ξ的矩形面积（图4-1）。问题在于点ξ的具体位置一般是不知道的，因而难以准确算出()f ξ。我们将()f ξ称为区间[,]a b 上的平均高度。这样，只要对平均高度()f ξ提供一种算法，相应地便获得一种数值求积分方法。
如果我们用两端的算术平均作为平均高度()f ξ的近似值，这样导出的求积公式
[()()]2
b a
T f a f b -=
+ (4-1) 便是我们所熟悉的梯形公式（图4-2）。而如果改用区间中点2
a b
c +=的“高度”()f c 近似地取代
平均高度()f ξ，则可导出所谓中矩形公式（简称矩形公式）
()2a b R b a f +??
=- ???
(4-2)
更一般地，我们可以在区间[,]a b 上适当选取某些节点k x ，然后用()k f x 加权平均得到平均高度()f ξ的近似值，这样构造出的求积公式具有下列形式：
()()n
b
k k a
k f x dx A f x =≈∑?
(4-3)
式中k x 称为求积节点；k A 成为求积系数，亦称伴随节点k x 的权。权k A 仅仅与节点k x 的选取有关，而不依赖于被积函数()f x 的具体形式。
这类由积分区间上的某些点上处的函数值的线性组合作为定积分的近似值的求积公式通常称为机械求积公式，它避免了Newton-Leibnitz 公式寻求原函数的困难。对于求积公式(4-3)，关键在于确定节点{}k x 和相应的系数{}k A 。
1.2 代数精度的概念
由Weierstrass 定理可知，对闭区间上任意的连续函数，都可用多项式一致逼近。一般说来，多项式的次数越高，逼近程度越好。这样，如果求积公式对m 阶多项式精确成立，那么求积公式的误差仅来源于m 阶多项式对连续函数的逼近误差。因此自然有如下的定义
定义4.1 如果某个求积公式对于次数不超过m 的多项式均准确地成立，但对于1m +次多项式就不准确成立，则称该求积公式具有m 次代数精度。
例1 判断求积公式
1
1
1
()[58(0)5(9
f x dx f f f -≈++? 的代数精度。
解 记
111
()()()[58(0)5(9
I f f x dx I f f f f -==++?，
因为
1
1
1
(1)2(1)(585)29
I dx I -===
++=?， 1
11
()()[58
05(09
I x xdx I x -==??+?=?=0，
图4-1 图4-2
1222112
()()(50.6805
0.6)93I x x
dx I x -==?+?+?=?2=，3
13
333311()()[505(]09I x x dx I x -==?++?=?=0，
1444112
()()(50.36
050.36)95I
x x dx I x -==?++?=?2=，5
15555511
()()[505(]09I x x dx I x -==?++?=?=0，
16
6633112()()[5(0.6)05(0.6)]0.2497
I x x dx I x -==?++?=≠?2=，7
所以求积公式具有5次代数精度。
1.3插值型的求积公式
最直接自然的一种想法是用()f x 在[,]a b 上的插值多项式()n x ?代替()f x ，由于代数多项式的原函数是容易求出的，我们以()n x ?在[,]a b 上的积分值作为所求积分()I f 的近似值，即
()()b
n a
I f x dx ?≈?
这样得到的求积分公式称为插值型求积公式。通常采用Lagrange 插值。
设[,]a b 上有1n +个互异节点01,,
,n x x x ，()f x 的n 次Lagrange 插值多项式为
()()()n n k k k L x l x f x ==∑
其中0()n i
k j k i
j k
x x l x x x =≠-=
-∏，插值型求积公式为 0
()()()n
b
n k k a
k I f L x dx A f x =≈=∑? (4-4)
其中(), 0,1,
,b
k k
a A l x dx k n =
=?
。可看出，{}k A 仅由积分区间[,]a b 与插值节点{}k x 确定，与
被积函数()f x 的形式无关。求积公式(4-4)的截断误差为
(1)1()
()()()()(1)!
n b
b
b
n n n a
a
a
f R f f x dx L x dx x dx n ξω++=-=+???
(4-5)
定义4.2 求积公式
()()n
b
k k a
k f x dx A f x =≈∑?
如其系数()b
k k
a A l x dx =
?
，则称此求积公式为插值型求积公式。
定理4.1 形如(4-3)的求积公式至少有n 次代数精度的充分必要条件是插值型的。
证明 如果求积公式(4-3)是插值型的，由公式(4-5)可知，对于次数不超过n 的多项式()f x ，其余项[]R f 等于零，因而这时求积公式至少具有n 次代数精度。
反之，如果求积公式(4-3)至少具有n 次代数精度，那么对于插值基函数()k l x 应准确成立，并注意到()k j jk l x δ=，即有
()()n
b
k
j k j k a j l x dx A l x A ===∑?
所以求积公式(4-3)是插值型的。
1.4 求积公式的收敛性与稳定性
定义4.3 在求积公式(4-3)中，若
0lim ()()n
b
k k a
n k h A f x f x dx →∞=→=∑?
其中11max()i i i n
h x x -≤≤=-，则称求积公式(4-3)是收敛的。
实际使用任何求积公式时，除截断误差外，还有舍入误差，因此我们必须研究其数值稳定性。在求积公式(4-3)中，由于计算()k f x 可能产生误差k δ，实际得到k f ，即()k k k f x f δ=+，记
()(),()n
n
n k k n k k k k I f A f x I f A f ====∑∑
如果对任给正数0ε>，只要误差k δ充分小就有
()()()n
n n k
k k k I f I f A
f x f ε=??-=
-≤??
∑ (4-6) 它表明求积公式(4-3)计算是稳定的，由此给出：
定义4.4 对任给0ε>，若存在0δ>，只要() (0,1,,)k k f x f k n δ-≤=就有(4-6)成立，
则称求积公式(4-3)是稳定的。
定理4.2 若求积公式(4-3)中系数0 (0,1,,)k A k n >=，则此求积公式是稳定的；若k A 有正有
负，计算可能不稳定。
证明 对任给0ε>，若取b a
ε
δ=
-，对0,1,
,k n =都有()k k f x f δ-≤，则有
()()()()n
n
n
n n k
k k k k k k k k k I f I f A
f x f A f x f A δ===??-=
-≤-≤??
∑∑∑
注意对任何代数精度0≥的求积公式均有
(1)1n
b
k
n a
k A
I dx b a ====-∑?
可见0k A >时，有
()()()n
n
n n k k k k I f I f A A b a δδδε==-≤==-=∑
∑
由定义4.4可知求积公式(4-3)是稳定的。
若k A 有正有负时，假设(())0k k k A f x f ->，且()k k f x f δ-=，有
00
()()()()()n
n
n n k
k k k k k
k k n n
k k k k I f I f A
f x f A f x f A A b a δδδ
====??-=
-=-??
=>=-∑∑∑∑
它表明初始数据的误差可能会引起计算结果误差的增大，即计算可能不稳定。
2 Newton-Cotes 公式
2.1 Cotes 系数
被积函数在积分区间内变化平缓，可用等距节点插值公式近似。将积分区间[,]a b 划分为n 等分，步长b a
h n
-=
，等距节点,0,1,,k x a kh k n =+=。此时求积公式(4-4)中的积分系数可得到简化 00()()n n b b b j k k a a a
j j k j j k
j k
x x x a jh
A l x dx dx dx x x k j h
==≠≠---===--∏∏???
作变换x a th =+，则有
000000
()(1)(1)()()()()!()!!()!n k n k n n n
n
n n k j j j j k
j k
j k
t j h h b a A hdt t j dt t j dt k j h k n k k n k n --===≠≠≠----==-=----∏∏∏??? 令
()
00
(1)()!()!n k
n n n k
j j k
C
t j dt k n k n -=≠-=--∏? 则()
()n k k A b a C =-，求积公式(4-4)可简化为
()0
()()()n
n k k k I f b a C f x =≈-∑ (4-7)
称为n 阶Newton-Cotes 公式，简记为N-C 公式，{}
()
n k C 称为Cotes 系数。
由()n k C 的表达式可看出，它不但与被积函数无关，而且与积分区间也无关。因此可将Cotes 系
数事先列成表格供查用（见表4-1）。
N-C 公式的截断误差为
(1)2(1)
000
()()()()()(1)!(1)!n n n n
b
n n n j a
j j f h R f x x dx f t j dt n n ξξ+++===-=-++∏∏?
? (4-8) 1n =时
[]11()()()()()()222b a
I f b a f a f b f a f b -??=-+=+????
(4-9)
为梯形公式
2n =时
141()()()()()()4()()662662a b b a a b I f b a f a f f b f a f f b +-+????
=-++=++????????
(4-10) 为辛普生公式。
4n =时
3()()7()32()12()32()7()90424b a b a b a b a I f f a f a f f a f b --+-??
=
++++++????
(4-11)
为Cotes 公式。
表4-1
从表4-1可看出，当8n =时出现了负系数，由定理4.2可知，实际计算中将使舍入误差增大，并且往往难以估计。从而Newton-Cotes 公式的收敛性和稳定性得不到保证，因此实际计算中不用高阶Newton-Cotes 公式。
2.2 偶阶求积公式的代数精度
作为插值型的求积公式，n 阶的牛顿-柯特斯公式至少具有n 次的代数精度。求积公式的代数精度能否进一步提高呢？
定理4.3 当阶为偶数时，N C -公式(4-7)至少具有1n +次代数精度。 证明 我们只要验证，当n 为偶数时，N-C 公式对1
()n f x x +=的余项为零。按余项公式(4-8)，
由于这里(1)
()(1)!n f
x n +=+，从而
()()n
b
n j
a
j R f x x dx ==-∏?
引进变换x a th =+，并注意到j x a jh =+，有
2
00
()()n
n
n n j R f h
t j dt +==-∏?
当n 为偶数，则
2n 为整数，再令2
n
t u =+，进一步有 2
2222
022
()()()02n n
n n n n n n
n j j n n R f h u j du h u j du ++--==-=+-=-=∏∏??
因为被积函数为奇函数。 □
2.3 几种低阶求积公式的余项
梯形求积公式的余项为
()
()()2!
b
T a
f R I T x a x b dx ξ''=-=--?
由于()()x a x b --在[,]a b 上不变号，利用积分中值定理有
3
()()()()()(,)212
b T a f f R x a x b dx b a a b ηηη''''=
--=--∈?， (4-12) Simpson 公式的余项为
[]()()4()()6
b
S a
b a
R I S f x dx f a f c f b -=-=-
++? 这里2
a b
c +=
。构造次数不超过3的多项式()H x ，使满足 ()(),()(),()(),()()H a f a H c f c H c f c H b f b ''====
由于Simpson 公式具有三次代数精度，它对于这样构造的三次式()H x 是准确的，即
[]()()4()()6
b
a
b a
H x dx H a H c H b -=
++?
所以
[]()()b
S a
R f x H x dx =-?
由第二章的例6，可知
(4)
21()()()()()()4!
f x H x f x a x c x b ξ-=
--- 因2()()()x a x c x b ---在[,]a b 上保号，应用积分中值定理有
(4)
24
(4)1()()()()4!
()(,)1802b s a R f x a x c x b dx
b a b a f a b ηηη=
-----??=-
∈ ???
?， (4-13)
3 复化求积公式
前面导出的误差估计式表明，用N-C 公式计算积分近似值时，步长越小，截断误差越小。但缩小步长等于增加节点数，亦即提高插值多项式的次数，Runge 现象表明，这样并不一定能提高精度。理论上已经证明，当n →∞时，N-C 公式所求得的近似值不一定收敛于积分的准确值，而且随着n 的增大，N-C 公式是不稳定的。因此，实际中不采用高阶N-C 公式，为提高计算精度，可考虑对被积函数用分段低次多项式插值，由此导出复化求积公式。
3.1 复化梯形公式
将区间[,]a b 划分为n 等分，分点,,0,1,,k b a
x a kh h k n n
-=+=
=，在每个区间1[,]k k x x +
(0,1,
,1)k n =-上采用梯形公式，则得
[]1
1
1
10
()()()()()2k k
n n b
x k k n a
x k k h I f x dx f x dx f x f x R f +--+=====++∑∑??
(4-14)
记
[]11
101()()()2()()22n n n k k k k k h h T f x f x f a f x f b --+==??
=+=++????
∑∑ (4-15)
称为复化梯形公式，其余项为
3211
100
()1()()()(,)1212n n n n k k k k k k k h b a h R f I T f f x x n ηηη--+==-''''=-=--∈∑∑=， 由于2()[,]f x C a b ∈，且
1
01010
1min ()()max ()n k k k k n k n k f f f n ηηη-≤≤-≤≤-=''''''≤≤∑ 所以存在(,)a b η∈使
1
1()()n k k f f n ηη-=''''=∑
于是复化梯形公式余项为
2
()()()12
n b a R f h f η-''=-
(4-16) 从式(4-16)可以看出，余项误差是2
h 阶，所以当2()[,]f x C a b ∈，有
lim ()b
n a
n T f x dx →∞
=?，
即复化梯形公式是收敛的。事实上只要()[,]f x C a b ∈，则可得收敛性，因为由(4-15)得
1011()()() ()2n n b n k k a k k b a b a T f x f x f x dx n n n -==--??
=+→→∞????
∑∑? 所以复化梯形公式(4-15)收敛。此外，n T 的求积系数为正，由定理4.2知复化梯形公式是稳定的。
3.2 复化辛普森公式
将区间[,]a b 划分为n 等分，在每个区间1[,]k k x x +上采用辛普森公式，记121
2
k k x x h +=+
则得 1
1
1
110
()()()4()()()6k k
n n b
x k k k n a
x k k h I f x dx f x dx f x f x f x R f +--++==??===+++??∑∑??
(4-17) 记
11
1201()4()2()()6n n n k k k k h S f a f x f x f b --+==??
=+++????
∑∑ (4-18)
称为复化辛普森求积公式，其余项由(4-13)得
4
1(4)
10
1()()
(,)1802n n n k k k k k h R f I S h f x x ηη-+=??=-=-∈ ???∑， 于是当4
()[,]f x C a b ∈时，与复化梯形公式相似有
4
(4)
()()(,)1802n n b a h R f I S f a b ηη-??=-=-∈ ???
， (4-19)
可以看出误差阶是4
h ，收敛性是显然的。事实上，只要()[,]f x C a b ∈，则有
111200114()()()6() ()
n n n n k k
k k k k b
a b a b a b a S f x f x f x n n n f x dx n --+===---??=++????→→∞∑∑∑? 此外，由于n S 中求积系数均为正数，故知复化辛普森求积公式计算稳定。
例2 根据函数表4-1
表4-1
用复化梯形公式和复化辛普森公式计算1
0sin x
I dx x =
?的近似值，并估计误差。
解 由复化梯形公式
7
1
1[(0)(1)2()]0.945691168k k
I f f f =≈++=∑
由复化辛普森公式
34
11
121[(0)(1)2()4()]0.9460841648k k k k I f f f f ==-≈+++=∑∑
与准确值0.9460831
I =比较，显然用复化Simpson 公式计算精度较高。
为了利用余项公式估计误差，要求sin ()x
f x x =
的高阶导数，由于 10sin ()cos()x
f x xt dt x
==?
所以有
1
1()
00()cos()cos()2
k k k
k d k f
x xt dt t xt dt dx π==+?? 于是
1
1()
001
1
max ()cos()21k k
k x k f
x t xt dt t dt k π≤≤=+≤=
+?
? 由复化梯形误差公式(4-16)得
2
28801111
()max ()0.000434121283
x h R f I T f x ≤≤??''=-≤≤= ???
由复化辛普森误差公式(4-19)得
4
6
44111()0.2711018085
R f I S -??=-≤=? ???
例3 若用复化求积分公式计算积分
1
x I e dx -=?
的近似值，要求计算结果有四位有效数字，n 应取多大？
解 因为当01x ≤≤时，有
10.31x e e -≤≤≤
于是
1
0.31x e dx -<
要求计算结果有四位有效数字，即要求误差不超过
41
102
-?。又因为 ()() 1 [0,1]k x
f x e x -=≤∈
由式(4-16)得
22411()1012122
T h R h f ξ-''=≤=?
即41
106
n ≥
?，开方得40.8n ≥。因此若用复化梯形公式求积分，n 应等于41才能达到精度。 若用复化Simpson 公式，由式(4-19)
44
4(4)41111()10180218016180162
S h h R f n ξ-????=≤=≤? ? ?
?????? 即得 1.62n ≥。故应取2n =。
4龙贝格求积公式
4.1梯形公式的逐次分半算法
如前所述，复化求积公式的截断误差随着步长的缩小而减少，而且如果被积函数的高阶导数容易计算和估计时，由给定的精度可以预先确定步长，不过这样做常常是很困难的，一般不值得推崇。实际计算时，我们总是从某个步长出发计算近似值，若精度不够可将步长逐次分半以提高近似值，直到求得满足精度要求的近似值。
设将区间[,]a b 分为n 等分，共有1n +个分点，如果将求积区间再二分一次，则分点增至21n +个，我们将二分前后两个积分值联系起来加以考虑。注意到每个子区间1[,]k k x x +经过二分只增加了
一个分点1
12
1
()2k k k x
x x ++=
+，用复化梯形公式求得该子区间上的积分值为 121()2()()4
k k k h f x f x f x ++??++?? 注意，这里b a
h n
-=代表二分前的步长，将每个子区间上的积分值相加得
[]11
211200
()()()42n n n k k k k k h h T f x f x f x --++===++∑∑
即
1
2120
1()22n n n k k h T T f x -+==+∑ (4-20)
这表明，将步长由h 缩小为
2
h
时，2n T 等于n T 的一半再加新增加节点处的函数值乘以当前步长。 算法4.1
1．输入,,(),a b f x ε
2．置01,,[()()]2
b a
m h T h f a f b -==
=+ 3．置0F =，对11,2,,2m k -=
((21))F F f a k h =++-
4．01
2
T T hF =+
5．若03T T ε-<，输出I T ≈，停机；否则01,,2
h
m m h T T +???，转3。
4.2 李查逊（Richardson ）外推法
假设用某种数值方法求量I 的近似值，一般地，近似值是步长h 的函数，记为1()I h ，相应的误
差为
12112()k p p p k I I h h h h ααα-=++
++
(4-21)
其中12(1,2,
),0i k i p p p α=<<<<<
是与h 无关的常数。若用h α代替(4-21)中的h ，则得
121
1
22
11212()()()()k k
k
p p p k p p p p p p k I I h h h h h h h ααααααααααααα-=++
++=++
++
(4-22)
式(4-22)减去式(4-21)乘以1
p α
，得
13
3
2
1
2
1
11123()[()]
()()()k
k
p p p p p p p p p p k I I h I I h h h h αααααααααα
α---=-+-+
+-+
取α满足1α≠，以11p
α-除上式两端，得
1321
1123()()1k p p p p k p I h I h I b h b h b h ααα
--=++++- (4-23)
其中112())(2,3,)i p p p
i b i αααα=--=仍与h 无关。令
11
112()()
()1p p I h I h I h ααα
-=- 由式(4-23)，以2()I h 作为I 的近似值，其误差至少为2()p
O h ，因此2()I h 收敛于I 的速度比1()I h 快。
不断重复以上作法，可以得到一个函数序列
11
11()()
()2,3,1m m p m m m p I h I h I h m ααα-----==
-， (4-24) 以()m I h 近似I ，误差为()()m
p m I I h O h -=。随着m 的增大，收敛速度越来越快，这就是Richardson
外推法。
4.3 龙贝格求积公式
由前面知道，复化梯形公式的截断误差为2()O h 。进一步分析，我们有如下欧拉—麦克劳林（Euler-Maclaurin ）公式：
定理4.4 设()[,]f x C a b ∞∈，则有
24212()k k I T h h h h ααα-=++
++
其中系数(1,2,
)k k α=与h 无关。
把李查逊外推法与欧拉—麦克劳林公式相结合，可以得到求积公式的外推算法。特别地，在外推算法式(4-24)中，取1
,22
k p k α=
=，并记0()()T h T h =，则有 114()()
2(),1,2,41
j m m m m
h
T T h T h m ---==
- (4-25)
经过(1,2,
)m m =次加速后，余项便取下列形式：
2(1)2(2)12()m m m T h I h h δδ++=++
+
(4-26)
上述处理方法通常称为李查逊（Richardson ）外推加速方法。
为研究Romberg 求积方法的机器实现，引入记号：以()0k T 表示二分k 次后求得的梯形值，且以
()
k m
T 表示序列{}()0k T 的m 次加速值，则依以上递推公式得到 ()
(1)()
1141(),1,2,4141
m k k k m
m m m m
T
h T T k +--=-=-
-
称为龙贝格求积算法。
Romberg 公式的计算过程见下表4-2
表4-2
算法4.2
(1) 输入,,(),a b f x ε
(2) 置[](0)
0,()(),12
h
h b a T f a f b k =-=
+= (3) 计算
1
2
()
(1)001
11()[(())]22k k k i T h T h f a i h --==++-∑
对1,,j k =
(1)()
11
()441
j k j k j j j k j j
j
T T T
-+-----=
-
(4) (0)
(0)
1
k
k T
T
ε--<，输出(0)()b
k
a
T
f x dx ≈?停机；否则
,12
h
h k k ?+?，返回(3)。 例4 用Romberg 算法计算积分1
320
I x dx =
?
。
解 利用逐次分半算法(4-20)和Romberg 算法(4-25)，计算结果见表4-3。
(0)01
[(0)(1)]0.5000002T f f =+=
(1)(0)0011
0.5()0.42677722T T f =+?=
(2)(1)00113
0.25[()()]0.407018244
T T f f =+?+=
(3)
(2)0010.251357[()()()()]0.401812228888
T T f f f f =+?+++=
表4-3
k
()0k T ()1k T ()2k T ()3k T ()4k T ()5k T
0 0.500000
1 0.426777 0.402369
2 0.407018
3 0.401812 0.400077 0.40005
4 0.400050 4 0.400463 0.400009 0.400009 5
0.400118 0.400002 0.400002 0.400002
5 高斯求积公式
5.1 一般理论
等距节点的插值型求积公式，虽然计算简单，使用方便，但是这种节点等距的规定却限制了求积公式的代数精度。试想如果对节点不加限制，并适当选择求积系数，可能会提高求积公式的精度。Gauss 型求积公式的思想也正如此，亦即在节点数n 固定时，适当地选取节点{}k x 与求积系数{}k A ，使求积分公式具有最高精度。
设有1n +个互异节点01,,
,n x x x 的机械求积分公式
()()()n
b
k k a
k x f x dx A f x ρ=≈∑?
(2-27)
具有m 次代数精度，那么有取() (0,1,
,)l f x x l m ==，式(1)精确成立，即
()() (0,1,,)n
b l
l
j j a
j A x x x dx l m ρ===∑? (2-28)
式(2)构成1m +阶的非线性方程组，且具有22n +个未知数, (0,1,,)k k x A k n =，
所以当()x ρ给定后，只要122m n +≤+，即21m n ≤+时，方程组有解。这表明式1n +个节点的求积公式的代
数精度可达到21n +。
另一方面，对式(1)，不管如何选择{}k x 与{}k A ，最高精度不可能超过21n +。事实上，对任意的互异节点0{}n k k x =，令
222
222101()()()()
()n n n p x x x x x x x x ω++==---
有
220
()0n
k n k k A P x +==∑，然而21()()0b
n a
x P x dx ρ+>?。
定义4 如果求积分公式(4-27)具有21n +次代数精度，则称这组节点{}k x 为Gauss 点，相应公式(4-27)称为带权()x ρ的高斯求积公式。
定理5 插值型求积公式的节点01n a x x x b ≤<<<≤是高斯点的充分必要条件是以这些节
点为零点的多项式
101()()()
()n n x x x x x x x ω+=---
与任何次数不超过n 的多项式()P x 带权正交，即
1()()()0b
n a
x P x x dx ρω+=?
(4-29)
证明 必要性。设()n P x H ∈，则121()()n n P x x H ω++∈，因此，如果01,,,n x x x 是高斯点，则
式(1)对于1()()()n f x P x x ω+=精确成立，即有
110
()()()()()0n
b
n k k n k a
k x P x x dx A P x x ρωω++===∑?
故式(4-29)成立。
再证充分性。对于21()n f x H +?∈，用1()n x ω+除()f x ，记商为()P x ，余式为()q x ，即
1()()()()n f x P x x q x ω+=+，其中(),()n P x q x H ∈，由式(4-29)可得
()()()()b
b
a
a
x f x dx x q x dx ρρ=?
? (4-30)
由于所给求积公式(4-27) 是插值型的，它对于()n q x H ∈是精确成立的，即
()()()n
b
k k a
k x f x dx A q x ρ==∑?
再注意到1()0 (0,1,
,)n k x k n ω+==，知()() (0,1,
,)k k q x f x k n ==，从而由(4-30)有
()()()()()n
b
b
k k a
a
k x f x dx x q x dx A f x ρρ===∑?
?
可见求积公式(4-27)对一切次数不超过21n +的多项式精确成立，因此(0,1,
,)k x k n =为高斯点。
证毕。
定理表明在[,]a b 上关于权()x ρ的正交多项式系中的1n +次多项式的零点就是求积公式(4-27)
的高斯点。因此，求Gauss 点等价于求[,]a b 上关于权()x ρ的1n +次正交多项式的1n +个实根。有了求积节点(0,1,
,)k x k n =后，可如下确定求积系数
()()()n
b
k j k j k a
j x l x dx A l x A ρ===∑?
其中0()n
j
k j k
j
j k
x x l x x
x =≠-=
-∏。
下面讨论高斯求积公式的余项。设在节点(0,1,,)k x k n =上()f x 的21n +次Hermite 插值多
项式为()H x ，即
212
1()(),()(),0,1,,n k k n k k H x f x H x f x k n ++''===
由Hermite 余项公式
(222
1()()()()(22)!
n n f f x H x x n ξω++-=+
有
[]00(222
1()()()()
()()()
()()()()()()()()()()(22)!
n
b
k k a
k n
b
k k a
k b
b
a
a
b
a
n b
n a
R f x f x dx A f x x f x dx A H x x f x dx x H x dx x f x H x dx
f x x dx
n ρρρρρξρω==++=-=-=-=-=+∑?∑?????
定理6 高斯求积公式的求积系数(0,1,,)k A k n =全是正的。
证明 由于具有高斯节点(0,1,,)k x k n =的高斯求积公式具有21n +次代数精度，所以对于多
项式0() (0,1,,)n
j
k j k
j
j k
x x l x k n x
x =≠-=
=-∏，公式准确成立，即
220
()()() (0,1,,)n
b
k
j k j k a
j x l x dx A l x A k n ρ====∑?
推论 高斯求积公式是稳定的。
定理7 设()[,]f x C a b ∈，高斯求积公式是收敛的，即
lim ()()()n
b
k k a
n k A f x x f x dx ρ→∞
==∑?
6 数值微分
在微积分学里，求函数()f x 的导数()f x '一般来讲是容易办到的，但若所给函数()f x 由表格给出，则()f x '就不那么容易了，这种对列表函数求导数通常称为数值微分。
6.1 利用差商求导数
最简单的数值微分公式是用差商近似代替导数 1）向前差商
000()()
()f x h f x f x h +-'≈
(4-31)
2）向后差商
000()()
()f x f x h f x h
--'≈
(4-32)
2）中心差商
000()()
()2f x h f x h f x h
+--'≈
(4-33)
在几何图形上，这三种差商分别表示弦AB 、AC 和BC 的斜率，将这三条弦同过A 点的切线AT 相比较，从图1可以看出，一般地说，以BC 的斜率更接近于切线AT 的斜率0()f x '，因此就精度而言，(3)式更为可取，称
00()()
()2f x h f x h G h h
+--=
(4-34)
为求0()f x '的中点公式。
现在考察式(4-34)计算近似导数0()f x '所产生的截断误差，首先分别将0()f x h ±在0x 处作Talor 展开，有
0()f x h ±
0 0 0
x
图1 差商示意图
1、解：将)(x V n 按最后一行展开，即知)(x V n 是n 次多项式。 由于 n i i i n n n n n i n x x x x x x x x x x V ...1...1... ......... ...... 1 )(21110 20 0---= ，.1,...,1,0-=n i 故知0)(=i n x V ，即110,...,,-n x x x 是)(x V n 的根。又)(x V n 的最高 次幂 n x 的系数为 )(...1...1... ...... .........1),...,,(101 1 21 11 2 2221 02001101j n i j i n n n n n n n n n n n x x x x x x x x x x x x x x V -== ∏-≤<≤-----------。 故知).)...()()(,...,,()(1101101------=n n n n x x x x x x x x x V x V 6、解：（1）设 .)(k x x f =当n k ,...,1,0=时，有.0)()1(=+x f n 对 )(x f 构造Lagrange 插值多项式， ),()(0 x l x x L j n j k j n ∑== 其 0)()! 1() ()()()(1)1(=+=-=++x w n f x L x F x R n n n n ξ， ξ介于j x 之间，.,...,1,0n j = 故 ),()(x L x f n =即 .,...,1,0,)(0 n k x x l x k j n j k j ==∑= 特别地，当0=k 时， 10) (=∑=n j x j l 。 （2） 0)()1(1) ()1()()(0000=-=??? ? ??-??? ? ??-=--=-===∑∑∑∑k j j i j i k j k i i j i i k j n j k i i j k n j j x x x x i k x l x x i k x l x x ）利用（。 7、证明：以b a ,为节点进行线性插值，得 )()()(1 b f a b a x a f b a b x x P --+--= 因 0)()(==b f a f ，故0)(1=x P 。而 ))()(("2 1 )()(1b x a x f x P x f --= -ξ，b a <<ξ。 故)("max )(8 122)("max )(max 2 2 x f a b a b x f x f b x a b x a b x a ≤≤≤≤≤≤-=??? ??-≤。 14、解：设 ))...()(()(21n n x x x x x x a x f ---=， k x x g =)(，记)() (1 ∏=-=n j j n x x x w ，则 ),()(x w a x f n n =).()(' j n n j x w a x f = 由差商的性质知 [])! 1()(1,..,,1) (' 1 )(')('1 211 11 -== ==-===∑∑∑ n g a x x x g a x w x a x w a x x f x n n n n n j j n k j n n j j n n k j n j j k j ξ， ξ介于n x x ,...,1之间。 当20-≤≤ n k 时，0)()1(=-ξn g ， 当 1-=n k 时，)!1()(1-=-n g n ξ， 故 ???-=-≤≤=-= --=∑1,,20,0)!1()(1) ('1 11 n k a n k n g a x f x n n n n j j k j ξ 16、解：根据差商与微商的关系，有 [] 1! 7! 7!7)(2,...,2,2)7(7 10===ξf f ， [ ] 0! 80 !8)(2,...,2,2)8(8 1 ===ξf f 。 （ 13)(47+++=x x x x f 是7次多项式， 故 ,!7)()7(=x f 0)()8(=x f ）。 25、解：（1） 右边= [][]dx x S x f x S dx x S x f b a b a ??-+-)(")(")("2)(")("2 = [] d x x S x f x S x S x S x f x f b a ?-++-)("2)(")("2)(")(")("2)(" 222 = [] d x x S x f b a ?-)(")(" 22 = [][]dx x S dx x f b a b a 2 2 )(")("??- =左边。 （2）左边= ? -b a dx x S x f x S ))(")(")(("
数值分析 第二章 2．当1,1,2x =-时，()0,3,4f x =-,求()f x 的二次插值多项式。 解： 0120121200102021101201220211,1,2, ()0,()3,()4;()()1 ()(1)(2)()()2()()1 ()(1)(2) ()()6 ()()1 ()(1)(1) ()()3 x x x f x f x f x x x x x l x x x x x x x x x x x l x x x x x x x x x x x l x x x x x x x ==-===-=--==-+-----==------= =-+-- 则二次拉格朗日插值多项式为 2 20 ()()k k k L x y l x ==∑ 0223()4() 14 (1)(2)(1)(1)23 537623 l x l x x x x x x x =-+=---+ -+= +- 6．设,0,1,,j x j n =L 为互异节点，求证： （1） 0()n k k j j j x l x x =≡∑ (0,1,,);k n =L （2）0 ()()0n k j j j x x l x =-≡∑ (0,1,,);k n =L 证明 （1） 令()k f x x = 若插值节点为,0,1,,j x j n =L ，则函数()f x 的n 次插值多项式为0 ()()n k n j j j L x x l x == ∑。 插值余项为(1)1() ()()()()(1)! n n n n f R x f x L x x n ξω++=-= + 又,k n ≤Q
(1)()0 ()0 n n f R x ξ+∴=∴= 0()n k k j j j x l x x =∴=∑ (0,1,,);k n =L 0 000 (2)()() (())()()(()) n k j j j n n j i k i k j j j i n n i k i i k j j i j x x l x C x x l x C x x l x =-==-==-=-=-∑∑∑∑∑ 0i n ≤≤Q 又 由上题结论可知 ()n k i j j j x l x x ==∑ ()()0 n i k i i k i k C x x x x -=∴=-=-=∑原式 ∴得证。 7设[]2 (),f x C a b ∈且()()0,f a f b ==求证： 21 max ()()max ().8 a x b a x b f x b a f x ≤≤≤≤''≤- 解：令01,x a x b ==，以此为插值节点，则线性插值多项式为 10 101010 ()() ()x x x x L x f x f x x x x x --=+-- =() () x b x a f a f b a b x a --=+-- 1()()0()0 f a f b L x ==∴=Q 又 插值余项为1011 ()()()()()()2 R x f x L x f x x x x x ''=-= -- 011 ()()()()2 f x f x x x x x ''∴= --
《计算机数学基础》数值部分第五单元辅导 14 常微分方程的数值解法 一.重点内容 1. 欧拉公式： ）心知1）a 儿+1 =儿 + hfg ，儿） m 1、 伙=0丄2,…川一 1） I 无=x Q +kh 局部截断误差是0（*）。 2. 改进欧拉公式： 预报一校正公式： 预报值 _v*+1 =儿+ hf （x k ,儿） - h - 校正值 y M = y k +-［f （x kt y k ） + /（x A+1, y M ）］ 即 儿+1 =儿+ ￡ "（忑'儿）+心+「儿+ hfg ,儿））］ 或表成平均的形式： 儿=儿+ hfg ,儿） '儿=儿+"（无+】，儿） +K ） 改进欧拉法的局部截断误差是0（2） 3. 龙格一库塔法 二阶龙格一库塔法的局部截断误差是0（爪） 三阶龙格一库塔法的局部截断误差是0（护） 四阶龙格F 塔法公式：儿计=儿+ 2（匕+ 2心+ 2? + ?） 四阶龙格一库塔法的局部截断误差是0（爪）。 二实例 y' = — y — xv f2（0 < x < 0.6） 例1用欧拉法解初值问题｛ ' ? -取步长/匸02计算过程保留 b （o ）= 1 4位小数。 解/i=0.2. f （x ）= —y —xy 2<,首先建立欧拉迭代格式 y*+i =儿+ hf g,y k ） = y k -hy k -hx k y ； =0.2 儿（4 - x k y k ）（k = 0,1,2） K 2=f(x n +^h, yk+-hK\)t gg+舟人，>'n +y/?A3）；
当k=0, xi=0.2 时，已知x（）=0,y（）=l,有 y(0?2)今i=0?2X l(4-0X 1)=0.8000 当k=\. M=0?4时，已知“=0?2」尸0?8,有 y(0?4)今2=0.2 X 0.8X(4-0.2X0.8)=0.614 4 当k=2, xs=0.6 时，已知x2=0.4,y2=0.6144,有 y(0?6)今3=0.2 X0.6144X (4-0.4 X 0.4613)=0.8000 「J, ,2 ?_ ZX 例2用欧拉预报一校正公式求解初值问题\y + v +V sinx=,取步长/?=0.2,计算 .y ⑴=1 y(0.2),y(0.4)的近似值，计算过程保留5位小数。 解步长力=0.2,此时/(x,y)=—y—fsiiu 欧拉预报一校正公式为： 预报值兀I = y k + hfg y k) - I J_ 校正值)3=儿+尹(忑，儿)+ fg,儿+1)] 有迭代格式] 预报值儿+] = y k 4-h(-y k -y； sin x k) =y k (0?8-0?2儿sin x k) < h 、—— 2 校止值y如]=儿 +尸[(一片一力sinxJ + LN+i-yl sin.v I+1)] ——?> =儿(°?9一0?1儿sin心)一0?1(儿+| +y；j sin心利) 、"M=0.別=1」)=1 时，Xj=1.2> 有 儿=yo(°?8-O?2yo sinx0) = 1 x (0.8-02x lsin 1) = 0.63171 y(1.2) ?= lx(0.9-0.1xlxsinl)-0.1(0.63171+0.631712sinl.2) = 0.71549 当 T xi=1.2, yi=0.71549 时,x2=1.4,有 y2 =儿(0.8-0?2儿sinXj) = 0.71549x(0.8-02x0.71549sinl.2) =0.47697 y(14) z y2 = 0.71549x(0.9-0.1x0.71549xsin 1.2)-0.1(0.47697+ 0.476972 sin 1.4) =0.52608 V = 8 — 3y 例3写出用四阶龙格一库塔法求解初值问题^ ‘的计算公式，取步长/匸0.2计 b(0) = 2 算y(0.4)的近似值。讣算过程保留4位小数。 解此处.心,刃=8 —3”四阶龙格一库塔法公式为 艰=儿 + % + 2? + 2勺 + ?) 1 h, y n+ y/?A3)： 本例计算公式为: 0 2 呱严儿+三(32?+2?+心
数值分析习题集 （适合课程《数值方法A 》和《数值方法B 》） 长沙理工大学 第一章 绪 论 1. 设x >0,x 的相对误差为δ,求ln x 的误差. 2. 设x 的相对误差为2％,求n x 的相对误差. 3. 下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指出 它们是几位有效数字: *****123451.1021,0.031,385.6,56.430,7 1.0.x x x x x =====? 4. 利用公式(3.3)求下列各近似值的误差限: ********12412324(),(),()/,i x x x ii x x x iii x x ++其中**** 1234 ,,,x x x x 均为第3题所给的数. 5. 计算球体积要使相对误差限为1％,问度量半径R 时允许的相对误差限是多少? 6. 设028,Y =按递推公式 1n n Y Y -=( n=1,2,…) 计算到100Y .27.982(五位有效数字),试问计算100Y 将有多大误差? 7. 求方程2 5610x x -+=的两个根,使它至少具有四位有效数字27.982). 8. 当N 充分大时,怎样求2 1 1N dx x +∞+?? 9. 正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝2 ? 10. 设 212S gt = 假定g 是准确的,而对t 的测量有±0.1秒的误差,证明当t 增加时S 的绝对 误差增加,而相对误差却减小. 11. 序列 {}n y 满足递推关系1101n n y y -=-(n=1,2,…),若0 1.41y =≈(三位有效数字), 计算到 10y 时误差有多大?这个计算过程稳定吗? 12. 计算6 1)f =, 1.4≈,利用下列等式计算,哪一个得到的结果最好? 3 -- 13. ()ln(f x x =,求f (30)的值.若开平方用六位函数表,问求对数时误差有多大?若
第四章 数值积分与数值微分 1.确定下列求积公式中的特定参数，使其代数精度尽量高，并指明所构造出的求积公式所具有的代数精度： 101210121 12120 (1)()()(0)(); (2)()()(0)(); (3)()[(1)2()3()]/3; (4)()[(0)()]/2[(0)()]; h h h h h f x dx A f h A f A f h f x dx A f h A f A f h f x dx f f x f x f x dx h f f h ah f f h -----≈-++≈-++≈-++''≈++-?? ?? 解： 求解求积公式的代数精度时，应根据代数精度的定义，即求积公式对于次数不超过m 的多项式均能准确地成立，但对于m+1次多项式就不准确成立，进行验证性求解。 （1）若101(1) ()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1012h A A A -=++ 令()f x x =，则 110A h Ah -=-+ 令2 ()f x x =，则 3 221123 h h A h A -=+ 从而解得 011431313A h A h A h -?=?? ? =?? ?=?? 令3 ()f x x =，则 3()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++=
令4()f x x =，则 455 1012()5 2 ()(0)()3 h h h h f x dx x dx h A f h A f A f h h ---== -++=? ? 故此时， 101()()(0)()h h f x dx A f h A f A f h --≠-++? 故 101()()(0)()h h f x dx A f h A f A f h --≈-++? 具有3次代数精度。 （2）若 21012()()(0)()h h f x dx A f h A f A f h --≈-++? 令()1f x =，则 1014h A A A -=++ 令()f x x =，则 110A h Ah -=-+ 令2 ()f x x =，则 3 2211163 h h A h A -=+ 从而解得 1143 8383A h A h A h -?=-?? ? =?? ?=?? 令3 ()f x x =，则 22322()0h h h h f x dx x dx --==? ? 101()(0)()0A f h A f A f h --++=
数值分析 2?当x=1,—1,2时，f(x)=O, 一3,4,求f(x)的二次插值多项式。解： X 0 =1,x j = — 1,x 2 = 2, f(X。)= 0, f (xj = -3, f (x2)= 4; l o(x)=(x-xi^~x2\=-1(x 1)(x-2) (x o -X/X o _x2) 2 (x -x0)(x -x2) 1 l i(x) 0 2(x-1)(x-2) (x i ~x0)(x i ~x2) 6 (x—x0)(x—x,) 1 l2(x) 0 1(x-1)(x 1) (X2 -X°)(X2 - X i) 3 则二次拉格朗日插值多项式为 2 L 2(X)= ' y k 1 k ( x) kz0 = -3l°(x) 4l2(x) 1 4 =(x_1)(x—2) 4 (x-1)(x 1) 2 3 5 2 3 7 x x - 6 2 3 6?设Xj, j =0,1,||(,n 为互异节点，求证： n (1 )7 x：l j(x) =x k(k =0,1川,n); j=0 n (2 )7 (X j -x)k l j(x)三0 (k =0,1川,n); j ￡ 证明 (1)令f(x)=x k
n 若插值节点为X j, j =0,1,|l(, n，则函数f (x)的n次插值多项式为L n(x)八x k l j(x)。 j=0 f (n 十)(?) 插值余项为R n(X)二f(X)-L n(X) n1(X) (n +1)!
.f(n1)( ^0 R n(X)=O n 二瓦x k l j(x) =x k(k =0,1川,n); j ：o n ⑵、(X j -x)k l j(x) j卫 n n =為(' C?x j(—x)k_L)l j(x) j =0 i =0 n n i k i i =為C k( -x) (、X j l j(x)) i =0 j=0 又70 _i _n 由上题结论可知 n .原式二''C k(-x)k_L x' i=0 =(X -X)k =0 -得证。 7设f (x) c2 la,b 1且f (a) =f (b)二0,求证: max f(x)兰一(b-a) max a $至小一*丘f (x). 解：令x^a,x^b，以此为插值节点，则线性插值多项式为 L i(x^ f(x o) x x f (xj X o —人x -X o X —X o x-b x-a ==f(a) f(b)- a - b x -a 又T f (a) = f (b)二0 L i(x) = 0 1 插值余项为R(x)二f (x) - L,(x) f (x)(x - X Q)(X - xj 1 f(x) = 2 f (x)(x -X g)(X -xj
（适合课程《数值方法A 》和《数值方法B 》） 第一章 绪 论 1. 设x >0,x 的相对误差为δ,求ln x 的误差. 2. 设x 的相对误差为2％,求n x 的相对误差. 3. 下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指出它们是几位 有效数字: ***** 123451.1021,0.031,385.6,56.430,7 1.0.x x x x x =====? 4. 利用公式(3.3)求下列各近似值的误差限: * * * * * * * * 12412324(),(),()/,i x x x ii x x x iii x x ++其中* * * * 1234,,,x x x x 均为第3题所给的数. 5. 计算球体积要使相对误差限为1％,问度量半径R 时允许的相对误差限是多少? 6. 设028,Y =按递推公式 11783 100 n n Y Y -=- ( n=1,2,…) 计算到100Y .若取783≈27.982(五位有效数字),试问计算100Y 将有多大误差? 7. 求方程2 5610x x -+=的两个根,使它至少具有四位有效数字(783≈27.982). 8. 当N 充分大时,怎样求 2 11N dx x +∞+?? 9. 正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝2 ? 10. 设2 12S gt = 假定g 是准确的,而对t 的测量有±0.1秒的误差,证明当t 增加时S 的绝对误差增加, 而相对误差却减小. 11. 序列{}n y 满足递推关系1101 n n y y -=-(n=1,2,…),若02 1.41y =≈(三位有效数字),计算到10 y 时误差有多大?这个计算过程稳定吗? 12. 计算6 (21)f =-,取 2 1.4≈,利用下列等式计算,哪一个得到的结果最好? 3 6 3 11,(322), ,9970 2. (21) (322) --++ 13. 2 ()ln(1)f x x x =- -,求f (30)的值.若开平方用六位函数表,问求对数时误差有多大?若改用另一等 价公式 2 2 ln(1)ln(1)x x x x - -=-+ + 计算,求对数时误差有多大? 14. 试用消元法解方程组{ 10 10 12121010; 2. x x x x +=+=假定只用三位数计算,问结果是否可靠? 15. 已知三角形面积 1sin , 2 s ab c = 其中c 为弧度, 02c π << ,且测量a ,b ,c 的误差分别为,,.a b c ???证 明面积的误差s ?满足 . s a b c s a b c ????≤ ++ 第二章 插值法 1. 根据( 2.2)定义的范德蒙行列式,令
习 题 一 解 答 1．取3.14，3.15， 227，355113 作为π的近似值，求各自的绝对误差，相对误差和有效数字的位数。 分析：求绝对误差的方法是按定义直接计算。求相对误差的一般方法是先求出绝对误差再按定义式计算。注意，不应先求相对误差再求绝对误差。有效数字位数可以根据定义来求，即先由绝对误差确定近似数的绝对误差不超过那一位的半个单位，再确定有效数的末位是哪一位，进一步确定有效数字和有效数位。有了定理2后，可以根据定理2更规范地解答。根据定理2，首先要将数值转化为科学记数形式，然后解答。 解：（1）绝对误差: e(x)=π－3.14＝3.14159265…－3.14＝0.00159…≈0.0016。 相对误差： 3()0.0016 ()0.51103.14r e x e x x -==≈? 有效数字： 因为π＝3.14159265…=0.314159265…×10，3.14＝0.314×10，m=1。 而π－3.14＝3.14159265…－3.14＝0.00159… 所以│π－3.14│＝0.00159…≤0.005=0.5×10－2＝21311 101022 --?=? 所以，3.14作为π的近似值有3个有效数字。 （2）绝对误差: e(x)=π－3.15＝3.14159265…－3.14＝－0.008407…≈－0.0085。 相对误差： 2()0.0085 ()0.27103.15r e x e x x --==≈-? 有效数字： 因为π＝3.14159265…=0.314159265…×10，3.15＝0.315×10，m=1。 而π－3.15＝3.14159265…－3.15＝－0.008407… 所以│π－3.15│＝0.008407……≤0.05=0.5×10－1 ＝11211101022 --?=? 所以，3.15作为π的近似值有2个有效数字。 （3）绝对误差: 22 () 3.14159265 3.1428571430.0012644930.00137e x π=-=-=-≈- 相对误差：
i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
第四版 数值分析习题 第一章 绪 论 1. 设x >0,x 的相对误差为δ,求ln x 的误差. 2. 设x 的相对误差为2％,求n x 的相对误差. 3. 下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指 出它们是几位有效数字: *****123451.1021,0.031,385.6,56.430,7 1.0.x x x x x =====? 4. 利用公式求下列各近似值的误差限: ********12412324(),(),()/,i x x x ii x x x iii x x ++其中**** 1234 ,,,x x x x 均为第3题所给的数. 5. 计算球体积要使相对误差限为1％,问度量半径R 时允许的相对误差限是多少? 6. 设028,Y =按递推公式 1n n Y Y -=…) 计算到100Y .(五位有效数字),试问计算100Y 将有多大误差? 7. 求方程2 5610x x -+=的两个根,使它至少具有四位有效数字. 8. 当N 充分大时,怎样求 2 11N dx x +∞ +? ? 9. 正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝2 ? 10. 设 212S gt = 假定g 是准确的,而对t 的测量有±秒的误差,证明当t 增加时S 的绝对误 差增加,而相对误差却减小. 11. 序列 {}n y 满足递推关系1101n n y y -=-(n=1,2,…),若0 1.41y =≈(三位有效数字), 计算到 10y 时误差有多大?这个计算过程稳定吗? 12. 计算61)f =, 1.4≈,利用下列等式计算,哪一个得到的结果最好? 3 -- 13. ()ln(f x x =,求f (30)的值.若开平方用六位函数表,问求对数时误差有多大?若改用另一等价公式 ln(ln(x x =- 计算,求对数时误差有多大?
第四版 数值分析习题 第一章绪论 1.设x>0,x得相对误差为δ,求得误差、 2.设x得相对误差为2％,求得相对误差、 3.下列各数都就是经过四舍五入得到得近似数,即误差限不超过最后一位得半个单位,试指 出它们就是几位有效数字: 4.利用公式(3、3)求下列各近似值得误差限: 其中均为第3题所给得数、 5.计算球体积要使相对误差限为1％,问度量半径R时允许得相对误差限就是多少? 6.设按递推公式 ( n=1,2,…) 计算到、若取≈27、982(五位有效数字),试问计算将有多大误差? 7.求方程得两个根,使它至少具有四位有效数字(≈27、982)、 8.当N充分大时,怎样求? 9.正方形得边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝? 10.设假定g就是准确得,而对t得测量有±0、1秒得误差,证明当t增加时S得绝对误差增 加,而相对误差却减小、 11.序列满足递推关系(n=1,2,…),若(三位有效数字),计算到时误差有多大?这个计算过程 稳定吗? 12.计算,取,利用下列等式计算,哪一个得到得结果最好? 13.,求f(30)得值、若开平方用六位函数表,问求对数时误差有多大?若改用另一等价公式 计算,求对数时误差有多大? 14.试用消元法解方程组假定只用三位数计算,问结果就是否可靠? 15.已知三角形面积其中c为弧度,,且测量a ,b ,c得误差分别为证明面积得误差满足 第二章插值法 1.根据(2、2)定义得范德蒙行列式,令 证明就是n次多项式,它得根就是,且 、 2.当x= 1 , -1 , 2 时, f(x)= 0 , -3 , 4 ,求f(x)得二次插值多项式、 3. 4., 研究用线性插值求cos x 近似值时得总误差界、
第一章 题12 给定节点01x =-，11x =，23x =，34x =，试分别对下列函数导出拉格朗日插值余项： （1） （1） 3 ()432f x x x =-+ （2） （2） 4 3 ()2f x x x =- 解 （1）(4) ()0f x =， 由拉格朗日插值余项得(4)0123() ()()()()()()0 4!f f x p x x x x x x x x x ξ-=----=； （2）(4) ()4!f x = 由拉格朗日插值余项得 01234! ()()()()()() 4! f x p x x x x x x x x x -= ----(1)(1)(3)(4)x x x x =+---. 题15 证明：对于()f x 以0x ，1x 为节点的一次插值多项式()p x ，插值误差 012 10()()()max () 8x x x x x f x p x f x ≤≤-''-≤. 证 由拉格朗日插值余项得 01() ()()()()2!f f x p x x x x x ξ''-= --，其中01x x ξ≤≤， 01 0101max ()()()()()()()() 2!2!x x x f x f f x p x x x x x x x x x ξ≤≤''''-=--≤-- 01210()max () 8x x x x x f x ≤≤-''≤. 题22 采用下列方法构造满足条件(0)(0)0p p '==，(1)(1)1p p '==的插值多项式 ()p x ： （1） （1） 用待定系数法； （2） （2） 利用承袭性，先考察插值条件(0)(0)0p p '==，(1)1p =的插值多项式 ()p x . 解 （1）有四个插值条件，故设230123()p x a a x a x a x =+++，2 123()23p x a a x a x '=++， 代入得方程组001231123010231 a a a a a a a a a =? ?+++=?? =? ?++=? 解之，得01230 021 a a a a =??=?? =??=-?
实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60;
x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下： t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27
第一章 绪论（12） 1、设0>x ，x 的相对误差为δ，求x ln 的误差。 [解]设0*>x 为x 的近似值，则有相对误差为δε=)(*x r ，绝对误差为**)(x x δε=，从而x ln 的误差为δδεε=='=* ****1)()(ln )(ln x x x x x ， 相对误差为* * ** ln ln ) (ln )(ln x x x x r δ εε= = 。 2、设x 的相对误差为2%，求n x 的相对误差。 [解]设*x 为x 的近似值，则有相对误差为%2)(*=x r ε，绝对误差为**%2)(x x =ε，从而n x 的误差为n n x x n x n x x n x x x ** 1 *** %2%2) ()()()(ln * ?=='=-=εε， 相对误差为%2) () (ln )(ln *** n x x x n r == εε。 3、下列各数都是经过四舍五入得到的近似数，即误差不超过最后一位的半个单位，试指出它们是几位有效数字： 1021.1*1=x ，031.0*2=x ，6.385*3=x ，430.56*4=x ，0.17*5 ?=x 。 [解]1021.1*1 =x 有5位有效数字；0031.0* 2=x 有2位有效数字；6.385*3=x 有4位有效数字；430.56* 4 =x 有5位有效数字；0.17*5?=x 有2位有效数字。 4、利用公式（3.3）求下列各近似值的误差限，其中* 4*3*2*1,,,x x x x 均为第3题所给 的数。 （1）* 4*2*1x x x ++； [解]3 334* 4*2*11** *4*2*1*1005.1102 1 10211021)()()()()(----=?=?+?+?=++=? ??? ????=++∑x x x x x f x x x e n k k k εεεε； （2）* 3*2 *1x x x ；
1第一章 习题解答 1 设x >0，x 的相对误差限为δ，求 ln x 的误差。 解：设 x 的准确值为x *，则有 ( | x – x * | /|x *| ) ≤ δ 所以 e (ln x )=| ln x – ln x * | =| x – x * | ×| (ln x )’|x=ξ·≈ ( | x – x * | / | x *| ) ≤ δ 另解： e (ln x )=| ln x – ln x * | =| ln (x / x *) | = | ln (( x – x * + x *)/ x *) | = | ln (( x – x * )/ x * + 1) |≤( | x – x * | /|x *| ) ≤ δ 2 设 x = – 2.18 和 y = 2.1200 都是由准确值经四舍五入而得到的近似值。求绝对误差限ε( x ) 和 ε( y ) 。 解：| e (x ) | = |e (– 2.18)|≤ 0.005，| e (y ) | = |e ( 2.1200)|≤ 0.00005，所以 ε( x )=0.005， ε( y ) = 0.00005。 3 下近似值的绝对误差限都是 0.005，问各近似值有几位有效数字 x 1=1.38，x 2= –0.0312，x 3= 0.00086 解：根据有效数字定义，绝对误差限不超过末位数半个单位。由题设知，x 1，x 2， x 3有效数末位数均为小数点后第二位。故x 1具有三位有效数字，x 2具有一位有效数字，x 3具有零位有效数字。 4 已知近似数x 有两位有效数字，试求其相对误差限。 解：| e r (x ) | ≤ 5 × 10– 2 。 5 设 y 0 = 28，按递推公式 y n = y n-1 – 783/ 100 ( n = 1，2，…) 计算到y 100。若取≈78327.982 （五位有效数字），试问，计算 y 100 将有多大的误差？ 解：由于初值 y 0 = 28 没有误差，误差是由≈78327.982所引起。记 x = 27.982，783?=x δ。则利用理论准确成立的递推式 y n = y n-1 – 783/ 100 和实际计算中递推式 Y n = Y n-1 – x / 100 (Y 0 = y 0) 两式相减，得 e ( Y n ) = Y n – y n = Y n-1 – y n-1 – ( x – 783)/ 100 所以，有 e ( Y n ) = e ( Y n-1) – δ / 100 利用上式求和 δ?=∑∑=?=100111001)()(n n n n Y e Y e 化简，得 e ( Y 100) = e ( Y 0) – δ = δ 所以，计算y 100 的误差界为 4100105001.05.0)(?×=×=≤δεY 6 求方程 x 2 – 56x + 1 = 0的两个根，问要使它们具有四位有效数字，D=ac b 42 ?至少要取几位有效数字？ 如果利用韦达定理，D 又应该取几位有效数字？ 解：在方程中，a = 1，b = – 56，c = 1，故D=4562?≈55.96427，取七位有效数字。
第一章绪论 1.1 引言 常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具。微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力工具，在17到18世纪，在力学、天文、科学技术、物理中，就已借助微分方程取得了巨大的成就。1864年Leverrer根据这个方程预见了海王星的存在，并确定出海王星在天空中的位置。现在，常微分方程在许多方面获得了日新月异的应用。这些应用也为常微分方程的进一步发展提供了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。 研究常微分方程常用数值解是数学工作者的一项基本的且重要的工作。在国内外众多数学家的不懈努力，使此学科基本上形成了一套完美的体系。微分方程的首要问题是如何求一个给定方程的通解或特解。到目前为止，人们已经对许多微分方程得出了求解的一般方法。由于在生产实际和科学研究中所遇到的微分方程问题比较复杂，使这些问题的解即使能求出解析表达式，也往往因计算量太大而难于求出，而对于一些典型的微分方程则可以运用基本方法求出其解析解，并可以根据初值问题的条件把其中的任意常数确定下来。 由于求通解存在许多困难，人们就开始研究带某种定解条件的特解。首先是Cauchy对微分方程初始解的存在惟一性进行了研究。目前解的存在惟一性、延拓性、大范围的存在性以及解对初始解和参数的延续性和可微性等理论问题都已发展成熟。与此同时，人们开始采取各种近似方法来求微分方程的特解，例如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程的近似解。最后，由于当代高科技的发展为数学的广泛应用和深入研究提供了更好的手段。用计算机结合Matlab软件求方程的精确解、近似解，对解的性态进行图示和定性、稳定性研究都十分方便有效。 本章先介绍常微分的一般概念、导出微分方程的一些典型例子及求解微分方程的思路分析。从而得到常微分方程的常用数值解法。