实验名称 数值积分111
探索实验7 数值积分法
一、 实验目的
了解求积公式及代数精度概念，理解并掌握求定积分的求积公式的算法构造和计算，学习用计算机求定积分的一些科学计算方法和简单的编程技术和能用程序实现这些算法。
二、概念与结论
1. 求积公式：
计算定积分的如下形式的近似公式：
称为求积公式。
2．代数精度：
若求积公式
对一切不高于m 次的 多项都准确成立，而对于m+1次多项式等号不成立，则称此求积公式 的代数精度为m 。
代数精度越高，求积公式越好。
3.求积余项：
4.Newton-Cotes 求积公式的代数精度 n 点Newton-Cotes 求积公式的代数精度至少可以达到n-1,且当n 为奇数时，可以达到n 。
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]区间上无穷次可微，则有如下展开式：
?∑=≈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)()()
10(1)
()()(11112<<--=q q h F q qh F h F p p ?∑=-=b a n
k k k x f A dx x f f R 1)
()()(
?=b a dx
x f I )(T(h)=I+a 1h 2+a 2h 4+a 3h 6+…+a m h 2m +…
式中T(h)是函数f(x)在[a,b]区间上的复化梯形值Tn,
三、程序中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 组成的一个表，元素可以是任何内容。
如：{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 个元素应该也是一个表。
四、方法与程序
在实际问题中，往往会遇到被积函数f(x)的原函数无法用初等函数来表示，或函数只能用表格表示，或有的虽然能用初等函数表示，但过分复杂，所以这些情形都需要去建立定积分的近似计算公式来做积分计算。数值积分是进行定积分计算的一种方法，它可以解决不能用定积分基本公式计算的所有定积分问题。数值积分涉及很多计算公式，这里主要介绍Newton-Cotes 求积公式、复合求积公式、Romberg 求积方法和Monte-Carlo 方法的构造过程和算法程序。
1. n 点 Newton —Cotes 求积公式
n 点 Newton —Cotes 求积公式又称为等距节点求积公式，它是利用被积函数f(x)在积分区间[a,b]的n 个等分节点上的函数值构造的插值函数?(x)代替f(x)做定积分计算所构造求积公式。这个求积公式是通常做定积分近似计算的梯形公式和抛物线公式的推广，主要在理论上用的多些。
1.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 ）公式：
Simpson 求积公式求积余项为
如果想得到其他的 Newton —Cotes 求积公式只要在有关书中查出Cotes 系数表就可以马上得到相应的Newton —Cotes 求积公式。
?∑=-≈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 +-≈?))()2
(4)((6)(b f b a f a f a b dx x f b a +++-≈??∏?∏∑?∏∑??-≠=-≠==≠==--+--=-+---=--=≈10,1)(10,11
,111111)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 ∈''--=ηη],{)()(2880
1)()4(5b a f a b f R ∈--=ηη
1.2 n点 Newton—Cotes求积公式算法：
1．输入被积函数f(x)及积分上下限a,b
2．选择Cotes系数构造求积公式
3．用求积公式求定积分
1.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},
{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系数。
1.4 例题与实验
例1. 用n=3和n=4的Newton-Cotes 求积公式求定积分
的近似值。
解：执行n 点 Newton —Cotes 求积公式程后，在输入的窗口中按提示分别输入
1、3、Exp[-x/2]、3，每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出用n=3的Newton-Cotes 求积公式计算出的定积分结果：
1 2 1
2 (------ + --- + --------- )
6E 3/2 3 E 6 Sqrt[E]
定积分=0.76705953
再执行n 点 Newton —Cotes 求积公式程后，在输入的窗口中按提示分别输入
1、3、Exp[-x/2]、4，每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出用n=4的Newton-Cotes 求积公式计算出的定积分结果：
1 3 3 1
2 (------ + ------ + ------ + ----------)
8E 3/2 8E 7/6 8E 5/6 8 Sqrt[E]
定积分=0.76691628
因此用n=3和n=4的Newton-Cotes 求积公式求本题定积分近似值分别为
0.76705953和0.76691628
注意到本题的精确值为0.766800999….,可见n=4的Newton-Cotes 求积公式计算结果较好。
2. 复化求积公式
复化求积公式是把积分区间分成若干小区间，在每个小区间上采用次数不高的插值公式，如梯形公式或抛物线公式，构造出相应的求积公式，然后再把它们加起来得到整个区间上的求积公式。复化求积公式克服了高次Newton-Cotes 公式计算不稳定的问题，其运算简单且易于在计算机上实现。常用的复化求积公式是复化梯形公式和复化抛物线公式，下面分别讨论。
2.1复化梯形公式的构造过程:
把区间[a,b] n 等分，取节点x k =a+kh,k=0,1,...n, h=(b-a)/n,对每个小区间[x k ,x k+1]用梯形求积公式，再累加起来得：
dx
e x ?-31
2
公式
就是复化梯形公式。它的求积余项为：
由复化梯形公式的余项，可以看到，当n 不断变大时， T n 无限接近定积分值。因此复化梯形公式可以使定积分的计算达到任意精度。为了计算简单，提高效率，常用|T 2n – T n |<ε来得到满足精度要求的定积分值。
2.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 ∑
∑?∑?-=+-=-=++-=+≈=+1111010
))(2)()((2))()((2)()(1n 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)(),(21
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 ∈''--=≤≤''-=-=+-=∑?ηηηη2,))(4)(2)()((6))()(4)((6)()(211110*********
1h
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 +=+++-=++≈=+-=-=+++-=-=∑∑∑
?∑?+],[)(2180)()(2
180)(),()4(441
10)4(45
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 =+++-≈∑∑
?-=-=+111021))(4)(2)()((6)(
由复化抛物线公式的余项，可以看到，当n不断变大时，S n无限接近定积分值。因此复化抛物线公式也可以使定积分的计算达到任意精度，且收敛的速度比复化梯形公式更快。为了计算简单，提高效率，常用| S2n–S n|<ε来得到满足精度要求的定积分值。
2.3 复化梯形求积公式算法：
1.输入被积函数f(x),积分上下限a,b,和求积精度ε
2. n?1, 计算T n
3.计算T2n
4．判断|T2n– T n|<ε是否成立，如果成立，输出定积分近似值，停止
5．否则，T n? T2n , n?2n,转3
2.4 复化抛物线求积公式算法：
1.输入被积函数f(x),积分上下限a,b,和求积精度ε
2. n?1, 计算S n
3.计算S2n
4.判断|S2n– S n|<ε是否成立，如果成立，输出定积分近似值，停止
5.否则，S n? S2n , n?2n,转3
2.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]
说明：本程序从梯形公式T1开始，用复合梯形求积公式求[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 时使用了公式
2.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}]);
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}上的定积分近似值
∑=??
? ??--+-+=n i n n i n a b a f n a b T T 12)12(2221
2.7例题与实验
例2. 用复合梯形公式求定积分
的近似值,要求误差小于10-7。
解：执行复化梯形求积公式程后，在输入的窗口中按提示分别输入：0、3、Exp[-x^2]、10^(-7)，每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出的计算结果如下： n=2 定积分值为0.7313702518 误差=0.0474305
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
n=8
定积分值为-12.06995132
n=16
定积分值为-12.07032146
dx e x ?
-10
2dx
x e x ?π0
cos
定积分值为-12.07034476。
注意到本题的定积分值为-12.07034631639，从中可以看到复合抛物线求积公式收敛还是很快的。你能通过修改复合抛物线求积公式一次得出如上所有计算结果吗？
3.Romberg 求积公式
Romberg 又称逐次分半加速收敛法，它是对复合梯形求积公式采用Richardson 加速技术得到的一种数值求积方法。
3.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) 进行加速，可以有
继续加速下去，可以得到加速序列
式中m 代表对I 进行的加速次数。
此外，注意到T 1(h/2)=T 2n ，且直接计算可以知道T 2(h)就是复合抛物线公式S n ,于是有计算关系式
类似地，还有如下计算格式：
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
442--=Tn
T S n n 1
4)()2(4)(22223--=h T h T h T 1
4)()2(4)(11--=--m m m m m h T h T h T 1
44222--=Sn
S C n 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 值作为积分值。
3.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、
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? T 16n ,
s2? (4t2-t1)/3, t1?t2,
c2? (16s2-s1)/15, s1?s2 ,
R2? (64c2-c1)/63, c1?c2
n ?2n,转4
3.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;
144323--=
Cn
C R n n ε
≤--R R n n 122
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数列的值：T1、T2、T4、T8、S1、S2、S4、C1、C2使，其中t=Table[0,{4}]提供一个4元素的表用于存放T1、T2、T4、T8。
3.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
因此得到满足精度要求的定积分值为3.141592653589794 ，误差=-2.358113704?10-13
4. Monte-Carlo 求积方法
Monte-Carlo 求积方法以概率论大数定理为依据式，借助随机数来求定积分的近似值的一种方法，该方法特别适用于求多重积分和积分区域复杂的重积分，不足之处是收敛较慢。
4.1 Monte-Carlo 求积方法的构造过程:
假设重积分
是有限的，式中Ω?R n ，G(x)是定义在Ω上的多元函数, P(x)是定义在Ω上的分布函数，则有
只要式中x (k)是Ω上关于分布P(x)的n 个任意独立选取的随机点列即可，这个结论可以有概率论中的大数定律：
以概率1成立。
利用以上结果求重积分的方法就称为Monte-Carlo 求积方法。对具体的积分计算有如下计算公式：
1）求定积分的计算公式
式中x k 是[a,b]均匀分布的n 个任意独立选取的随机数列。
2）求二重积分的计算公式
dx x ?
+102
14∑=≈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)()(∑??=≈n k k k D y x f n D d y x f 1
),(||),(σ
式中(x k,y k )是 平面区域D 中均匀分布的n 个任意独立选取的随机点列，|D|表示区域D 的面积。
3）求三重积分的计算公式
式中(x k,y k ,z k )是空间区域Ω中均匀分布的n 个任意独立选取的随机点列，|Ω|表示区域Ω的体积。
Monte-Carlo 求积方法的收敛阶为O(n 1/2).
4.2 Monte-Carlo 求积方法算法：
1.输入被积函数f(x),积分区域边界和随机点的个数n
2.在积分区域中产生n 个随机点列
3.计算被积函数f(x)在随机点列上的函数值并相加求平均
4.用平均值与积分区域的测度之积作为积分近似值，这里测度就是长度、面积、体积之类的度量。
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 后，计算机则给出定积分近似值。
程序中变量说明：
a ：存放积分下限
b: 存放积分上限
f[x]: 存放被积函数f(x)
n: 存放随机点个数n
（*求二重积分*）
Clear[a,b,x,y,n,f];
a=Input["a="]
b=Input["b="]
c=Input["c="]
d=Input["d="]
f[x_,y_]=Input["被积函数f(x,y)="]
n=Input["随机点个数n="];
s=(b-a)*(d-c)*Sum[f[Random[Real,{a,b}], Random[Real,{c,d}]],{n}]/n
Print["n=",n," 积分≈",N[s,10]];
∑???=ΩΩ≈n k k k k z y x f n dv z y x f 1
),,(||),,(
说明：本程序从用Monte-Carlo 求积方法求[a,b]?[c,d]上的定积分近似值。程序执行后，按要求通过键盘输入积分变量x 的下限a 、积分上限b 和积分变量y 的下限c 、积分上限d 、被积函数f(x)和随机点个数n 后，计算机则给出定积分近似值。
程序中变量说明：
a ：存放积分变量x 下限
b: 存放积分变量x 上限
c ：存放积分变量y 下限
d: 存放积分变量y 上限
f[x]: 存放被积函数f(x)
n: 存放随机点个数n
4.4例题与实验
例5.用Monte-Carlo 求积方法求定积分
的近似值,随机点个数取为200。
解：执行Monte-Carlo 求积方法求定积分程序后，在输入的窗口中按提示分别输入：0、1、Sin[x]/x 、200后，每次输入后用鼠标点击窗口的“OK ”按扭，计算机在屏幕上给出的计算结果如下：
n=200 积分=0.9463437539
本题积分准确值为0.94608307036….
例6.用Monte-Carlo 求积方法求定积分
的近似值,用随机点个数分别取为50，60，70，80，…，200的计算值进行计算实验来观察计算结果特点并提出一种从一组Monte-Carlo 求积计算值推断的较好积分近似值的方法。
解：修改Monte-Carlo 求积方法求二重积分程序为：
f[x_,y_]:=Exp[-(x+y)];
Do[s=Sum[f[Random[Real,{1,2}],y=Random[Real,{2,3}]],
{i,1,n},{j,1,n}]/n^2//N;
Print["n=",n," 积分=",N[s,10]];,{n,50,200,10}]
执行后，计算机在屏幕上给出的计算结果如下：
n=50 积分=0.020*******
dx x
Sinx ?
10??--2
132dy
e dx y x
n=60 积分=0.020********
n=70 积分=0.0197868345
n=80 积分=0.01983216355
n=90 积分=0.01993588289
n=100 积分=0.0199346424
n=110 积分=0.01979465142
n=120 积分=0.01987866262
n=130 积分=0.01985729255
n=140 积分=0.01984628761
n=150 积分=0.01991800975
n=160 积分=0.01985555608
n=170 积分=0.01992147037
n=180 积分=0.01982700398
n=190 积分=0.01988543485
n=200 积分=0.01989785087
本题积分准确值为0.0198937375894….,观察计算结果看到积分值并不是严格按照实验点数的增加而单调收敛的，因此有时对某个n值计算得出的结果可能没有比n小的值精确。但它的收敛性是肯定的，不过收敛速度很慢。
观察这16个计算结果值与准确值做比较，看到这样的特点：积分结果从小数点第二位后出现不同，但这16个数在小数点第二位有14个为1、在小数点第三位有14个为9，在小数点第四位有8个为8，它们在相应位数上取值最多，而这些值正好与准确值的相应数相同。因此，我们提出用统计一组计算数相应位数出现的数字规律确定积分结果的较好近似值，即，如果这组树的同一数字位中的某个数字出现的最多，则该数字应为准确值的有效数字，把所有有效数字组合在一起得到的树就是该积分值较好的近似值。
五、思考
1.什么样的求积公式可以求出具有给定精度的定积分值？
2.在给定求积精度ε后，用复合梯形求积公式的余项来事先估计复合梯形求积公式的小区间数n与用判别式|T2n-T n|<ε得出的小区间数n是否一样？用实验说明你的解答。
3.分析n点Newto-Cotes求积程序的编程优缺点，你能给出哪些改进？
4. 分析复合梯形求积程序的编程优缺点，你能给出哪些改进？
5.分析Romber方法程序的编程优缺点，你能给出哪些改进？
6.分析Monte-Carlo方法的编程优缺点，你能给出哪些改进？
六、练习
1分别用n=3，4，8的Newton-Cotes求积公式求定积分
2．求定积分
dx
x
?-
6
2
sin
4
π
dx
e x
?-
1
2
2
π
的近似值,要求误差小于10-7。
1）用复合梯形公式
2）用复合抛物线求积公式
3. .用Romberg 求积方法求定积分
的近似值,要求误差小于10-12。
4． 用Monte-Carlo 求积方法求定积分求定积分
的近似值。用随机点个数分别取为20，30，40，50进行计算。
dx x x dx x x ?
?10312sin )2cos )1dx
x e x ?π0cos
数值分析实验报告1
实验一误差分析 实验1.1（病态问题） 实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出：考虑一个高次的代数多项式 显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动 其中ε（1.1）和（1.221,,,a a 的输出b ”和“poly ε。 （1（2 （3）写成展 关于α solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考Matlab 的帮助。 实验过程： 程序： a=poly(1:20); rr=roots(a); forn=2:21 n form=1:9 ess=10^(-6-m);
ve=zeros(1,21); ve(n)=ess; r=roots(a+ve); -6-m s=max(abs(r-rr)) end end 利用符号函数：（思考题一）a=poly(1:20); y=poly2sym(a); rr=solve(y) n
很容易的得出对一个多次的代数多项式的其中某一项进行很小的扰动，对其多项式的根会有一定的扰动的，所以对于这类病态问题可以借助于MATLAB来进行问题的分析。 学号：06450210 姓名：万轩 实验二插值法
数值积分实验报告
数值分析实验报告 实验四数值积分 一、用复合辛普森和龙贝格算法计算： 复合辛普森主函数xps： function xps(a,b,eps) n=0;Sd=0; S=(f(a)+f(b))*(b-a)/2; while abs(Sd-S)>eps Sd=S; n=n+1; h=(b-a)/n; for i=1:n+1 x(i)=a+(i-1)*h; end S1=f(x(1))+f(x(n+1)); S2=0; S3=0; for i=2:n S2=S2+f(x(i)); end S2=2*S2; for i=1:n S3=S3+f((x(i)+x(i+1))/2); end S3=4*S3; S=(S1+S2+S3)*h/6; end fprintf('%.15f\n',S); 龙贝格主函数romberg2： function romberg2 (a,b,eps) %a,b为区间，eps为精度 Rd=0; R=(b-a)/2*(f(a)+f(b)); N=0; while abs(Rd-R)>eps Rd=R; N=N+1; for k=1:2 if k==1 n=N*2;
else; n=N; end h=(b-a)/n; for i=1:n+1 x(i)=a+(i-1)*h; end C=0; for i=1:n C1=7*f(x(i))+32*f(x(i)+1/4*h)+12*f(x(i)+2/4*h)+32*f(x(i)+3/4*h)+7*f(x(i+1)); C=C+C1*h/90; end if k==1 R=C*64/63; else R=R-C/63; end end end fprintf('结果为：%.15f',R); 1、建立被积函数文件f.m function y=f(x) y=exp(-x^2); 2、调用xps.m、romberg2.m求定积分. >> xps(0,0.5,0.0000001) 0.461281071728228 >>romberg2 (0,0.5,0.0000001) 结果为： 0.461281006413932
数值稳定性验证实验报告
实验课程：数值计算方法专业：数学与应用数学班级：08070141 学号：37 姓名：汪鹏飞 中北大学理学院
实验1 赛德尔迭代法 【实验目的】 熟悉用塞德尔迭代法解线性方程组 【实验内容】 1.了解MATLAB 语言的用法 2.用塞德尔迭代法解下列线性方程组 1234123412341234 54 1012581034 x x x x x x x x x x x x x x x x ---=-??-+--=?? --+-=??---+=? 【实验所使用的仪器设备与软件平台】 计算机，MATLAB7.0 【实验方法与步骤】 1.先找出系数矩阵A ，将前面没有算过的x j 分别和矩阵的(,)A i j 相乘,然后将累加的和赋值给sum ，即(),j s u m s u m A i j x =+?.算 出()/(,) i i x b sum A i i =-,依次循环,算出所有的i x 。 2.若i x 前后两次之差的绝对值小于所给的误差限ε，则输出i x .否则重复以上过程，直到满足误差条件为止. 【实验结果】 (A 是系数矩阵，b 是右边向量，x 是迭代初值，ep 是误差限) function y=seidel(A,b,x,ep) n=length(b); er=1; k=0; while er>=ep
k=k+1; for i=[1:1:n] q=x(i); sum=0; for j=[1:1:n] if j~=i sum=sum+A(i,j)*x(j); end end x(i)=(b(i)-sum)/A(i,i); er=abs(q-x(i)); end end fprintf('迭代次数k=%d\n',k) disp(x') 【结果分析与讨论】 >> A=[5 -1 -1 -1;-1 10 -1 -1;-1 -1 5 -1;-1 -1 -1 10]; b=[-4 12 8 34]; seidel(A,b,[0 0 0 0],1e-3) 迭代次数k=6 0.99897849430002 1.99958456867649 2.99953139743435 3.99980944604109
数值计算实验报告
(此文档为word格式，下载后您可任意编辑修改！) 2012级6班###(学号)计算机数值方法 实验报告成绩册 姓名:宋元台 学号： 成绩：
数值计算方法与算法实验报告 学期： 2014 至 2015 第 1 学期 2014年 12月1日课程名称: 数值计算方法与算法专业:信息与计算科学班级 12级5班 实验编号： 1实验项目Neton插值多项式指导教师：孙峪怀 姓名：宋元台学号：实验成绩： 一、实验目的及要求 实验目的： 掌握Newton插值多项式的算法，理解Newton插值多项式构造过程中基函数的继承特点，掌握差商表的计算特点。 实验要求： 1. 给出Newton插值算法 2. 用C语言实现算法 二、实验内容 三、实验步骤(该部分不够填写.请填写附页)
1.算法分析： 下面用伪码描述Newton插值多项式的算法： Step1 输入插值节点数n，插值点序列{x(i),f(i)},i=1,2,……,n,要计算的插值点x. Step2 形成差商表 for i=0 to n for j=n to i f(j)=((f(j)-f(j-1)(x(j)-x(j-1-i)); Step3 置初始值temp=1,newton=f(0) Step4 for i=1 to n temp=(x-x(i-1))*temp*由temp(k)=(x-x(k-1))*temp(k-1)形成 (x-x(0).....(x-x(i-1)* Newton=newton+temp*f(i); Step5 输出f（x）的近似数值newton(x)=newton. 2.用C语言实现算法的程序代码 #include
计算方法-数值积分实验
实验二数值积分实验 一. 实验目的 （1）熟悉数值积分与数值微分方法的基本思想，加深对数值积分与数值微分方法的理解。 （2）熟悉Matlab编程环境，利用Matlab实现具体的数值积分与数值微分方法。 二. 实验要求 用Matlab软件实现复化梯形方法、复化辛甫生方法、龙贝格方法和高斯公式的相应算法，并用实例在计算机上计算。 三.实验内容 1. 实验题目 已知x e x f x4 sin 1 ) (- + =的数据表 分别编写用复化梯形法、复化辛甫生公法、龙贝格法、三点高斯法求积分?=1 ) (dx x f I 近似值的计算机程序。 A.复化梯形法: a．编写文件Trapezoid.m,代码如下所示：
b．编写文件f2.m: c.运行： B．复化辛甫生公法 a．编写文件FSimpson.m,代码如下所示：
b．编写文件f2.m: function f=f2(x) f=1+exp(-x).*sin(4*x); c.运行： C．龙贝格法
a．编写文件Romberg.m,代码如下所示： b.运行：
D.三点高斯法 a.编写文件TGauss.m文件，如下所示：
b．运行： 2. 设计思想 要求针对上述题目，详细分析每种算法的设计思想。 总体的思想是化复杂为简单的重复 A．复化梯形法使用直接法，通过递归，缩减规模； B．复化辛甫生也是使用直接法，根据公式直接进行编程，通过递归缩减规模； C．龙贝格算法应该在做了的几个中最体现了“化复杂为简单的重复”的思想，多个循环通过变量的适当递增，和一个for循环语句来实现，循环主体只有一句话，但确是整个程序中的亮点和难点； D．三点高斯法直接通过一条简单的公式来编写程序，难度不大； 四.实验体会 对实验过程进行分析总结，对比不同方法的精度，指出每种算法的设计要点及应注意的事项，以及自己通过实验所获得的对数值积分方法的理解。
原料药稳定性试验报告
L- 腈化物稳定性试验报告 一、概述 L-腈化物是L- 肉碱生产过程中的第一步中间体（第二步中间体： L-肉碱粗品；第三步中间体：L-肉碱潮品），由于L- 肉碱生产工艺为 间歇操作，即每生产一步中间体，生产完毕并出具合格检测报告后，存 入中间体仓库，以备下一步生产投料所需。根据本公司L- 肉碱产品的 整个生产周期，L- 腈化物入库后可能存放的最长时间为4 周（约28 天）。以此周期为时间依据制定了L- 腈化物稳定性试验方案，用于验 证L-腈化物在再试验期限内的各项质量指标数据的稳定性，并且能否符 合L- 腈化物的质量标准，此次稳定性试验的整个周期为28 天，具体 的稳定性试验方案以ICH 药物稳定性指导原则为基础制定，以确保L- 腈化化物稳定性试验的可操作性。 二、验证日期 2010 年1 月13 日- 2010 年2 月10 日 三、验证方案 1）样品储存和包装： 考虑到L- 腈化物今后的贮藏、使用过程，本次用于稳定性试验的样品 批次与最终规模生产所用的L- 腈化物的包装和放置条件相同。 2）样品批次选择：此次稳定性试验共抽取三批样品，且抽取样品的批次与 最终规模生产时的合成路线和生产工艺相同
3）抽样频率和日期：从2010.1.13 起，每隔7 天取样一次，共取五次，具体日期为：2010.1.13 、2010.1.20 、2010.1.27 、 2010.2.3 、2010.2.10 ，以确保试验次数足以满足L- 腈化物的稳 定性试验的需要。。 4）检测项目：根据L- 腈化物的质量标准的规定，此次稳定性试验的检测项目共五项，分别为外观、氯含量、熔点、比旋度、干燥失重。这 些指标在L- 腈化物的储存过程中可能会发生变化，且有可能影响 其质量和有效性。 5）试样来源和抽样：L- 腈化物由公司102 车间生产，经检测合格后储存于中间体仓库，本次稳定性试验的L- 腈化物均取自于该中间体仓 库，其抽样方法和抽样量均按照L- 腈化物抽样方案进行抽样。抽 样完毕后直接进行检测分析，并对检测结果进行登记，保存，作为稳 定性数据评估的依据。 四、稳定性试验数据变化趋势分析及评估 通过对三批L- 腈化物的稳定性试验，对其物理、化学方面稳定性资料进行评价，旨在建立未来相似情况下，大规模生产出的L- 腈化物是否适用 现有的再试验期（28天）。批号间的变化程度是否会影响未来生产的
几种定积分的数值计算方法
几种定积分的数值计算方法 摘要:本文归纳了定积分近似计算中的几种常用方法,并着重分析了各种数值方法的计 算思想,结合实例,对其优劣性作了简要说明. 关键词:数值方法;矩形法;梯形法;抛物线法;类矩形;类梯形 Several Numerical Methods for Solving Definite Integrals Abstract:Several common methods for solving definite integrals are summarized in this paper. Meantime, the idea for each method is emphatically analyzed. Afterwards, a numerical example is illustrated to show that the advantages and disadvantages of these methods. Keywords:Numerical methods, Rectangle method, Trapezoidal method, Parabolic method, Class rectangle, Class trapezoid
1. 引言 在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数 )(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式 ?-=b a a F b F x f ) ()()( 求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用. 在科学研究和实际生产中,经常遇到求积分的计算问题,由积分学知识可知,若函数)(x f 在区间],[b a 连续且原函数为)(x F ,则可用牛顿-莱布尼茨公式 ?-=b a a F b F x f ) ()()( 求得积分.这个公式不论在理论上还是在解决实际问题中都起到了很大的作用.另外,对于求导数也有一系列的求导公式和求导法则.但是,在实际问题中遇到求积分的计算,经常会有这样的情况: (1)函数)(x f 的原函数无法用初等函数给出.例如积分 dx e x ?-1 02 , ? 1 sin dx x x 等,从而无法用牛顿-莱布尼茨公式计算出积分。 (2)函数)(x f 使用表格形式或图形给出,因而无法直接用积分公式或导数公式。 (3)函数)(x f 的原函数或导数值虽然能够求出,但形式过于复杂,不便使用. 由此可见,利用原函数求积分或利用求导法则求导数有它的局限性,所以就有了求解数值积分的很多方法,目前有牛顿—柯特斯公式法,矩形法,梯形法,抛物线法,随机投点法,平均值法,高斯型求积法,龙贝格积分法,李查逊外推算法等等,本文对其中部分方法作一个比较． 2.几何意义上的数值算法 s 在几何上表示以],[b a 为底,以曲线)(x f y =为曲边的曲边梯形的面积A ,因此,计 算s 的近似值也就是A 的近似值,如图1所示.沿着积分区间],[b a ,可以把大的曲边梯形分割成许多小的曲边梯形面积之和.常采用均匀分割,假设],[b a 上等分n 的小区间 ,x 1-i h x i +=b x a x n ==,0,其中n a b h -= 表示小区间的长度. 2.1矩形法
《数值计算方法》上机实验报告
《数值计算方法》上机实验报告华北电力大学 实验名称数值il?算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一. 各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程 *对于非线性方程，若已知根的一个近似值，将在处展开成一阶 xxfx ()0, fx ()xkk 泰勒公式 "f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2! 忽略高次项，有 ，fxfxfxxx 0 ()()()，，, kkk 右端是直线方程，用这个直线方程来近似非线性方程。将非线性方程的 **根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkk fx 0 fx 0 0,
解出 fX 0 *k XX,, k' fx 0 k 水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ik fx ()k 八XX, Ikk* fx()k 这就是牛顿迭代公式。 ,2,计算机程序框图:，见, ,3,输入变量、输出变量说明: X输入变量:迭代初值，迭代精度，迭代最大次数，\0 输出变量:当前迭代次数，当前迭代值xkl ,4,具体算例及求解结果: 2/16 华北电力大学实验报吿 开始 读入 l>k /fx()0?,0 fx 0 Oxx，，01* fx ()0 XX,,，?10 kk, ,1，kN, ?xx, 10 输出迭代输出X输出奇异标志1失败标志
,3,输入变量、输出变量说明: 结束 例：导出计算的牛顿迭代公式，并il ?算。（课本P39例2-16） 115cc （0）, 求解结果: 10. 750000 10.723837 10. 723805 10. 723805 2、列主元素消去法求解线性方程组，1,算法原理: 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角 3/16 华北电力大学实验报告方程组求解。 列选主元是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝 aakkkkkk 对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的 两行（包括常ekk 数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结 ,2,计算机程序框图:，见下页, 输入变量:系数矩阵元素，常向量元素baiji 输出变量:解向量元素bbb,,12n
数值分析实验指导 - 7 积分
数值分析实验指导 潘志斌 2014年3月
实验七 数值积分 数值实验综述：通过数值积分实验掌握数值积分的实现，理解各种数值积分公式的特性，并能用数值积分求解积分方程和微分方程。 基础实验 7.1 Newton-cotes 型求积公式 实验目的：学会Newton-cotes 型求积公式,并应用该算法于实际问题. 实验内容：求定积分 ? π cos xdx e x 实验要求：选择等分份数n ,用复化Simpson 求积公式求上述定积分的误差不超过810-的近似值,用MATLAB 中的内部函数int 求此定积分的准确值,与利用复化Simpson 求积公式计算的近似值进行比较。 7.2 Romberg 算法 实验目的：学会数值求积的Romberg 算法,并应用该算法于实际问题. 实验内容：求定积分 ? 1 5 .0dx x 实验要求： （1）要求程序不断加密对积分区间的等分,自动地控制Romberg 算法中的加速收敛过程,直到定积分近似值的误差不超过610-为止,输出求得的定积分近似值。 （2）可用MATLAB 中的内部函数int 求得此定积分的准确值与Romberg 算法计算的近似值进行比较。 7.3 Gauss 型求积公式 实验目的：学会Gauss 型求积公式,并应用该算法于实际问题. 实验内容：求定积分 ? -+4 42 1x dx 实验要求： （1）把Gauss 点的表格存入计算机,以Gauss-Legendre 求积公式作为本实验的例子,要求程序可以根据不同的阶数n ,自动地用n 阶Gauss-Legendre 求积
公式计算上述定积分的近似值.体会Gauss型求积公式是具有尽可能高的代数精度的数值求积公式。 （2）可用MATLAB中的内部函数int求得此定积分的准确值与Gauss型求积公式求得的值进行比较。
计算方法算法的数值稳定性实验报告
专业 序号 姓名 日期 实验1 算法的数值稳定性实验 【实验目的】 1．掌握用MATLAB 语言的编程训练，初步体验算法的软件实现； 2．通过对稳定算法和不稳定算法的结果分析、比较，深入理解算法的数值稳定性及其重要性。 【实验内容】 1．计算积分 ()dx a x x I n ?+=1 0) (n (n=0,1,2......，10) 其中a 为参数，分别对a=0.05及a=15按下列两种方案计算，列出其结果，并对其可靠性，说明原因。 2．方案一 用递推公式 n aI I n 1 1n + -=- (n=1,2,......,10) 递推初值可由积分直接得)1 ( 0a a In I += 3. 方案二 用递推公式 )1 (11-n n I a I n +-= (n=N,N-1,......,1) 根据估计式 ()()()11111+<<++n a I n a n 当1 n a +≥n 或 ()()n 1 111≤<++n I n a 当1 n n a 0+< ≤ 取递推初值为 ()()()() 11212])1(1111[21N +++=++++≈N a a a N a N a I 当1 a +≥ N N 或
()()]1111[21N N a I N +++= 当1 a 0+< ≤N N 计算中取N=13开始 【解】:手工分析怎样求解这题。 【计算机求解】:怎样设计程序？流程图？变量说明？能否将某算法设计成具有形式参数的函数 形式？ 【程序如下】: % myexp1_1.m --- 算法的数值稳定性实验 % 见 P11 实验课题(一) % function try_stable global n a N = 20; % 计算 N 个值 a =0.05;%或者a=15 % %-------------------------------------------- % % [方案I] 用递推公式 %I(k) = - a*I(k-1) + 1/k % I0 =log((a+1)/a); % 初值 I = zeros(N,1); % 创建 N x 1 矩阵(即列向量),元素全为零 I(1) =-a*I0+1; for k = 2:N I(k) =-a*I(k-1)+1/k; end % %--------------------------------------------
数值分析实验报告模板
数值分析实验报告模板 篇一：数值分析实验报告(一)(完整) 数值分析实验报告 1 2 3 4 5 篇二：数值分析实验报告 实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。利用二分法求解给定非线性方程的根，在给定的范围内，假设f(x,y)在[a,b]上连续，f(a)xf(b) 直接影响迭代的次数甚至迭代的收敛与发散。即若x0 偏离所求根较远，Newton法可能发散的结论。并且本实验中还利用利用改进的Newton法求解同样的方程，且将结果与Newton法的结果比较分析。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。掌握二分法的原理，验证二分法，在选对有根区间的前提下，必是收
敛，但精度不够。熟悉Matlab语言编程，学习编程要点。体会Newton使用时的优点，和局部收敛性，而在初值选取不当时，会发散。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b) Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式xk?1?xk?f(xk) f'(xk) 产生逼近解x*的迭代数列{xk}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 xk?1?xk?rf(xk) 'f(xk) 其中r为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x);
数值积分的matlab实现
实验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） 表 10-1 一般的，步长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); 在MATLAB 指令窗中输入指令： x=[1.0,1.1,1.2,1.3,1.4];y=[0.2500,0.2268,0.2066,0.1890,0.1736];diff3(x,y) 运行得各点的导数值为：-0.2470，-0.2170，-0.1890，-0.1650，-0.0014。所以()y f x =在=x 1.0,1.2,1.4处的导数值分别为-0.2470，-0.1890和-0.0014。 对于高阶导数，MATLAB 提供了几个指令借助于样条函数进行求导，详细使用步骤如下： step1：对给定数据点（x,y ），利用指令pp=spline(x,y)，获得三次样条函数数据pp ，供后面ppval 等指令使用。其中，pp 是一个分段多项式所对应的行向量，它包含此多项式的阶数、段数、节点的横坐标值和各段多项式的系数。 step2：对于上面所求的数据向量pp ，利用指令[breaks,coefs,m,n]=unmkpp(pp)进行处理，生成几个有序的分段多项式pp 。 step3：对各个分段多项式pp 的系数，利用函数ppval 生成其相应导数分段多项式的系数，再利用指令mkpp 生成相应的导数分段多项式 step4：将待求点xx 代入此导数多项式，即得样条导数值。 上述过程可建立M 函数文件ppd.m 实现如下： function dy=ppd(pp) [breaks,coefs,m]=unmkpp(pp);
稳定性试验办法
附件3 特殊医学用途配方食品稳定性研究要求（试行） 一、基本原则 特殊医学用途配方食品稳定性研究是质量控制研究的重要组成部分，其目的是通过设计试验获得产品质量特性在各种环境因素影响下随时间 稳定性研究用样品应在满足《特殊医学用途配方食品良好生产规范》要求及商业化生产条件下生产，产品配方、生产工艺、质量要求应与注册申请材料一致，包装材料和产品包装规格应与拟上市产品一致。 影响因素试验、开启后使用的稳定性试验等采用一批样品进行；加速试验和长期试验分别采用三批样品进行。 （二）考察时间点和考察时间
稳定性研究目的是考察产品质量在确定的温度、湿度等条件下随时间变化的规律，因此研究中一般需要设置多个时间点考察产品的质量变化。考察时间点应基于对产品性质的认识、稳定性趋势评价的要求而设置。加速试验考察时间为产品保质期的四分之一，且不得少于3个月。长期试验总体考察时间应涵盖所预期的保质期，中间取样点的设置应当考虑产品的稳定性特点和产品形态特点。对某些环境因素敏感的产品，应适当增加考 3.检验方法：稳定性试验考察项目原则上应当采用《食品安全国家标准特殊医学用途配方食品通则》（GB 29922）、《食品安全国家标准特殊医学用途婴儿配方食品通则》（GB 25596）规定的检验方法。国家标准中规定了检验方法而未采用的，或者国家标准中未规定检验方法而由申请人自行提供检验方法的，应当提供检验方法来源和（或）方法学验证资料。检验方法应当具有专属性并符合准确度和精密度等相关要求。
四、试验方法 （一）加速试验 加速试验是在高于长期贮存温度和湿度条件下，考察产品的稳定性，为配方和工艺设计、偏离实际贮存条件产品是否依旧能保持质量稳定提供依据，并初步预测产品在规定的贮存条件下的长期稳定性。加速试验条件由申请人根据产品特性、包装材料等因素确定。 ％。如在6 温度 ％， 25℃±2℃ 长期试验是在拟定贮存条件下考察产品在运输、保存、使用过程中的稳定性，为确认贮存条件及保质期等提供依据。长期试验条件由申请人根据产品特性、包装材料等因素确定。 长期试验考察时间应与产品保质期一致，取样时间点为第一年每3个月末一次，第二年每6个月末一次，第3年每年一次。 如保质期为24个月的产品，则应对0、3、6、9、12、18、24月样品进行
matlab计算方法实验报告5(数值积分)
计算方法实验报告(5) 学生姓名杨贤邦学号指导教师吴明芬实验时间2014.4.16地点综合实验大楼203 实验题目数值积分方法 实验目的●利用复化梯形、辛普森公式和龙贝格数值积分公式计算定积分的 近似植。 实验内容●梯形、辛普森、柯特斯法及其Matlab实现； ●变步长的梯形、辛普森、柯特斯法及其Matlab实现。 ●题目由同学从学习材料中任意选两题 算法分析梯形：function y=jifeng_tixing(a,b,n,fun) fa=feval(fun,a); fb=feval(fun,b); s=0; h=(b-a)/n; for k=1:n-1 xk=a+k*h; s=feval(fun,xk)+s; end y=(h/2)*(fa+fb+2*s); 辛普生：function y=jifeng_xingpu(a,b,n,fun) fa=feval(fun,a); fb=feval(fun,b); h=(b-a)/n; s=0; s2=feval(fun,a+0.5*h); for k=1:n-1 xk=a+k*h; s=feval(fun,xk)+s; s2=feval(fun,xk+(h/2))+s2; end
与源程序y=(h/6)*(fa+fb+2*s+4*s2); 龙贝格：function r2=jifeng_long(fun,a,b,e) h=b-a; t1=(h/2)*(feval(fun,a)+feval(fun,b)); k=1; r1=10; r2=0; c2=0; while abs(r2-r1)>e; s=0; x=a+h/2; while x=3 r1=r2; c2=s2+(1/15)*(s2-s1); r2=c2+(1/63)*(c2-c1); k=k+1;h=h/2; t1=t2;s1=s2; c1=c2; end end
数值计算实验报告
2012级6班###(学号)计算机数值方法 实验报告成绩册 姓名：安元龙 学号：2012060501 成绩：
数值计算方法与算法实验报告 学期： 2014 至___2015 第 1 学期 2014年 10月26日课程名称:__数值计算方法与算法 __ 专业:信息与计算科学 12级5班实验编号： 1实验项目Neton插值多项式指导教师__孙峪怀姓名：安元龙学号： 2012060501 实验成绩： 一、实验目的及要求 实验目的： 掌握Newton插值多项式的算法，理解Newton插值多项式构造过程中基函数的继承特点，掌握差商表的计算特点。 实验要求： 1. 给出Newton插值算法 2. 用C语言实现算法 二、实验内容 三、实验步骤(该部分不够填写.请填写附页) 1.算法分析： 下面用伪码描述Newton插值多项式的算法： Step1 输入插值节点数n，插值点序列{x(i),f(i)},i=1,2,……,n,要计算的插值点x. Step2 形成差商表 for i=0 to n for j=n to i f(j)=((f(j)-f(j-1)/(x(j)-x(j-1-i)); Step3 置初始值temp=1,newton=f(0) Step4 for i=1 to n temp=(x-x(i-1))*temp/*由temp(k)=(x-x(k-1))*temp(k-1)形成 (x-x(0).....(x-x(i-1)*/ Newton=newton+temp*f(i); Step5 输出f（x）的近似数值newton(x)=newton. #include
稳定性数据评价
稳定性数据评价 1.介绍 1.1 指南的目的 该指南的目的是为了提供如何使用根据ICH指南Q1A（R）里详述的“新原料药和制剂稳定性试验”原则（以后提到即作为总指导原则）而产生的稳定性数据的介绍来建议再试验期或货架期。该指南描述了何时及如何使用有限外推法来建议关于原料药的再试验期或超出来自长期储存条件的数据的观测范围的原料药货架期。 1.2 背景 总指导原则提供的关于稳定性数据的评价和统计分析的指南是性质上简要和范围上有限制。尽管总指导原则指出回归分析是可接收的方法来分析关于再试验期或货架期评价的定量稳定性数据，并建议用0.25显著性水平操作合并批的统计测试，它很少包括细节。另外，总指导原则不包括当复合因素包含在全面或折合-设计调查的情况。当到该方针的第4步，总指导原则的评价部分将会重复，因此删去。 1.3 指南的范围 该指南，总指导原则的附件，目的是当基于定量和定性测试性质的稳定性数据评价而建议再试验期或货架期和贮存条件时提供预期值的清晰解释。该指南概括了基于单个或复合因素和全面或折合-设计调查得出的稳定性数据以确定再试验期或货架期的介绍。ICH Q6A 和Q6B提供了关于调整和证实认可标准的指南。 2. 指南 2.1 一般原则 正规稳定性调查的设计和实行应符合总指导原则列出的原则。稳
定性调查的目的是，在测试最少三批原料药或制剂基础上，确立适用于将来在相似环境下生产和包装批的再试验期或货架期和标签贮存说明。 在稳定性资料的说明和评价里应采用系统性方法，其中应包括，视情况而，从物理、化学、生物和微生物试验，包括从那些与剂型有关的特定性质（例如，固体口服剂型的溶解速率）的结果。如果合适，应注意回顾质量平衡的合适性。应该考虑能引起质量平衡明显不足的因素，例如，降解机理和稳定性-显示能力和分析方法内在可变性。单批的变化程度作用以后生产批次在其再试验期或货架期间仍保留在其认可标准内的信心。 该指南里关于统计法的介绍不意味着当统计计算被证明是多余时，用统计计算仍可取。但在一些情况下统计分析在再试验期或货架期的外推法里是有用的且在其它情况可能提倡将次用于核实再试验期或货架期。 稳定性数据测定的基本原则同于单个-与多个-因素调查和全面-与折合-设计调查。正规稳定性调查里的数据测定，并视情况而定，使用支持数据来确定可能作用原料药或制剂的质量和性能的关键质量性质。应各自评估每个性质和为了建议再试验期或货架期而由调查结果构成的全面评估。所提议的再试验期或货架期不应超过任何单个性质的预测。 附录A里提供的流程图和附录B里提供的关于如何分析和评价从多因素或折合设计得到的关于适当的定量试验性质的长期稳定性数据。用于数据分析的统计方法应该考虑稳定性调查为估计再试验期或货架期而提供有效统计结论。附录B也应该提供关于如何使用再试验
数值分析实验报告总结
数值分析实验报告总结 随着电子计算机的普及与发展，科学计算已成为现代科 学的重要组成部分，因而数值计算方法的内容也愈来愈广泛和丰富。通过本学期的学习，主要掌握了一些数值方法的基本原理、具体算法，并通过编程在计算机上来实现这些算法。 算法算法是指由基本算术运算及运算顺序的规定构成的完 整的解题步骤。算法可以使用框图、算法语言、数学语言、自然语言来进行描述。具有的特征：正确性、有穷性、适用范围广、运算工作量少、使用资源少、逻辑结构简单、便于实现、计算结果可靠。 误差 计算机的计算结果通常是近似的，因此算法必有误差， 并且应能估计误差。误差是指近似值与真正值之差。绝对误差是指近似值与真正值之差或差的绝对值；相对误差：是指近似值与真正值之比或比的绝对值。误差来源见表 第三章泛函分析泛函分析概要 泛函分析是研究“函数的函数”、函数空间和它们之间 变换的一门较新的数学分支，隶属分析数学。它以各种学科
如果 a 是相容范数，且任何满足 为具体背景，在集合的基础上，把客观世界中的研究对象抽 范数 范数，是具有“长度”概念的函数。在线性代数、泛函 分析及相关的数学领域，泛函是一个函数，其为矢量空间内 的所有矢量赋予非零的正长度或大小。这里以 Cn 空间为例， Rn 空间类似。最常用的范数就是 P-范数。那么 当P 取1, 2 ,s 的时候分别是以下几种最简单的情形: 其中2-范数就是通常意义下的距离。 对于这些范数有以下不等式: 1 < n1/2 另外，若p 和q 是赫德尔共轭指标，即 1/p+1/q=1 么有赫德尔不等式: II = ||xH*y| 当p=q=2时就是柯西-许瓦兹不等式 般来讲矩阵范数除了正定性，齐次性和三角不等式之 矩阵范数通常也称为相容范数。 象为元素和空间。女口：距离空间，赋范线性空间， 内积空间。 1-范数: 1= x1 + x2 +?+ xn 2-范数: x 2=1/2 8 -范数: 8 =max oo ，那 外，还规定其必须满足相容性: 所以
数值积分与数值微分实验报告
实验三 数值积分程序设计算法 1）实验目的 通过本次实验熟悉并掌握各种数值积分算法及如何在matlab 中通过设计程序实现这些算法，从而更好地解决实际中的问题。 2）实验题目 给出积分 dx x I ? -= 3 2 2 1 1 1.用Simpson 公式和N=8的复合Simpson 公式求积分的近似值. 2.用复合梯形公式、复合抛物线公式、龙贝格公式求定积分，要求绝对误差为 7 10*2 1-= ε，将计算结果与精确解做比较，并对计算结果进行分析。 3）实验原理与理论基础 Simpson 公式 )]()2 ( 4)([6 b f b a f a f a b S +++-= 复化梯形公式 将定积分? = b a dx x f I )(的积分区间],[b a 分隔为n 等分，各节点为 n j jh a x j ,,1,0, =+= n a b h -= 复合梯形(Trapz)公式为 ])()(2)([21 1 ∑-=++-= n j j n b f x f a f n a b T 如果将],[b a 分隔为2n 等分，而n a b h /)(-=不变， 则 )]()(2)(2)([41 2 111 2b f x f x f a f n a b T n j j n j j n +++-= ∑∑-=+-= 其中 h j a h x x j j )2 1(2 12 1+ +=+ =+ ，)]()(2)(2)([41 2 11 1 2b f x f x f a f n a b T n j j n j j n +++-= ∑∑-=+ -= ∑ -=-++-+ =1 )2) 12((22 1n j n n a b j a f n a b T n=1时，a b h -=，则)]()([2 1b f a f a b T +-= )0(0T = )2 1(2 2 112h a f a b T T + -+ =)1(0T = 若12-=k n ，记)1(0-=k T T n ， ,2,1=k 1 2 --= k a b h jh a x j +=1 2 --+=k a b j a h x x j j 2 12 1+ =+ k a b j a 2 ) 12(-++=，则可得如下递推公式