第7章 数值积分法
7.1 实验目的
了解求积公式及代数精度概念，理解并掌握求定积分的求积公式的算法构造和计算，学习用计算机求定积分的一些科学计算方法和简单的编程技术和能用程序实现这些算法。 7.2 概念与结论
1. 求积公式
计算定积分的如下形式的近似公式：
称为求积公式。
2．代数精度 若求积公式
对一切不高于m 次的 多项都准确成立，而对于m+1次多项式等号不成立，则称此求积公式 的代数精度为m 。
代数精度越高，求积公式越好。
3.求积余项
4.Newton-Cotes 求积公式的代数精度
n 点Newton-Cotes 求积公式的代数精度至少可以达到n-1,且当n 为奇数时，可以达到n 。
?∑=≈b
a
n
k k k x f A dx x f 1
)
()(?∑=≈b
a
n
k k k x f A dx x f 1
)
()(?∑=-=b
a
n
k k k x f A dx x f f R 1
)
()()(
?
=b
a
dx
x f I )(5.Richardson 外推定理
设函数F 1(h)逼近量F*的余项为：
F*-F 1(h)=a 1h p1
+a 2h p2
+····+a k p k
+···
式中p k >p k-1>···>p 2>p 1>0, F*和a i (i=1,2, ···)都是与h 无关的常数,且k ≥1时，a k ≠0，则由：
定义的函数F 2(h)也逼近F*，且有
F*-F 2(h)= b 2h p2
+····+b k p k
+···
6. 关于复合梯形公式的展开定理
设f(x)在[a,b]区间上无穷次可微，则有如下展开式：
T(h)=I+a 1h 2
+a 2h 4
+a 3h 6
+…+a m h 2m
+…
式中T(h)是函数f(x)在[a,b]区间上的复化梯形值Tn,
7.3 程序中Mathematica 语句解释 1. 随机函数
Random[] 随机给出闭区间[0,1]内的一个实数 Random[Real, xmax] 随机给出闭区间[0,xmax]内的一个实数 Random[Real, {xmin, xmax}] 随机给出闭区间[xmin,xmax]内的一个实数 Random[Integer] 随机给出整数0或1
Random[Integer, {xmin, xmax}] 随机给出xmin 到xmax 之间的一个整数 Random[Complex] 随机给出单位正方形内的一个复数 2．{a1,a2,…,an}
表示由元素a1,a2,…,an 组成的一个表，元素可以是任何内容。
)
10(1)
()()(1
1112<<--=
q q
h F q qh F h F p p
如：{1，3，4，5}，{1,x,{2,3},x+y}，{{1，3}，{1，2，3}，{3，2，4}}等 3．list[[k]] 表list 中的第k 个元素 4．list[[i,j]]
表list 中第i 个元素中的第j 个元素,此时list 中的第i 个元素应该也是一个表。 7.4 方法、程序、实验
在实际问题中，往往会遇到被积函数f(x)的原函数无法用初等函数来表示，或函数只能用表格表示，或有的虽然能用初等函数表示，但过分复杂，所以这些情形都需要去建立定积分的近似计算公式来做积分计算。数值积分是进行定积分计算的一种方法，它可以解决不能用定积分基本公式计算的所有定积分问题。数值积分涉及很多计算公式，这里主要介绍Newton-Cotes 求积公式、复合求积公式、Romberg 求积方法和Monte-Carlo 方法的构造过程和算法程序。
1. n 点 Newton —Cotes 求积公式
n 点 Newton —Cotes 求积公式又称为等距节点求积公式，它是利用被积函数f(x)在积分区间[a,b]的n 个等分节点上的函数值构造的插值函数?(x)代替f(x)做定积分计算所构造求积公式。这个求积公式是通常做定积分近似计算的梯形公式和抛物线公式的推广，主要在理论上用的多些。 1) n 点 Newton —Cotes 求积公式的构造过程
将积分区间[a,b] 分为n-1等分，其中n 个节点 x i =a+(i-1}h, i=1,2,…,n ，h=(b -a)/(n-1)，然后用f(x)在这n 个节点上建立插值于f(x)的n-1次代数多项式P n-1(x)，引入变换
x=a+th, 0≤t ≤n-1
则有
)
1
)(())(()(,11,111∏∑∏∑≠==≠==--+-=--=n i
k k n i i n
i k k k i k n
i i n k i k t x f x x x x x f x P
带入定积分，有：
C k (n)
称为Cotes(柯特斯)系数, 则得到n 点 Newton —Cotes 求积公式:
n 点 Newton —Cotes 求积公式的求积余项为
当n=2时，2点的 Newton —Cotes 求积公式就是如下梯形公式：
梯形求积公式求积余项为
当n=3时，3点的Newton —Cotes 求积公式就是如下抛物线（Simpson ）公式：
?∑=-≈b
a
n
k k n k
x f C a b dx x f 1
)
()
()()())()((2
)(b f a f a
b dx x f b
a
+-≈
?
?∏
?∏
∑?
∏
∑??
-≠=-≠==≠==--+--=
-+---=--=≈1
0,1)
(1
0,11
,11
11
1
1)1
)((1)
)(()()(n n
i
k k n i n n
i
k k n i i b
a n
i
k k k
i k
n
i i b
a
n b a
dt k
i k t n c dt k
i k t x f n a b dx x x x x x f dx x P dx x f 令
dx
x x x x x x n f f R n b a n )())((!
)
()(21)(---=?
ξ],{)()(12
1
)(3b a f a b f R ∈''--
=ηη
Simpson 求积公式求积余项为
如果想得到其他的 Newton —Cotes 求积公式只要在有关书中查出Cotes 系数表就可以马上得到相应的Newton —Cotes 求积公式。
2) n 点 Newton —Cotes 求积公式算法
1． 输入被积函数f(x)及积分上下限a,b 2． 选择Cotes 系数构造求积公式 3． 用求积公式求定积分
3) n 点 Newton —Cotes 求积公式程序
Clear[a,b,x,n,s]; a=Input["a="] b=Input["b="]
f[x_]=Input["被积函数f(x)="] n= Input["求积节点个数n="]; c={{1/2,1/2}, {1/6,4/6,1/6}, {1/8,3/8,3/8,1/8},
{7/90,16/45,2/15,16/45,7/90},
{19/288,25/96,25/144,25/144,25/96,19/288}, {41/840,9/35,9/280,34/105,9/280,9/35,41/840},
))()2
(4)((6)(b f b a f a f a b dx x f b a +++-≈?
],{)()(2880
1
)()4(5b a f a b f R ∈--
=ηη
{751/17280,3577/17280,1323/17280,2989/17280,1323/17280,3577/17280,751/17280},
{989/28350,5888/28350,-928/28350,10496/28350,-4540/28350,
10496/28350,-928/28350,5888/28350,989/28350}};
h=(b-a)/(n-1);
x=Table[a+k*h,{k,0,n-1}];
s=(b-a)*Sum[c[[n-1,k]]*f[x[[k]]],{k,1,n}]
Print["定积分=",N[s,8]];
说明本程序用n（n=2,3,4,5,6，7，8，9）点 Newton—Cotes求积公式求[a,b]上的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和求积节点个数n后，计算机则给出定积分的近似值。
程序中变量说明
a：存放积分下限
b: 存放积分上限
f[x]: 存放被积函数f(x)
n: 存放求积节点个数
c: 存放Cotes系数
s: 存放定积分近似值
h: 存放节点步长
x：存放节点xi
注语句
c={{1/2,1/2},{1/6,4/6,1/6},{1/8,3/8,3/8,1/8},{7/90.16/45,2/15,16/45,7/90},
{19/288,25/96,25/144,25/144,25/96,19/288}}
的第i个分量表是具有i个节点的Cotes系数。
4)例题与实验
例1.用n=3和n=4的Newton-Cotes求积公式求定积分
的近似值。
解：执行n 点 Newton —Cotes 求积公式程后，在输入的窗口中按提示分别输入
1、3、Exp[-x/2]、3
每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出用n=3的Newton-Cotes 求积公式计算出的定积分结果：
)613261(
23
E
E E ++
定积分=0.76705953
再执行n 点 Newton —Cotes 求积公式程后，在输入的窗口中按提示分别输入
1、3、Exp[-x/2]、4
每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出用n=4的Newton-Cotes 求积公式计算出的定积分结果：
)
81838381(
26
5
6
7
3
E
E
E
E
+
+
+
定积分=0.76691628
因此用n=3和n=4的Newton-Cotes 求积公式求本题定积分近似值分别为
0.76705953和0.76691628
注意到本题的精确值为0.766800999….,可见n=4的Newton-Cotes 求积公式计算结果较好。
2. 复化求积公式
复化求积公式是把积分区间分成若干小区间，在每个小区间上采用次数不高的插值公式，如梯形公式或抛物线公式，构造出相应的求积公式，然后再把它们加起来得到整个区间上的求积公式。复化求积公式克服了高次Newton-Cotes 公式计算不稳定的问题，其运算简单且易于在计算机上实现。常用的复化求积公式是复化梯形公式和复化抛物线公式，下面分别讨论。
dx e x ?
-31
2
1) 复化梯形公式的构造过程
把区间[a,b] n 等分，取节点x k =a+kh,k=0,1,...n, h=(b-a)/n,对每个小区间[x k ,x k+1]用梯形求积公式，再累加起来得：
公式
就是复化梯形公式。它的求积余项为：
由复化梯形公式的余项，可以看到，当n 不断变大时， T n 无限接近定积分值。因此复化梯形公式可以使定积分的计算达到任意精度。为了计算简单，提高效率，常用|T 2n – T n |<ε来得到满足精度要求的定积分值。 2) 复化抛物线公式的构造过程 在每个小区间[x k ,x k+1]上，有
把区间[a,b] n 等分，取节点x k =a+kh,k=0,1,...n, h=(b-a)/n,对每个小区间[x k ,x k+1] 用抛物线求积公式，再累加起来得：
∑?
-==++-≈11
))(2)()((2)(n k n
k b a T x f b f a f n a b dx x f ∑∑?
∑?
-=+-=-=++-=+≈=+1
1
11
010))
(2)()((2))()((2)()(1
n k k k k n k b a n k x x x f b f a f n a b x f x f h dx x f dx x f k k
]
,[)(12
)()
(12)(),(2
1
10
3
b a f h a b x x f h T dx x f T f R k k n k k b a n n ∈''--
=≤≤''-=-=+-=∑
?
ηηηη
公式
就是复化抛物线公式。它的求积余项为：
由复化抛物线公式的余项，可以看到，当n 不断变大时，S n 无限接近定积分值。因此复化抛物线公式也可以使定积分的计算达到任意精度，且收敛的速度比复化梯形公式更快。为了计算简单，提高效率，常用| S 2n –S n |<ε来得到满足精度要求的定积分值。 3) 复化梯形求积公式算法
1.输入被积函数f(x),积分上下限a,b,和求积精度ε
2. n ?1, 计算T n
3. 计算T 2n
4． 判断|T 2n – T n |<ε是否成立，如果成立，输出定积分近似值，停止 5． 否则， T n ? T 2n , n ?2n,转3 4) 复化抛物线求积公式算法
1.输入被积函数f(x),积分上下限a,b,和求积精度ε
2. n ?1, 计算S n
2
,))(4)(2)()((6))
()(4)((6)()(2
11
11
2
11211
101
h
x x
x f x f b f a f n a b x f x f x f h
dx x f dx x f k k n k n k k k k k k n k b a n k x x k k
+=+++-=++≈
=+
-=-=+++-=-=∑∑∑?
∑?
+]
,[)(2
180)()
(2
180)(),()4(4
41
10)4(4
5
b a f h a b x x f h S dx x f T f R k k n k k b a n n ∈?--=≤≤?-=-=+-=∑
?
ηηηηn n k n k k k b a S x f x f b f a f n a
b dx x f =+++-≈∑∑
?
-=-=+
11102
1))(4)(2)()((6)(
3.计算S2n
4.判断|S2n– S n|<ε是否成立，如果成立，输出定积分近似值，停止
5.否则， S n? S2n , n?2n,转3
5) 复化梯形求积公式程序
(*复合梯形公式*)
Clear[a,b,x,n,f];
a=Input["a="]
b=Input["b="]
f[x_]=Input["被积函数f(x)="]
eps=Input["精度要求eps="];
n=1;
h=b-a;
t1=(f[a]+f[b])*h/2;
h=h/2;
t2=t1/2+h*f[a+h];
er=t2-t1//N;
While[Abs[er]>eps,
Print["n=",2^n,"定积分值为",N[t2,10]];
Print["误差=",er];
h=h/2;
t1=t2;
n=n+1;
t2=t1/2+h*Sum[f[a+k*h],{k,1,2^n,2}];
er=t2-t1//N];
Print["n=",2^n,"定积分值为",N[t2,10]];
Print["误差=",er]
说明 本程序从梯形公式T 1开始，用复合梯形求积公式求[a,b]上满足精度小于eps 要求的定积分近似值。程序执行后，按要求通过键盘输入积分下限a 、积分上限b 、被积函数f(x)和精度要求eps 后，计算机则给出满足精度要求的定积分近似值及中间计算值和误差。 程序中变量说明 a ：存放积分下限 b: 存放积分上限 f[x]: 存放被积函数f(x) eps: 存放求积精度ε h: 存放节点步长 x ：为函数f[x]:提供变量 t1: 存放复合梯形值T n
t2: 存放复合梯形值T 2n 和定积分近似值 n: 存放复合梯形公式的节点加密次数 er:存放误差
注：为提高计算效率，计算T 2n 时使用了公式
6) 复化抛物线求积公式程序 (*复合抛物线公式*） Clear[a,b,x,n,f,s]; a=Input["a="] b=Input["b="]
f[x_]=Input["被积函数f(x)="] n=Input["分割区间数n="]; h=(b-a)/n; a1=a+h/2;
s=h/6*(f[a]+f[b]+
2Sum[f[a+k*h],{k,1,n-1}]+4Sum[f[a1+k*h],{k,0,n-1}]);
∑=??
? ??--+-+=n i n n i n a b a f n a b T T 12)12(2221
Print["n=",n];
Print["定积分值为",N[s,10]]
说明 本程序用复合抛物线求积公式求[a,b]上定积分近似值。程序执行后，按要求通过键盘输入积分下限a 、积分上限b 、被积函数f(x)和分割区间数n 后，计算机则给出定积分的近似值Sn 。 程序中变量说明 a ：存放积分下限 b: 存放积分上限 f[x]: 存放被积函数f(x) h: 存放节点步长 x ：为函数f[x]提供变量 n: 存放分割区间数
s: 存放复合抛物线值S n 和定积分近似值
注：Mathematica 数学软件中也有一个求数值积分的命令，命令形式为：
NIntegrate[f[x],{x,a,b}]
它可以求出[a,b}上的定积分近似值.
7) 例题与实验
例2. 用复合梯形公式求定积分
的近似值,要求误差小于10-7
。
解：执行复化梯形求积公式程后，在输入的窗口中按提示分别输入：
0、3、Exp[-x^2]、10^(-7)
每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出的计算结果如下： n=2 定积分值为0.7313702518 误差=0.0474305
dx e x ?
-10
2
n=4 定积分值为0.7429840978 误差=0.0116138 n=8 定积分值为0.7458656148 误差=0.00288152 n=16 定积分值为0.7465845968 误差=0.000718982 n=32 定积分值为0.7467642547 误差=0.000179658 n=64 定积分值为0.7468091636 误差=0.000044909 n=128 定积分值为0.7468203905 误差=0.0000112269 n=256 定积分值为0.7468231972 误差=2.8067?10
-6
n=512 定积分值为0.7468238989 误差=7.01676?10-7
n=1024 定积分值为0.7468240743 误差=1.75419?10
-7
n=2048 定积分值为0.7468241182 误差=4.38548?10-8
从计算结果可以得出经过11次加密分割积分区间得到满足精度要求的定积分近似值0.7468241182。此外，通过计算过程可以明显看到复合梯形求积公式的收敛性。
例3. 分别用n=1,2,4,8,16,32的复合抛物线求积公式求定积分
的近似值。
解：执行复化抛物线求积公式程后，在输入的窗口中按提示分别输入：
0、 Pi 、Exp[x]*Cos[x]、1
每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出的计算结果如下： n=1
定积分值为-11.59283955
再重复如上操作，将分割区间数分别输入2，4，8，16，32，得如下输出结果： n=2
定积分值为-11.98494402 n=4
定积分值为-12.06420896
dx
x e x ?π
cos
n=8
定积分值为-12.06995132 n=16
定积分值为-12.07032146 n=32
定积分值为-12.07034476。
注意到本题的定积分值为-12.07034631639，从中可以看到复合抛物线求积公式收敛还是很快的。你能通过修改复合抛物线求积公式一次得出如上所有计算结果吗？
3.Romberg 求积公式
Romberg 又称逐次分半加速收敛法，它是对复合梯形求积公式采用Richardson 加速技术得到的一种数值求积方法。
1) Romberg 求积公式的构造过程
由复合梯形公式的展开定理，有如下关系式：
T 1(h)-I=a 1h 2
+a 2h 4
+a 3h 6
+…+a m h 2m
+…
这里
利用Richardson 外推定理对T 1(h)进行加速，注意到这里p 1=2,取q=1/2有
于是得到收敛于I 的加速公式T 2(h),如果再T 2(h) 进行加速，可以有
继续加速下去，可以得到加速序列
Tn
h T dx x f I b a
==?)(,)(11
4)
()2(4211)()21
()2()(1121212--=
??
? ??--=
h T h
T h T h T h T 1
4)
()2(4)(22223--=
h T h
T h T
式中m 代表对I 进行的加速次数。
此外，注意到T 1(h/2)=T 2n ，且直接计算可以知道T 2(h)就是复合抛物线公式S n ,于是有计算关系式 类似地，还有如下计算格式：
和
于是我们至少可以得到如下3个收敛于定积分I 的数列： T 型值数列：T 1、T 2、T 4、T 8 、T 16 、…T 2n
、… Simpson 数列：S 1、S 2、S 4、S 8 、S 16 、…S 2n 、… Cotes 数列：C 1、C 2、C 4、C 8 、C 16 、…C 2n
、… Romberg 数列：R 1、R 2、R 4、R 8 、R 16 、…R 2n 、…
利用Romberg 数列求定积分I 的方法称为Romberg 求积方法。
如果给定求积精度ε，可以用
作为终止计算的条件，并取最后算出的R 值作为积分值。
2) Romberg 求积算法
1.输入被积函数f(x),积分上下限a,b,和求积精度ε
2.按顺序计算T 1、T 2、T 4、T 8 并做 t1? T 1、t2? T 2、t3? T 4、t4? T 8、
1
442--=
Tn T S n n 1
4)
()2(4)(11--=
--m m m m m h T h
T h T 1
44222--=
Sn S C n n 1
44323--=
Cn C R n n ε
≤--R
R n n
1
22
s1? (4t2-t1)/3,s2? (4t3-t2)/3,s3? (4t4-t3)/3;
c1? (16s2-s1)/15,c2? (16s3-s2)/15,R2? (64c2-c1)/63;
3.n?1,R1? c2 , t1?t4,s1?s3,c1?c2
4.判断|R2 – R1|<ε是否成立，如果成立，输出定积分近似值R2，停止
5.否则， R1? R2 ,
t2? T16n,
s2? (4t2-t1)/3, t1?t2,
c2? (16s2-s1)/15, s1?s2,
R2? (64c2-c1)/63, c1?c2
n?2n,转4
3) Romberg求积方法程序
Clear[a,b,x,n,f];
a=Input["a="]
b=Input["b="]
f[x_]=Input["被积函数f(x)="]
eps=Input["精度要求eps="];
m=1;h=b-a;
t=Table[0,{4}];
t[[1]]=(f[a]+f[b])*h/2;
Do[m=2m;
h=h/2;
t[[k+1]]=t[[k]]/2+h*Sum[f[a+i*h],{i,1,m,2}]//N
,{k,1,n-1}]
s1=(4t[[2]]-t[[1]])/3;
s2=(4t[[3]]-t[[2]])/3;
s3=(4t[[4]]-t[[3]])/3;
c1=(16s2-s1)/15;
c2=(16s3-s2)/15;
r1=c2;
r2=(64c2-c1)/63;
t1=t[[4]];s1=s3;c1=c2;er=r2-r1;
k=1;
While[Abs[er]>eps,
r1=r2;
Print["r(",k,")=",N[r2,20]," eps=",N[er,10]];
h=h/2;m=2m;
t2=t1/2+h*Sum[f[a+i*h],{i,1,m,2}]//N;
s2=(4t2-t1)/3;
c2=(16s2-s1)/15;
r2=(64c2-c1)/63;
t1=t2;
s1=s2;
er=r2-r1;
k=k+1;
c1=c2];
Print["r(",k,")=",N[r2,20]," eps=",N[er,10]];
说明本程序从用Romberg求积方法求[a,b]上满足精度小于eps要求的定积分近似值。程序执行后，按要求通过键盘输入积分下限a、积分上限b、被积函数f(x)和精度要求eps 后，计算机则给出满足精度要求的定积分近似值及中间计算值和误差。
程序中变量说明
a：存放积分下限
b: 存放积分上限
f[x]: 存放被积函数f(x)
eps: 存放求积精度
h: 存放节点步长 x ：为函数f[x]:提供变量 t,t1 ,t2: 存放复合梯形值T n s1 ,s2: 存放Simpson 数列值S n c1 ,c2: 存放Cotes 数列值C n R1 ,R2: 存放Romberg 数列值R n m: 存放复合梯形公式的节点加密次数 er:存放误差
注：为编程方便，程序中先算出了产生Romberg 数列的值：T 1、T 2、T 4、T 8 、S 1、S 2、S 4、C 1、C 2使，其中t=Table[0,{4}]提供一个4元素的表用于存放T 1、T 2、T 4、T 8。 4) 例题与实验
例4.用Romberg 求积方法求定积分
的近似值,要求误差小于10-12
。
解：执行Romberg 求积方法程序后，在输入的窗口中按提示分别输入：
0、1、 4/(1+x^2）、10^(-12)
每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出的计算结果如下： r(1)=3.141585783761874 eps=-8.310364015?10
-6
r(2)=3.141592638396796 eps=6.854634922?10-6
r(3)=3.141592653590029 eps=1.519323334?10-8
r(4)=3.141592653589794 eps=-2.358113704?10-13
dx
x ?+10214
因此得到满足精度要求的定积分值为3.141592653589794 ，误差=-2.358113704?10-13
4. Monte-Carlo 求积方法
Monte-Carlo 求积方法以概率论大数定理为依据式，借助随机数来求定积分的近似值的一种方法，该方法特别适用于求多重积分和积分区域复杂的重积分，不足之处是收敛较慢。 1) Monte-Carlo 求积方法的构造过程
假设重积分
是有限的，式中Ω?R n
，G(x)是定义在Ω上的多元函数, P(x)是定义在Ω上的分布函数，则有
只要式中x (k)
是Ω上关于分布P(x)的n 个任意独立选取的随机点列即可，这个结论可以有概率论中的大数定律：
以概率1成立。
利用以上结果求重积分的方法就称为Monte-Carlo 求积方法。对具体的积分计算有如下计算公式：
1） 求定积分的计算公式
式中x k 是[a,b]均匀分布的n 个任意独立选取的随机数列。
∑=≈n
k k x G n G I 1
)()
(1)(x
d x P x G G I ?Ω
=)()()()()(11
)(G I x G n Lim n
k k n =∑=∞→∑
?
=-≈n k k b a x f n a b dx x f 1
)()(
2） 求二重积分的计算公式
式中(x k,y k )是 平面区域D 中均匀分布的n 个任意独立选取的随机点列，|D|表示区域D 的面积。
3） 求三重积分的计算公式
式中(x k,y k ,z k )是空间区域Ω中均匀分布的n 个任意独立选取的随机点列，|Ω|表示区域Ω的体积。 Monte-Carlo 求积方法的收敛阶为O(n 1/2
). 2) Monte-Carlo 求积方法算法
1.输入被积函数f(x),积分区域边界和随机点的个数n
2.在积分区域中产生n 个随机点列
3.计算被积函数f(x)在随机点列上的函数值并相加求平均
4.用平均值与积分区域的测度之积作为积分近似值，这里测度就是长度、面积、体积之类的度量。 3) Monte-Carlo 求积方法程序 （*求定积分*） Clear[a,b,x,n,f]; a=Input["a="] b=Input["b="]
f[x_]=Input["被积函数f(x)="] n=Input["随机点个数n="];
s=(b-a)*Sum[f[Random[Real,{a,b}]],{n}]/n Print["n=",n," 积分≈",N[s,10]];
说明 本程序从用Monte-Carlo 求积方法求[a,b]上的定积分近似值。程序执行后，按要求通过键盘输入积分下限a 、积分上限b 、被积函数f(x)和随机点个数n 后，计算机则给出定积分近似值。 程序中变量说明
∑??
=≈n
k k k D
y x f n D d y x f 1
)
,(||),(σ∑???
=Ω
Ω≈n
k k k k z y x f n dv z y x f 1
)
,,(||),,(
东北大学秦皇岛分校 数值计算课程设计报告 数值积分算法及MATLAB实现 学院数学与统计学院 专业信息与计算科学 学号5133201 姓名陈悦 指导教师姜玉山张建波 成绩 教师评语： 指导教师签字： 2015年07月14日
1 绪论 数值分析是计算数学的一个主要部分，计算数学是数学科学的一个分支，它研究用计算机求解各种数学问题的数值检索方其理论与软件的实现.而数值分析主要研究数值计算. 现科学技术的发展与进步提出了越来越多的复杂的数值计算问题，这些问题的圆满解决已远人工手算所能胜任，必须依靠电子计算机快速准确的数据处理能力.这种用计算机处理数值问题的方法，成为科学计算.今天，科学计算的应用范围非常广泛，天气预报、工程设计、流体计算、经济规划和预测以及国防尖端的一些科研项目，如核武器的研制、导弹和火箭的发射等，始终是科学计算最为活跃的领域. 1.1 数值积分介绍 数值积分是数值分析的重要环节，实际问题当中常常需要计算积分，有些数值方法，如微分方程和积分方程的求解，也都和积分计算相联系. 求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来，因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的.另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解.由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题.对微积分学做出杰出贡献的数学大师，如I.牛顿、L.欧拉、C.F.高斯、拉格朗日等人都在数值积分这个领域作出了各自的贡献，并奠定了这个分支的理论基础. 构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式.特别在节点分布等距的情形称为牛顿-科特斯公式，例如梯形公式(Trapezoidal Approximations)与抛物线公式(Approximations Using Parabolas)就是最基本的近似公式.但它们的精度较差.龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式(Rhomberg Integration).当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分.数值积分还是微分方程数值解法的重要依据.许多重要公式都可以用数值积分方程导出.现探讨数值积分算法以及运用MATLAB软件的具体实现
实验１面向方程的数值积分方法仿真 1. 实验目的 运用CSS01.C 仿真程序解题,培养阅读及修改仿真程序的能力,学习并了解仿真程序的结构及特点.通过实验,加深理解4阶龙格-库塔法的原理及其稳定域. 2. 实验设备：装有BC 语言的PC 机一台 3. 实验内容 修改CSS01.C 仿真程序,对如下系统进行仿真. (1) 线性定常系统 ??????????321x x x =????????? ?---00600120000120??????????321x x x +??????????60000u y=[]00 1[]T x x x 32 1 ??????????)0()0()0(321x x x =???? ? ?????000, u=1(t) a)、对龙格库塔法进行分析：它是一种数值积分法，也就是微分方程初 值问题数值计算法，是对初值微分方程的离散化求解。对于数值积分法我们常用的是欧拉法以及二阶和四阶龙格库塔法其原理分别如下： 对于形如 dt dy =),(y t f 的微分方程 用欧拉法仿真其迭代公式为 ),(*1k k k k y t f h y y +=+，其中h 为仿真步长（下同）；
二阶龙格库塔法其迭代公式为 y k+1= y k +h/2*(k1+k2)，其中),(1k k y t f k =， );*,(21h k y h t f k k k ++= 四阶龙格库塔法其迭代公式为：)432221(*6/1k k k k h y y k k ++++=+， 其中),(1k k y t f k =，k2= f(t k +h/2,y k +k1*h/2)，k3= f(t k +h/2,y k +k2*h/2)， k4= f(t k +h,y k +k3*h)； 本程序中大致分以下四次计算： 第一次计算： g[1]*h/2——k[1][2]，y[1]——k[1][1]，把y[1]的值赋给k[1][1] Y[1]= k[1][1]+ k[1][2]； dt dy =),(y t f ——g[1] ， （计算k1）， g[1]*h/2——k[1][2]，（把k1*h/2存储到k[1][2]） 第二次计算： Y[1]= k[1][1]+ k[1][3]； dt dy =),(y t f ——g[1] ， （计算k2）， g[1]*h/2——k[1][3]，（把k2*h/2存储到k[1][3]） 第三次计算： dt dy =),(y t f ——g[1]，（计算k3）， Y[1]= k[1][1]+ k[1][4]； k[1][4]= g[1]*h ；（把k3*h 存储到k[1][4]） 第四次计算： dt dy =),(y t f ——g[1]，（计算k4）， 此时 把k1,k2,k3,k4 的值代入下式： ) 432221(*6/1k k k k h y y k k ++++=+ 即为 6/)4*3**22**21*(1k h k h k h k h y y k k ++++=+ ， 由于h*k1=2*k[1][2]，2*h*k2=4*k[1][3]，2*h*k3=2*k[1][4]，h*k4=h*g[1] 所以可得如程序中的计算公式： Y[1]= k[1][1]+(2*k[1][2]+4* k[1][3]+2* k[1][4]+h*g[1]）； 可见该仿真程序为四阶龙格库塔法的数值积分仿真程序。
学科分类号110.3420 州 GUIZHOU NORMAL COLLEGE 本科毕业论文 题目—几种常用数值积分方法的比较_____________ 姓名潘晓祥学号1006020540200 院（系）数学与计算机科学学院 __________________ 专业数学与应用数学年级_____________2010级 指导教师雍进军职称______________________讲师 二O—四年五月
贵州师范学院本科毕业论文（设计）诚信声明本人郑重声明：所呈交的本科毕业论文（设计），是本人在指导老师的指导下，独立进行研究工作所取得的成果，成果不存在知识产权争议，除文中已经注明引用的内容外，本论文不含任何其他个人或集体已经发表或撰写过的作品成果。对本文的研究做出重要贡献的个人和集体均已在文中以明确方式标明。本人完全意识到本声明的法律结果由本人承担。 本科毕业论文作者签名： 年月曰
贵州师范学院本科毕业论文（设计）任务书
研究方法: 本论文主要通过对相关文献和书籍的参考，合自己的见解，复化求积公式，Newton —Cotes求积公式,Romberg求积公式，高斯型求积公式进行讨论并进行上机实验，从代数精度，求积公式误差等角度对这些方法进行分析比较完成期限和采取的主要措施： 本论文计划用6个月的时间完成，阶段的任务如下： (1) 7月份查阅相关书籍和文献； (2) 8月份完成开题报告并交老师批阅； (3) 9月份完成论文初稿并交老师批阅； (4) 10月份完成论文二搞并交老师批阅； (5) 11月份完成论文三搞； (6) 12月份定稿. 主要措施：考相关书籍和文献，合自己的见解，老师的指导下和同学的帮助下完成 主要参考文献及资料名称： [1] 关治?陆金甫?数学分析基础(第二版) [M].北京：等教育出版社.2010.7 [2] 胡祖炽.林源渠.数值分析[M]北京：等教育出版社.1986.3 [3] 薛毅.数学分析与实验[M] 北京：业大学出版社2005.3 [4] 徐士良.数值分析与算法[M].北京：械工业出版社2007.1 [5] 王开荣.杨大地.应用数值分析[M]北京:等教育出版社2010.7 [6] 杨一都.数值计算方法[M].北京:等教育出版社.2008.4 [7] 韩明.王家宝.李林.数学实验(MATLAB版[M].上海：济大学出版社2012.1 [8] 圣宝建.关于数值积分若干问题的研究[J].南京信息工程大学.2009.05.01. : 42 [9] 刘绪军.几种求积公式计算精确度的比较[J].南京职业技术学院.2009. [10] 史万明.吴裕树.孙新.数值分析[M].北京理工大学出版社.2010.4. 指导教师意见: 签名： 年月日
第7章 数值积分法 7.1 实验目的 了解求积公式及代数精度概念，理解并掌握求定积分的求积公式的算法构造和计算，学习用计算机求定积分的一些科学计算方法和简单的编程技术和能用程序实现这些算法。 7.2 概念与结论 1. 求积公式 计算定积分的如下形式的近似公式： 称为求积公式。 2．代数精度 若求积公式 对一切不高于m 次的 多项都准确成立，而对于m+1次多项式等号不成立，则称此求积公式 的代数精度为m 。 代数精度越高，求积公式越好。 3.求积余项 4.Newton-Cotes 求积公式的代数精度 n 点Newton-Cotes 求积公式的代数精度至少可以达到n-1,且当n 为奇数时，可以达到n 。 ?∑=≈b a n k k k x f A dx x f 1 ) ()(?∑=≈b a n k k k x f A dx x f 1 ) ()(?∑=-=b a n k k k x f A dx x f f R 1 ) ()()(
? =b a dx x f I )(5.Richardson 外推定理 设函数F 1(h)逼近量F*的余项为： F*-F 1(h)=a 1h p1 +a 2h p2 +····+a k p k +··· 式中p k >p k-1>···>p 2>p 1>0, F*和a i (i=1,2, ···)都是与h 无关的常数,且k ≥1时，a k ≠0，则由： 定义的函数F 2(h)也逼近F*，且有 F*-F 2(h)= b 2h p2 +····+b k p k +··· 6. 关于复合梯形公式的展开定理 设f(x)在[a,b]区间上无穷次可微，则有如下展开式： T(h)=I+a 1h 2 +a 2h 4 +a 3h 6 +…+a m h 2m +… 式中T(h)是函数f(x)在[a,b]区间上的复化梯形值Tn, 7.3 程序中Mathematica 语句解释 1. 随机函数 Random[] 随机给出闭区间[0,1]内的一个实数 Random[Real, xmax] 随机给出闭区间[0,xmax]内的一个实数 Random[Real, {xmin, xmax}] 随机给出闭区间[xmin,xmax]内的一个实数 Random[Integer] 随机给出整数0或1 Random[Integer, {xmin, xmax}] 随机给出xmin 到xmax 之间的一个整数 Random[Complex] 随机给出单位正方形内的一个复数 2．{a1,a2,…,an} 表示由元素a1,a2,…,an 组成的一个表，元素可以是任何内容。 ) 10(1) ()()(1 1112<<--= q q h F q qh F h F p p
计算机数值模拟是一项综合应用技术，它对教学、科研、设计、产生、管理、决策等部门都有很大的应用价值，为此世界各国均投入了相当多的资金和人力进行研究。其重要性具体体现在以下几个方面： a．从广义上讲，计算机模拟本身就可以看作一种基本试验。计算机计算弹体的侵彻与炸药爆炸过程以及各种非线性波的相互作用等问题，实际上是求解含有很多线性与非线性的偏微分方程、积分方程以及代数方程等的耦合方程组。利用解析方法求解爆炸力学问题是非常困难的，一般只能考虑一些很简单的问题。利用试验方法费用昂贵，还只能表征初始状态和最终状态，中间过程无法得知，因而也无法帮助研究人员了解问题的实质。而数值模拟在某种意义上比理论与试验对问题的认识更为深刻、更为细致，不仅可以了解问题的结果，而且可随时连续动态地、重复地显示事物的发展，了解其整体与局部的细致过程。 b．数值模拟可以直观地显示目前还不易观测到的、说不清楚的一些现象，容易为人理解和分析；还可以显示任何试验都无法看到的发生在结构内部的一些物理现象。如弹体在不均匀介质侵彻过程中的受力和偏转；爆炸波在介质中的传播过程和地下结构的破坏过程。同时，数值模拟可以替代一些危险、昂贵的甚至是难于实施的试验，如反应堆的爆炸事故，核爆炸的过程与效应等。 c．数值模拟促进了试验的发展，对试验方案的科学制定、试验过程中测点的最佳位臵、仪表量程等的确定提供更可靠的理论指导。侵彻、爆炸试验，费用是极其昂贵的，并且存在一定的危险，因此数值模拟不但有很大的经济效益，而且可以加速理论、试验研究的进程。 d．一次投资，长期受益。虽然数值模拟大型软件系统的研制需要花费相当多的经费和人力资源，但和试验相比，数值模拟软件是可以进行拷贝移植、重复利用，并可进行适当修改而满足不同情况的需求。 总之，数值模拟计算已经与理论分析、试验研究成为科学技术探索研究的三个相互依存、不可缺少的手段。正如美国著名数学家拉克斯(P. Lax)所说“科学计算是关系到国家安全、经济发展和科技进步的关键性环节，是事关国家命脉的大事。”
实验一面向微分方程的数值积分法仿真 一、实验目的 1．掌握数值积分法的基本概念、原理及应用； 2．用龙格-库塔法解算微分方程，增加编写仿真程序的能力； 3．分析数值积分算法的计算步长与计算精度、速度、稳定性的关系； 4. 对数值算法中的“病态问题”进行研究。 二、实验内容 1、已知系统微分方程及初值条件 ,(0)1y t y y =+= 取步长0.1h =，试分别用欧拉方程法和RK4法求2t h =时的y 值，并将求得的值与解析解 ()21t y t e t =--比较(将三个解绘于同一坐标中，且用数值进行比较)，说明造成差异的原 因。（①编程完成；②选用MATLAB ode 函数完成。） 程序代码如下： t0=0; tf=2; h=0.1; y1=1; y2=1; y3=1; t1=0; t2=0; t3=0 n=round(tf-t0)/h; for i=1:n y1(i+1)=y1(i)+h*(2*h+y1(i)); t1=[t1,t1(i)+h]; end for i=1:n k1=y2(i)+t2(i); k2=y2(i)+h*k1/2+t2(i)+h/2; k3=y2(i)+h*k2/2+t2(i)+h/2; k4=y2(i)+h*k3+t2(i)+h; y2(i+1)=y2(i)+h*(k1+2*k2+2*k3+k4)/6; t2=[t2,t2(i)+h]; end for i=1:n y3(i+1)=2*exp(t3(i))-t3(i)-1; t3=[t3,t3(i)+h]; end plot(t1,y1,'r',t2,y2,'g',t3,y3,'k') 实验结果如下；
实验10 数值积分 实验目的： 1．了解数值积分的基本原理； 2．熟练掌握数值积分的MATLAB 实现； 3．会用数值积分方法解决一些实际问题。 实验内容： 积分是数学中的一个基本概念，在实际问题中也有很广泛的应用。同微分一样，在《微积分》中，它也是通过极限定义的，由于实际问题中遇到的函数一般都以列表形式给出，所以常常不能用来直接进行积分。此外有些函数虽然有解析式，但其原函数不是初等函数，所以仍然得不到积分的精确值，如不定积分?1 0 d sin x x x 。这时我们一般考虑用数值方法计算其 近似值，称为数值积分。 10.1 数值微分简介 设函数()y f x =在* x 可导，则其导数为 h x f h x f x f h ) ()(lim )(**0* -+='→ （10.1） 如果函数()y f x =以列表形式给出（见表10-1），则其精确值无法求得，但可由下式求得其近似值 h x f h x f x f ) ()()(*** -+≈' （10.2） 一般的，步长h 越小，所得结果越精确。（10.2）式右端项的分子称为函数()y f x =在 *x 的差分，分母称为自变量在*x 的差分，所以右端项又称为差商。数值微分即用差商近似
代替微商。常用的差商公式为： 000()() ()2f x h f x h f x h +--'≈ （10.3） h y y y x f 243)(2 100-+-≈ ' （10.4） h y y y x f n n n n 234)(12+-≈ '-- （10.5） 其误差均为2 ()O h ，称为统称三点公式。 10.2 数值微分的MATLAB 实现 MATLAB 提供了一个指令求解一阶向前差分，其使用格式为： dx=diff(x) 其中x 是n 维数组，dx 为1n -维数组[]21321,, ,n x x x x x x ---,这样基于两点的数值导 数可通过指令diff(x)/h 实现。对于三点公式，读者可参考例1的M 函数文件diff3.m 。 例1 用三点公式计算()y f x =在=x 1.0,1.2,1.4处的导数值，()f x 的值由下表给出。 解：建立三点公式的M 函数文件diff3.m 如下： function f=diff3(x,y) n=length(x);h=x(2)-x(1); f(1)=(-3*y(1)+4*y(2)-y(3))/(2*h); for j=2:n-1 f(j)=(y(j+1)-y(j-1))/(2*h); end f(n)=(y(n-2)-4*y(n-1)+3*y(n))/(2*h);
实验名称：数值积分算法误差分析 1.实验原理 1)欧拉法原理 在数学和计算机科学中，欧拉方法（Euler method）命名自它的发明者莱昂哈德·欧拉，是一种一阶数值方法，用以对给定初值的常微分方程(即初值问题)求解。它是一种解决常微分方程数值积分的最基本的一类显型方法(Explicit method)。 微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值，这个过程称为离散化。实现离散化的基本途径是用向前差商来近似代替导数，这就是欧拉算法实现的依据。欧拉(Euler)算法是数值求解中最基本、最简单的方法，但其求解精度较低，一般不在工程中单独进行运算。所谓数值求解，就是求问题的解y(x)在一系列点上的值y(xi)的近似值yi。对于常微分方程： 0 ) ( ] , [ ), , ( y a y b a x y x f dx dy =∈ = 可以将区间[a,b]分成n段，那么方程在第x i点有y'(x i) = f(x i,y(x i))，再用向前差商近似代替导数则为 ： ，在这里，h是步长，即相邻两个结点间的距离。因此可以根据xi点和yi点的数值计算出y i+1来： 1
，i=0,1,2,L 这就是欧拉格式，若初值y i+ 1是已知的，则可依据上式逐步算出数值解y1,y2,L。 1)龙哥库塔法原理 数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。 龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。该方法主要是在已知方程导数和初值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。令初值问题表述如下。 则，对于该问题的RK4由如下方程给出： 其中 这样，下一个值(y n+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均：
目录 第一章数值积分计算的重述 (1) 1.1引言 (1) 1.2问题重述 (2) 第二章复化梯形公式 (3) 2.1 复化梯形公式的算法描述 (3) 2.2 复化梯形公式在C语言中的实现 (3) 2.3 测试结果 (4) 第三章复化simpson公式 (6) 3.1 复化simpson公式的算法描述 (6) 3.2 复化simpson公式在C语言中的实现 (6) 3.3 测试结果 (7) 第四章复化cotes公式 (8) 4.1 复化cotes公式的算法描述 (8) 4.2 复化cotes公式在C语言中的实现 (9) 4.3 测试结果 (10) 第五章Romberg积分法 (11) 5.1 Romberg积分法的算法描述 (11) 5.2 Romberg积分法在C中的实现 (12) 5.3 测试结果 (13) 第六章结果对比分析和体会 (144) 参考文献 (16) 附录 (16)
数值积分?-10 2 dx e x （一） 第一章 数值积分计算的重述 1.1引言 数值积分是积分计算的重要方法，是数值逼近的重要内容，是函数插值的最直接应用，也是工程技术计算中常常遇到的一个问题。在应用上，人们常要求算出具体数值，因此数值积分就成了数值分析的一个重要内容。在更为复杂的计算问题中，数值积分也常常是一个基本组成部分。 在微积分理论中，我们知道了牛顿-莱布尼茨(Newton-Leibniz)公式 ()() () b a f x d x F b F a =-? 其中()F x 是被积函数()f x 的某个原函数。但是随着学习的深入，我们发现一个问题: 对很多实际问题，上述公式却无能为力。这主要是因为:它们或是被积函数没有解析形式的原函数,或是只知道被积函数在一些点上的值,而不知道函数的形式,对此,牛顿—莱布尼茨(Newton-Leibniz)公式就无能为力了。此外,即使被积函数存在原函数,但因找原函数很复杂,人们也不愿花费太多的时间在求原函数上,这些都促使人们寻找定积分近似计算方法的研究,特别是有了计算机后,人们希望这种定积分近似计算方法能在计算机上实现,并保证计算结果的精度,具有这种特性的定积分近似计算方法称为数值积分。由定积分知识,定积分只与被积函数和积分区间有关,而在对被积函数做插值逼近时,多项式的次数越高,对被积函数的光滑程度要求也越高,且会出现Runge 现象。如7n >时,Newton-Cotes 公式就是不稳定的。因而,人们把目标转向积分区间,类似分段插值,把积分区间分割成若干小区间,在每个小区间上使用次数较低的Newton-Cotes 公式,然后把每个小区间上的结果加起来作为函数在整个区间上积分的近似,这就是复化的基本思想。本文主要
实验三利用数值积分算法的仿真实验 一、实验目的 1) 熟悉MATLAB的工作环境； 2) 掌握MATLAB的 .M文件编写规则，并在命令窗 口调试和运行程序； 3) 掌握利用欧拉法、梯形法、二阶显式Adams法 及四阶龙格库塔法构建系统仿真模型的方法, 并对仿真结果进行分析。 二、实验内容与要求 系统电路如图3所示。电路元件参数：直流电压源，电阻，电感，电容。电路元件初始值：电感电流，电容电压。系统输出量为电容电压。试利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建系统仿真模型，并求出离散系统的输出量响应曲线。连续系统输出响应的解析解为： 其中， 三、实验原理 在连续系统的数字仿真算法中,较常用的有欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法等。欧拉法、梯形法和二阶显式Adams法是利用离散相似原理构造的仿真算法，而显式四阶Runge-Kutta法是利用Taylor级数匹配原理构造的仿真算法。 对于线性系统，其状态方程表达式为: (3-1)式（3-1）中，是系统的n维状态向量,是系统的m维输入向量，是系统的r维输出向量。A为阶参数矩阵，又称动态矩
阵，B为阶输入矩阵，C为阶输出矩阵，D为阶交联矩阵。利用前向欧拉法构建线性系统的仿真模型为： （3-2）式中，为积分步长，为单位矩阵。利用后向欧拉法构建线性系统的仿真模型为： （3-3）利用梯形法构建线性系统的仿真模型为： （3-4）利用二阶显式Adams法构建线性系统的仿真模型为： （3-5） 式中： （3-6）二阶显式Adams法为多步计算方法，利用多步计算方法对系统进行仿真时，需要与之具有相同计算精度的单步计算方法辅助计算。二阶显式Adams法的计算精度为二阶，可以采用梯形法或改进的Euler法等辅助计算。利用改进的Euler法构建线性系统的仿真模型为： （3-7）利用显式四阶Runge-Kutta法构建线性系统的仿真模型为： （3-8） 四、实验步骤与方法 1. 建立系统数学模型 根据图3所示电路，系统状态方程模型: (3-9)式中，状态变量，输出变量，系数矩阵为：
实验一面向方程的数值积分法仿真 1、实验目的 学习用Matlab编写数值积分法仿真程序，深入理解一阶常微分方程组数值积分法的原理和程序结构，加深理解4阶龙格—库塔法的原理以及步长与算法稳定性的关系。 2、实验内容 线性定常系统仿真，二阶系统传递函数为 假设状态初值为零，系统输入为阶跃输入， =1。 （1）根据所学的原理，将传递函数转换为线性状态空间模型（能控标准型），写出状态空间模型表达式。 （2）对状态空间模型采用定步长 RK4 法进行仿真，取步长 h=0.1，仿真时间长度为20秒，做出阻尼比分别为ξ=0，0.1，0.5，1.0，1.5 时的响应曲线，并叠加显示。
答：由于本实验没有控制作用，所以只需要编写三个函数，即“.m”文件。他们分别是 线性系统数值积分法仿真框架函数 function [t,y]=w_LinearSimu(tstart,tstop,h,x0,u0,A,B,C,D) t=[tstart:h:tstop];%t数一个行序列 cntt=size(t,2);%返回列数 cnty=size(C,1);%返回y的维数 y=zeros(cnty,cntt);%构造一个空矩阵，用来存储结果 y0=C*x0+D*u0;% %eval([OutputFile,'(tstart,x0,u0)']);%计算初始输出 y(:,1)=y0 ;%将y0作为输出的第1列 curx=x0; %当前一步的x curu=u0; %当前一步的u cury=y0; %当前一步的y for i=1:1:cntt-1 curx=w_LinearStepIntegral(h,curx,curu,'RK4',A,B);%单步积分运算,针对线性模型
数值积分 第1章 理论依据 逼近论——构造一个简单函数p(x)近似表示f(x),然后对 p(x)求积分得到 f(x)的积分的近似值。基于插值原理，推导出数值积分的基本公式。 §1插值求积公式 为了用数值方法求 b a I(f)=f(x)dx ? ，对被积函数f(x)在给定的n+1个节点 上作Lagrange 插值，用插值函数Pn(x)代替f(x)，就可用I （Pn(x)）构造求积公式，近似地计算定积分I(f(x))。 §2Newton —Cotes 公式 §2.1Newton —Cotes 公式的推导 当§1.1插值求积公式的插值节点为等距节点时，就得到Newton —Cotes 公式。 将区间[a,b]n 等分， b a h n -= ，n+1个节点为 x k =a+kh (k=0,1,…,n) 在节点上对f(x)的Lagrange 插值多项式是： 0()()() n n j n k k j k j j k x x p x f x x x ==≠-=-∑∏ 用P n (x)代替f(x)构造求积公式： 0()()()n n b b j n n k a a k j k j j k x x I p x dx f x dx x x ==≠-==-∑∏?? 记，(k=0,1,…,n) 作代换x=a+th 带入上式，变为： () 00()n n n n k k j j k b a t j A dt b a C n k j =≠? --==--∏?
其中： (k=0,1,…,n) （1-1） 这个积分是有理多项式积分，它与被积函数f(x)和区间[a,b]无关。只要确定n 就能计算出系数。 于是得到称为Newton —Cotes 公式的求积公式： ()0 ()n n n k k k I b a C y ==-∑ （1-2） 其中称为Newton —Cotes 系数。如表1所示。 §2.2Newton —Cotes 公式误差和稳定性 在积分公式中用插值多项式Pn(x)代替f(x)的插值误差是 (1)0 ()()()()()(1)!n n n n k k f R x f x p x x x n ξ+==-=-+∏ 因此，Newton —Cotes 公式的截断误差是 (1)0 ()()()(1)!n n b k a k f R f x x dx n ξ+==-+∏? （1-3） 讨论舍入误差对计算结果产生的影响，设（1-2）式近似计算()b a f x dx ? 其中计算函数值f(xn)有误差值（k=0,1,2, …,n ）。在（1-2）式中令 ? 设计算无误差，舍入误差也忽略，则，由（1-2）式 计算时引式的误差为 () ()()() 0000()[()(())()(...) n n n n n n n k k k k n n n k k e b a C f x C f x b a C C εεε===--+=--++∑∑ 如果皆为正，并设，则 ，故 有
实验一：面向方程的数值积分方法仿真(线性定常系统) 1. 实验目的： 加深理解四阶龙格--库塔法的原理及其稳定性。 2. 实验内容： 对下列系统进行仿真 (1) 线性定常系统 []??????????=??????????+????????????????????---=????? ???????321001600032100300110000110321...x x x y u x x x x x x ， 其初值为：)] 2(1)(1[)(000)0(3)0(2)0(1--?=???? ??????=??????????t t t t u x x x (2) 非线性系统?????+-=-=)()()()() ()()()(t y t bx t sy dt t dy t y t ax t rx dt t dx ① r=0.001,a =2*104，s=0.015, b =10 - 4；x(0)=1200,y(0)=600 ② r=0.001,a =2*10 - 6，s=0.01, b =10 - 6；x(0)=12000,y(0)=600 3. 实验要求： (1) 为保证稳定性，分析系统(1)的最大仿真步长（方法自选，保留两位有效数字）； (2) 设计Matlab 或 C 程序，用四阶龙格库塔法进行仿真计算，改变参数及仿真步长，观察实验结果，寻找最适宜的仿真步长和临界仿真步长。 4. 实验报告： (1) 实验所用程序清单； (2) 实验结果及分析。 注：实验报告可以采用电子文档(标明图号)+书面形式，其中书面报告内容为： (1) 系统(1)最大步长的理论分析； (2) 仿真结果分析； 书面报告中不必列出实验题目与结果图（以下同）。
摘要 本文应用插值积分法和逼近论的思想，简单重述了推导Newton-Cotes公式和Gauss-Legendre求积公式的过程，以及这两个公式的系数、精度等问题。并以这两种数值积分的求解方法为基础，应用quad、guass函数编写具体Matlab程序，通过计算机软件计算出所给题目的近似数值积分。对二者所得的结果进行比较，从而研究了用Newton-Cotes和Gauss-Legendre公式求积分的方法和二者的精确度问题。得知，这两种求积公式所得的结果在精度上的确存在差异，结合理论部分更加充分地说明了，n相同时Gauss-Legendre公式比Newton-Cotes公式具有更高的代数精度，但当代数精度相同时，二者计算的结果仍存在细微的差异。 关键字：插值积分、Newton-Cotes公式、Gauss-Legendre公式
数值积分 1 理论依据 逼近论——构造一个简单函数p(x)近似表示f(x),然后对 p(x)求积分得到 f(x)的积分的近似值。基于插值原理，推导出数值积分的基本公式。 1.1插值求积公式 为了用数值方法求 b a I(f)=f(x)dx ，对被积函数f(x)在给定的n+1个节点上 作Lagrange插值，用插值函数Pn(x)代替f(x)，就可用I（Pn(x)）构造求积公式，近似地计算定积分I(f(x))。
2 Newton —Cotes 公式 2.1 Newton —Cotes 公式的推导 当§1.1插值求积公式的插值节点为等距节点时，就得到Newton —Cotes 公式。 将区间[a,b]n 等分， b a h n -= ，n+1个节点为 x k =a+kh (k=0,1,…,n) 在节点上对f(x)的Lagrange 插值多项式是： 0()()() n n j n k k j k j j k x x p x f x x x ==≠-=-∑∏ 用P n (x)代替f(x)构造求积公式： 0()()()n n b b j n n k a a k j k j j k x x I p x dx f x dx x x ==≠-==-∑∏?? 记错误！未找到引用源。，错误！未找到引用源。(k=0,1,…,n) 作代换x=a+th 带入上式，变为： () 00()n n n n k k j j k b a t j A dt b a C n k j =≠? --= =--∏? 其中：错误！未找到引用源。 (k=0,1,…,n) （1-1） 这个积分是有理多项式积分，它与被积函数f(x)和区间[a,b]无关。只要确定n 就能计算出系数错误！未找到引用源。。 于是得到称为Newton —Cotes 公式的求积公式： ()0()n n n k k k I b a C y ==-∑ （1-2） 其中错误！未找到引用源。称为Newton —Cotes 系数。如表1所示。 表1 Newton —Cotes 系数 n 1 1/ 2 1/2 2 1/6 4/6 1/6
利用Matlab 实现Romberg 数值积分算法 一、内容摘要 针对于某些多项式积分，利用Newton —Leibniz 积分公式求解时有困难，可以采用数值积分的方法，求解指定精度的近似解，本文利用Matlab 中的.m 文件编写了复化梯形公式与Romberg 的数值积分算法的程序，求解多项式的数值积分，比较两者的收敛速度。 二、数值积分公式 1.复化梯形公式求解数值积分的基础是将区间一等分时的Newton —Cotes 求积公式: I =(x)[f(a)f(b)]2 b a b a f dx -≈ +? 其几何意义是，利用区间端点的函数值、与端点构成的梯形面积来近似(x)f 在区间[a,b]上的积分值，截断误差为： 3" (b a)()12 f η-- (a,b)η∈ 具有一次的代数精度，很明显，这样的近似求解精度很难满足计算的要求，因而，可以采用将积分区间不停地对分，当区间足够小的时候，利用梯形公式求解每一个小区间的积分近似值，然后将所有的区间加起来，作为被求函数的积分，可以根据计算精度的要求，划分对分的区间个数，得到复化梯形公式：
I =1 1 (b a)(b a) (x)dx [f(a)f(b)2(a )]2n b a k k f f n n -=--≈+++∑? 其截断误差为： 2" (b a)h ()12 R f η--= (a,b)η∈ 数值积分算法 使用复化的梯形公式计算的数值积分，其收敛速度比减慢，为此，采用Romberg 数值积分。其思想主要是，根据I 的近似值2n T 加上I 与2 n T 的近似误差，作为新的I 的近视，反复迭代，求出满足计算精度的近似解。 用2n T 近似I 所产生的误差可用下式进行估算： 12221 ()3 n n n I T T T -?=-=- 新的I 的近似值： 122 n n j T T -=?+ j =（0 1 2 ….） Romberg 数值积分算法计算顺序 i=0 (1) 002T i=1 (2) 102T (3) 012T
实验三 利用数值积分算法的仿真实验 学号： 姓名： 学院： 一． 实验目的 1) 熟悉MATLAB 的工作环境； 2) 掌握在MATLAB 命令窗口调试运行程序； 3) 掌握M 文件编写规则及在MATLAB 命令窗口运行程序； 4) 掌握利用欧拉法、梯形法、二阶显式Adams 法及四阶龙格库塔法构建 系统仿真模型的方法,并对仿真结果进行分析。 二．实验内容 电路如图１所示电路进行仿真试验。电路元件参数：V E 1=，Ω=10R ， H L 01.0=，F C μ1=。电路元件初始值：A i L 0)0(=，V u c 0)0(=。系统输 出量为电容电压)(t u c 。 DC R C ) (t u c +- ) (t i L 图１ RLC 串联电路 E 三． 实验步骤 1. 求连续系统传递函数 根据所示电路图，我们利用电路原理建立系统的传递函数模型，根据系统的传递函数是在零初始条件下输出量的拉普拉斯变换与输入量的拉普拉斯变换之比，可得该系统的传递函数：
LC Ls R s LC s E s U s G C /1//1)()()(2 ++== 2. 离散系统仿真模型 在连续系统的数字仿真算法中,较常用的有欧拉法、梯形法、二阶显式Adams 法及显式四阶Runge-Kutta 法等。欧拉法、梯形法和二阶显式Adams 法是利用离散相似原理构造的仿真算法，而显式四阶Runge-Kutta 法是利用Taylor 级数匹配原理构造的仿真算法。对于线性系统，其状态方程表达式为: ()()() ()()() t t t t t t ?=+? =+?x Ax Bu y Cx Du 00)(x x =t 式（3-1）中，[]T n t x t x t x )()()(21 =x 是系统的n 维状态向 量,[] T m t u t u t u t )()()()(21 =u 是系统的m 维输入向量， []T r t y t y t y t )()() ()(21 =y 是系统的r 维输出向量。A 为n n ?阶参数矩阵， 又称动态矩阵，B 为m n ?阶输入矩阵，C 为n r ?阶输出矩阵，D 为m r ?阶交联矩阵。根据图所示电路，系统状态方程模型: ) ()()()(t t E t t Cx y B Ax x =+= 式中，状态变量T C L T u i x x ],[],[21==x ，输出变量C u t =)(y ，系数矩阵为： ??????--=0/1/1/C L L R A ，?? ? ???=0/1L B ，[]10=C 。 (1) 欧拉法 利用前向欧拉法构建线性系统的仿真模型为： ()1111 m m m m m m m m h h h ++++?=+=++?? =+??x x x 1A x Bu y Cx Du 式中，h 为积分步长，1为单位矩阵。利用后向欧拉法构建线性系统的仿真模型为：