![数值分析--第4章数值积分与数值微分[1]](https://img.doczj.com/imgec/1rojmonp82pa8rlwlouq098q9pplku36-c1.webp)
![数值分析--第4章数值积分与数值微分[1]](https://img.doczj.com/imgec/1rojmonp82pa8rlwlouq098q9pplku36-f2.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 -≈++? 的代数精度。
解 记
1
11()()()[58(0)5(9
I f f x dx I
f f f f -==++? ， 因为
1
11(1)2(1)(585)29
I dx I
-===++=? ，
1
11()()[5805(09
I x xdx I
x -==??+?=? =0，
图4-1 图4-2
1222112
()()(50.68050.6)93I x x dx I x -==?+?+?=? 2=，3
1
3333311()()[505(]09I x x dx I
x -==?++?=? =0， 1444112
()()(50.36050.36)95I x x dx I x -==?++?=? 2=，5
1
5555511()()[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 δ=+ ，记 0
()(),()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 δ-≤ ，则有 0
()()()()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.9460831I = 比较，显然用复化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 ，则得
12112211212()()()()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
1
1123()[()]
()()()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 =，令
2222
22101()()()()()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 差商示意图
第六章习题解答 2、利用梯形公式和Simpson 公式求积分2 1 ln xdx ? 的近似值，并估计两种方法计算值的最大 误差限。 解：①由梯形公式： 21ln 2 ()[()()][ln1ln 2]0.3466222 b a T f f a f b --= +=+=≈ 最大误差限 3''2 ()111 ()()0.0833******** T b a R f f ηη-=-=≤=≈ 其中，(1,2)η∈ ②由梯形公式： 13()[()4()()][ln14ln()ln 2]0.38586262 b a b a S f f a f f b -+= ++=++≈ 最大误差限 5(4)4()66 ()()0.0021288028802880 S b a R f f ηη-=-=≤≈， 其中，(1,2)η∈。 4、推导中点求积公式 3''()()()()() ()224 b a a b b a f x dx b a f f a b ξξ+-=-+<
《计算机数学基础》数值部分第五单元辅导 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?+心
第四章 数值积分与数值微分 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 --++=
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 。
第六章习题解答 2 2、利用梯形公式和 Simpson 公式求积分 ln xdx 的近似值， 并估计两种方法计算值的最大 1 误差限。 解：①由梯形公式： T ( f ) b a [ f (a) f (b)] 2 1 [ln1 ln 2] ln 2 0.3466 2 2 2 最大误差限 R ( f ) (b a)3 f '' ( ) 1 1 1 0.0833 T 12 12 2 12 12 其中， (1,2) ②由梯形公式： b a 4 f ( b a f (b)] 1 4ln( 3 ln 2] 0.3858 S( f ) [ f (a) ) [ln1 ) 6 2 6 2 最大误差限 R S ( f ) (b a)5 f (4) ( ) 6 6 0.0021， 2880 2880 4 2880 其中， (1,2) 。 4、推导中点求积公式 f ( x)dx (b a) f ( a b ) (b a) 3 (a b) b a 2 24 证明： 构造一次函数 P （ x ），使 P a 2 b f a b , P ' ( a b ) f ' ( a b ), P '' ( x) 0 2 2 2 则，易求得 P( x) f ' ( a b )( x a b ) f ( a b ) 2 2 2 且 P(x)dx f ' ( a b )( x a b ) f ( a b ) dx b b a a 2 2 2 f ( a b )dx (b a) f ( a b ) ，令 P(x)dx I ( f ) b b a 2 2 a 现分析截断误差：令 r ( x) f ( x) P(x) f ( x) f ' ( a b )( x a b ) f ( a b ) 2 2 2 由 r ' ( x) f ' (x) f ' ( a b ) 易知 x a 2 b 为 r (x) 的二重零点， 2 a b )2 ， 所以可令 r (x) ( x)( x 2
第六章 习 题 1. 计算下列矩阵的1A ，2A ，A ∞三种范数。 （1）1101A ???=????，（2）312020116A ????=??????? . 2. 用Jacobi 方法和Gauss-Seidel 迭代求解方程组 1231231 238322041133631236x x x x x x x x x ?+=??+?=??++=? 要求取(0)(0,0,0)T x =计算到(5)x ，并分别与精确解(3,2,1)T x =比较。 3. 用Gauss-Seidel 迭代求解 12312312 35163621122x x x x x x x x x ??=??++=???+=?? 以(0)(1,1,1)T x =?为初值，当(1)() 310k k x x +?∞?<时，迭代终止。 4. 已知方程组121122,2,x x b tx x b +=?? +=? （1）写出解方程组的Jacobi 迭代矩阵，并讨论迭代收敛条件。 （2）写出解方程组的Gauss-Seidel 迭代矩阵，并讨论迭代收敛条件. 5. 设有系数矩阵 122111221A ?????=?????? ， 211111112B ?????=??????? ， 证明：（1）对于系数矩阵A ，Jacobi 迭代收敛，而Gauss-Seidel 迭代不收敛. （2）对于矩阵B ，. 6. 讨论方程组 112233302021212x b x b x b ?????????????=??????????????????? 用Jacobi 迭代和Gauss-Seidel 迭代的收敛性；如果都收敛，比较哪种方法收敛更快.
实验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
习 题 六 解 答 1、在区间［0，1］上用欧拉法求解下列的初值问题，取步长h=0.1。 (1)210(1)(0)2y y y '?=--?=?(2)sin (0)0x y x e y -'?=+?=? 解：（1）取h=0.1，本初值问题的欧拉公式具体形式为 21(1)(0,1,2,)n n n y y y n +=--= 由初值y 0=y(0)=2出发计算，所得数值结果如下： x 0=0，y 0=2； x 1=0.1，2100(1)211y y y =--=-= x 2=0.2，2211(1)101y y y =--=-= 指出： 可以看出，实际上求出的所有数值解都是1。 （2）取h=0.1，本初值问题的欧拉公式具体形式为 21(sin )(0,1,2,)n x n n n y y h x e n -+=++= 由初值y 0=y(0)=0出发计算，所得数值结果如下： x 0=0，y 0=0； x 1=0.1， 02 1000 (sin )00.1(sin 0)00.1(01)0.1x y y h x e e -=++=+?+=+?+= x 2=0.2， 122110.1 (sin )0.10.1(sin 0.1)0.10.1(0.10.9)0.2 x y y h x e e --=++=+?+=+?+= 指出： 本小题的求解过程中，函数值计算需要用到计算器。 2、用欧拉法和改进的欧拉法(预测－校正法)求解初值问题，取步长h=0.1。 22(00.5) (0)1 y x y x y '?=-≤≤? =? 解：(1) 取h=0.1，本初值问题的欧拉公式具体形式为 2 1(2)(0,1,2,)n n n n y y h x y n +=+-= 由初值y 0=y(0)=1出发计算，所得数值结果如下：
第一章绪论 1.1 引言 常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具。微分方程的理论和方法从17世纪末开始发展起来，很快成了研究自然现象的强有力工具，在17到18世纪，在力学、天文、科学技术、物理中，就已借助微分方程取得了巨大的成就。1864年Leverrer根据这个方程预见了海王星的存在，并确定出海王星在天空中的位置。现在，常微分方程在许多方面获得了日新月异的应用。这些应用也为常微分方程的进一步发展提供了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。 研究常微分方程常用数值解是数学工作者的一项基本的且重要的工作。在国内外众多数学家的不懈努力，使此学科基本上形成了一套完美的体系。微分方程的首要问题是如何求一个给定方程的通解或特解。到目前为止，人们已经对许多微分方程得出了求解的一般方法。由于在生产实际和科学研究中所遇到的微分方程问题比较复杂，使这些问题的解即使能求出解析表达式，也往往因计算量太大而难于求出，而对于一些典型的微分方程则可以运用基本方法求出其解析解，并可以根据初值问题的条件把其中的任意常数确定下来。 由于求通解存在许多困难，人们就开始研究带某种定解条件的特解。首先是Cauchy对微分方程初始解的存在惟一性进行了研究。目前解的存在惟一性、延拓性、大范围的存在性以及解对初始解和参数的延续性和可微性等理论问题都已发展成熟。与此同时，人们开始采取各种近似方法来求微分方程的特解，例如求微分方程数值解的Euler折线法、Runge-Kutta法等，可以求得若干个点上微分方程的近似解。最后，由于当代高科技的发展为数学的广泛应用和深入研究提供了更好的手段。用计算机结合Matlab软件求方程的精确解、近似解，对解的性态进行图示和定性、稳定性研究都十分方便有效。 本章先介绍常微分的一般概念、导出微分方程的一些典型例子及求解微分方程的思路分析。从而得到常微分方程的常用数值解法。
第六章非线性方程的数值解法习题解答 填空题： 1. 求方程()x f x =根的牛顿迭代格式是__________________。 Ans:1()1()n n n n n x f x x x f x +-=- '- 2.求解方程 在(1, 2)内根的下列迭代法中， (1) (2) (3) (4) 收敛的迭代法是（A ）. A ．(1)和(2) B. (2)和(3) C. (3)和(4) D. (4)和(1) 3.若0)()(，故迭代发散。 以上三中以第二种迭代格式较好。 2、设方程()0f x =有根，且'0()m f x M <≤≤。试证明由迭代格式1()k k k x x f x λ+=- (0,1,2,)k =L 产生的迭代序列{}0k k x ∞ =对任意的初值0(,)x ∈-∞+∞，当2 0M λ<< 时，均收敛于方程的根。
31.（1）解：x (k+1)=（I-αA ）x (k)+αb 因此迭代矩阵B= I-αA=1?3α?2α?α1?2α 其特征值λ1=1-4α；λ2=1-α； 因为迭代收敛的充要条件是：ρ(B)<1，等价于|λ1|<1,且|λ2|<1 因此0<α<1/2时，迭代矩阵收敛。当谱半径最小时，收敛速度最大，等价于|λ1|和|λ2|在(0,1/2)上的最小值，求得当α=3/5时，谱半径最小，迭代速度最快。 (2)已知对称正定矩阵A 的最大和最小特征值分别为：λ1和λn ，且二者均大于0. 不妨设A 的特征值为λ，特征向量为x 。则(I-αA)x=x-αAx=(1-αλ)x 即迭代矩阵B= I-αA 的特征值为1-αλ,而且迭代矩阵收敛的充要条件是ρ(B)<1，即矩阵B 的特征值绝对值最大的要小于1，即|1-αλ1|<1且|1-αλn |<1，求得0<α<2(λ1)-1； 当谱半径最小，即在区间0<α<2(λ1)-1内，|1-αλn |和|1-αλ1|<1的最小值，求得α=2/(λ1+λn ) 32. 证明：因为SOR 方法的迭代矩阵L w =(D-wL)-1[(1-w)D+wU] 假设L w 存在一个特征值λ满足|λ|>1。 所以det(λI-L w )=0 所以det((D-wL)-1[(1-w)D+wU])=det(D-wL)-1*det[(1-λ-1(1-w))D-wL-λ-1wU]=0 (1) 由于矩阵A 是严格对角占优矩阵，所以a ii >0，因此(D-wL)-1的主对角
元素大于零，即det(D-wL)-1>0。 而由于|λ|>1，且0
本科生课程设计报告实习课程数值分析 学院名称管理科学学院 专业名称信息与计算科学 学生姓名 学生学号 指导教师 实验地点 实验成绩 二〇一六年六月二〇一六年六 摘要 常微分方程数值解法是计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.
关键词：数值解法；常微分
1. 实验内容和要求 常微分方程初值问题 有精确解2()cos(2)x y x x e x -=+。 要求：分别取步长h = ，，，采用改进的Euler 方法、4阶经典龙格-库塔R －K 方法和4阶Adams 预测-校正方法计算初值问题。以表格形式列出10个等距节点上的计算值和精确值，并比较他们的计算精确度。其中多步法需要的初值由经典R-K 法提供。运算时，取足以表示计算精度的有效位。 2. 算法说明 函数及变量说明 表1 函数及变量定义 1、欧拉方法： 1()()(,())i i i i y x y x hf x y x +=+ （i=0，1，2，3，......n-1） (0)y a = (其中a 为初值)
2、改进欧拉方法： ~ 1~111()()(,()) ()()[(,())(,())] 2 (0)i i i i i i i i i i y x y x hf x y x h y x y x f x y x f x y x y a ++++=+=++=（i=0，1，2......n-1） (其中a 为初值) 3、经典K-R 方法： 1 1213 243()6 (,)(,)22(,)22 (,) i i i i i i i i i i h y y K f x y h h K f x y K h h K f x y K K f x h y hK +? =+?? =???=++?? ? =++??=++?? 4、4阶adams 预测-校正方法 预测： 校正： Adsms 内插外插公式联合使用称为Adams 预测-校正系统，利用外插公式计算预测，用内插公式进行校正。计算时需要注意以下两点： 1、外插公式为显式，内插公式为隐式。故用内插外插公式时需要进行迭代。 2、这种预测-校正法是四步法，计算Yn+1时，不但用到前一步信息，而且要用到更前三步信息1-n f ，2-n f 3-n f ，因此它不是自动开始的，实际计算时必须借助某种单步法，用Runge-Kutta 格式为它提供初始值 3 21,,y y y ，依据局部截断误差公式得：
9.令*()(21),[0,1]n n T x T x x =-∈，试证*{()}n T x 是在[0,1] 上带权()x ρ=的正交多项式，并求****0123(),(),(),()T x T x T x T x . 证明： 1 1 * *0 1 1 * *011**0 ()()()(21)(21)211()()()()()2()()()()()()()()n m n m n m n m n m n n m n m x T x T x dx x T x dx t x x T x T x dx t T t dt t T t dt T x x T x T x dx t T t ρρρ---=--=-== = ???? ?令，则 由切比雪夫多项式1 01=02 m n dt m n m n ππ ≠??? =≠??==??? 所以*{()}n T x 是在[0,1] 上带权()x ρ= *00*11* 22 2 2*33233()(21)1()(21)21 ()(21)2(21)188()(21)4(21)3(21)3248181 T x T x T x T x x T x T x x x x T x T x x x x x x =-==-=-=-=--=-=-=---=-+- 14.已知实验数据如下： 用最小二乘法求形如2y a bx =+的经验公式，并求均方误差 解： 法方程为
22222(1,)(1,1)(1,)(,)(,1)(,)a y x b x y x x x ?????? =???? ?????? ?? 即 5 5327271.453277277699369321.5a b ??????=???????????? 解得 0.972579 0.050035a b =?? =? 拟合公式为20.9725790.050035y x =+ 均方误差 2 4 2 2 0[]0.015023i i i y a bx σ==--=∑ 21.给出()ln f x x =的函数表如下： 用拉格朗日插值求ln 0.54的近似值并估计误差（计算取1n =及2n =） 解：1n =时，取010.5,0.6x x == 由拉格朗日插值定理有 1 100.60.5 0.693147 0.510826 0.50.(60.60.51.82321)0 1.()6047()52 j j j x x x L x f x l x ==------=-=∑ 所以1ln0.54(0.54)0.620219L ≈=- 误差为ln 0.54(0.620219)= 0.004032ε=-- 2n =时，取0120.4,0.5,0.6x x x === 由拉格朗日插值定理有
13- 设函数()f x 具有二阶连续导数,{}(*)0,'(*)0,''(*)0,k f x f x f x x =≠≠是由牛顿迭代 法产生的序列,证明 121''(*)lim ()2'(*)k k k k k x x f x x x f x +→∞--=-- 解 牛顿迭代法为 1(),0,1,2,...'()k k k k f x x x k f x +=- = 故 1() '()k k k k f x x x f x +-=- 2112112 121212211()'()()'()()()(*)['()][()(*)]'() '()['()](*)'()['()](*)k k k k k k k k k k k k k k k k k k x x f x f x x x f x f x f x f x f x f x f x f x f f x x x f x f x x ξξ+--------??-=-=??-?? --=---- 其中k ξ介于k x 与*x 之间,1k ξ-介于1k x -与*x 之间,根据式(7.14)得 211222111'()['()]*lim lim () '()['()](*)1''(*)2'(*)k k k k k k k k k k k k x x f f x x x x x f x f x x f x f x ξξ+-→∞→∞-----=-=--- 公式7.14----12''(*)lim 2'(*)k k k e f x e f x +→∞= 15- 考虑下列修正的牛顿公式(单点斯蒂芬森方法) 21()(())()k k k k k k f x x x f x f x f x +=-+- 设()f x 有二阶连续导数,(*)0,'(*)0f x f x =≠,试证明该方法是二阶收敛的. 证明: 将(())k k f x f x +在k x 处作台劳展开,得 21(()) ()'()()''()()2k k k k k k f x f x f x f x f x f f x ξ+=++
桂林电子科技大学 数学与计算科学学院实验报告 实验室： 06406 实验日期： 2014 年 11 月 21 日 院（系） 数学与应用数学 年级、专业、班 姓名 成绩 课程 名称 数值分析实验 实验项目 名 称 实验积分 指导 教师 李光云 一 、实验目的 通过实验掌握利用Matlab 进行数值积分的操作，掌握Matlab 中的几种内置求积分函数，进一步理 解复化梯形，复化辛普生公式，并编程实现求数值积分 二、实验原理 Matlab 中，有内置函数计算积分： >> z = trapz(x,y) 其中，输入x ，y 分别为已知数据的自变量和因变量构成的向量，输出为积分值。 >> z = quad(fun,a,b) 这个命令是使用自适应求积的方法计算积分的命令。其中，fun 为被积函数，a ，b 为积分区间。 我们还可以利用复化梯形公式 ()()()()?? ? ??++-=∑-=11022n i i n x f x f x f n a b I 三、使用仪器，材料 电脑 MATLAB 四、实验内容与步骤 1. 编写复化辛普生公式的Matlab 的程序。 2. 利用复化梯形法程序计算1 2041I dx x =+?，记录下计算结果随着n 增加的变化情况，画图与复化梯形公式的情况比较收敛速度。 3. 积分 ?dx x x sin 的原函数无法用初等函数表达，结合Matlab 复化梯形程序，用描点法绘制其原函数?x dt t t 1sin 在区间[]50,1的图形。 五、实验过程原始记录（数据，图表，计算等）
一、 复化Simpson公式程序： function s=Simpson(a,b,n) %输出s为积分的数值解，输入(a,b)为积分区间，n为等分区间的个数. h=(b-a)/(n*2); s1=0; s2=0; s=h*(f(a)+f(b))/3;%先计算特殊两点相加. for k=1:n x1=a+h*(2*k-1); %利用循环计算其他点的相加. s1=s1+f(x1); end for k=1:(n-1) x2=a+h*2*k; %利用循环计算其他点的相加. s2=s2+f(x2); end s=s+h*(4*s1+2*s2)/3; 画图程序 format long; % k为等分区间个数，t存储积分值. k=2:1:40; for i=1:length(k) t(i)=Simpson(0,1,k(i)); disp([k(i),t(i)]); end plot(k,t,'.','MarkerSize',20) 二、 Simpson(0,1,26) ans = 3.14159265358779 >> Untitled6 2.00000000000000 3.14156862745098 3.00000000000000 3.14159178093604 4.00000000000000 3.14159250245871 5.00000000000000 3.14159261393922 6.00000000000000 3.14159264030538 7.00000000000000 3.14159264832065 8.00000000000000 3.14159265122482
数值分析第六章整合版（黑组） 一、填空题 1、已知 ()01 P x =，()1P x x =， () () 2 2312 x P x -= ，根据勒让德多项式的递推关系，则 求()3P x =（3532x x - ） 解：勒让德多项式的递推关系为()()()()()11121n n n n P x n xP x nP x +-+=+-，n=1，2……. 将()1P x x =，()() 2 2312 x P x -= 代入上式即可求出()3P x =3532 x x - 2、若)(x P 是],[)(b a C x f ∈的最佳3次逼近多项式，则)(x P 在],[b a 上存在5 个交替为 正、负偏差点。（考点：切比雪夫定理） 3、切比雪夫正交多项式可表示为(x)cos(narcosx)n T =，(x)n T 是最高次幂系数为12n - 的n 次多项式。（考点：切比雪夫多项式性质） 4、最佳一致问题同时存在正偏差点和负偏差点 （考点：最佳一致逼近定理3） 二、选择题 1、求函数3)1()(+=x x f 在区间[0,1]，],[,21b a x x ∈上的一次最佳一致逼近多项式（D ） A x +4358.0 B x 34358.0+ C x 54358.0+ D x 74358.0+ 2、设 的2-其中 为定义在[a,b]上的（A ） A 权函数 B 反函数 C 幂函数 D 函数 3、x e =）（x f ，-1≤x ≤1，且设= p(x)x a a 1 +，求a a 1 , 0使得)(x p 为)(x f 于[] 1,0上的最佳平方逼近多项式（A ） A:() 1021--=e e a ,311e a -= B:() e a e a e 111 03 1,2---== ) (x ρ],[)(b a C x f ∈()f x
思考题 10.设m T 是m 2+1个节点的复化梯形值，试证：143m m m T T S +-=就是m+12+1个节点的复化Simpson 值。 解：由复化梯形公式： 11()[()()2()]2n b a k h f x dx f a f b f a bh -=≈+++∑?，其中b a h n -=。 由题中条件可知：节点为m 2+1时 11 21112(2)[()()2()]2[()()2()]2(2)m m m m n k k T b a h f a f b f a bh b a f a f b f a b -==++=-+++=-+++?∑∑ 节点为m+12+1时： 11111121112(2)[()()2()]2[()()2()]2(2)m m m m n k k T b a h f a f b f a bh b a f a f b f a b ++++-==++=-+++=-+++?∑∑ m+12+1的复化Simpson 公式为： 21112111()[()()4()2()]3(2)i m m m b i a k k b a f x dx f a f b f x f x -+-==+-≈+++∑∑? 12111111211 2121 411433 11312(2)2(2)[()()4()2()]3(2)[()()2()]2(2)[()()2()]2(2)m m m i m m m m m m m m m i k k k k T T S b a b a b a f a f b f x f x b a f a f b f a b b a f a f b f a b +-++++-====++-==++-=+---+++-+++?-+++?∑∑∑∑
第6章 数值积分 --------学习小结 一、 本章学习体会 本章主要介绍了五种计算定积分的数值积分法，分别为：插值型求积公式、Newton-Cotes 求积公式、复化梯形公式与复化Simpson 公式、Gauss 型求积公式等。本章的重点在于掌握求积公式及其运用，并要学会求代数精度。而通过对求积公式进行比较，会发现其方法与以前所学习的解析方法有一定的不同，它并不需要求出定积分的原函数，而是去直接利用求积公式来求出所给定积分的近似值，使其达到一定的求解精度要求，从而根据不同的题型做出不同的解答，这对于我们今后的专业研究过程也有一定的作用。 例如：高阶Newton-Cotes 公式会出现数值不稳定，而低阶Newton-Cotes 公式有时又不能满足精度要求，可将积分区间[a ，b]分成若干小区间，在每个小区间上用低阶求积公式计算，然后求和，即运用复化求积法。通过运用matlab 软件，可以加深自己对各种求积公式的理解。根据求解要求，充分考虑已知条件，选择简便快捷的求积方法进行定积分求解，从而得出比较准确的结果。通过查阅相关书籍，加深对课本知识的理解，从而提高自己的自学能力。 二、本章知识梳理 1 求积公式及其代数精度： 求积公式的一般形式： () 0()()n b n k k a k f x dx f x λ=≈∑? 截断误差或余项：)()(0 k b a n k k n x f dx x f R ?∑=-=λ 代数精度：对于上面所列的求积公式，当()f x 为任何次数不高于m 的多项式时都成为等式，而当()f x 为某个m+1次多项式时不能成为等式，则称它具有m 次代数精度。 2 插值型求积公式： () ()()n b n k k a k f x dx f x λ =≈∑? 其中 ()()(0,1,...,) b n k k a l x dx k n λ ==?