数值积分上机实验报告_精品文档.docx
《数值积分上机实验报告_精品文档.docx》由会员分享，可在线阅读，更多相关《数值积分上机实验报告_精品文档.docx（17页珍藏版）》请在冰豆网上搜索。
数值积分上机实验报告
141110038桂贤进
题一：
数学上已经证明了
0141+x2dx=π
成立，所以可以通过数值积分来求π的近似值。
1.分别使用复合梯形、复合Simpson求积公式计算π的近似值。
选择不同的h,对每种求积公式，试将误差刻画为h的函数，并比较两方法的精度。
是否存在某个值，当低于这个值之后，再继续减小h的值，计算精度不再有所改进，为什么？
2.实现Romberg求积方法，并重复上面的计算；
3.实现自适应积分方法，并重复上面的计算。
解：
1.1公式分析：
（a）对于复合梯形公式
Tnf=h2[fa+fb+2i=1n-1fa+ih],h=b-an
（1）
离散误差为：
Enf=-nh312f
（2）ξ=-h2b-a12f
（2）ξ,a<ξ（2）
（b）对于复合Simpson公式
Smf=h3[fa+fb+4i=1mfa+2i-1h+2i=1m-1f（a+2ih）]（3）
h=b-a2m=b-an
离散误差为：
Emf=-mh590f（4）ξ=-h4b-a180f（4）ξ,a<ξ1.2实现算法结果：
分别利用梯形公式和Simpson公式计算结果如下：
（下表中E1f=π-Tf,E2f=π-Sf.此处π为MATLAB中的数,可以认为具有足够大的精度）
i
hi
T（f）
E1（f）
S（f）
E2（f）
1
1
3.000000000000000
0.1416
2
1/2
3.100000000000000
0.0416
3.133333333333333
0.0083
4
1/4
3.131176470588235
0.01042
3.141568627450980
2.4026e-05
6
1/6
3.136963066471264
0.00463
3.141591780936043
8.7265e-07
8
1/8
3.138988494491089
0.00260
3.141592502458707
1.5113e-07
10
1/10
3.139925988907159
0.00167
3.141592613939215
3.9651e-08
12
1/12
3.140435246846851
0.00116
3.141592640305380
1.3284e-08
20
1/20
3.141175986954129
4.1667e-04
3.141592652969785
6.2001e-10
30
1/30
3.141407468407330
1.8519e-04
3.141592653535359
5.4434e-11
40
1/40
3.141488486923612
1.0417e-04
3.141592653580105
9.6878e-12
50
1/50
3.141525986923254
6.6667e-05
3.141592653587253
2.5402e-12
100
1/100
3.141575986923129
1.6667e-05
3.141592653589753
3.9968e-14
200
1/200
3.141588486923130
4.1667e-06
3.141592653589793
0
从上表中可以看出：
复合Simpson公式比复合梯形公式精度高，误差收敛的速度快不少。
1.3误差下降速度对比：
从上图可以看出，复合Simpson公式误差的收敛速度比复合梯形公式的误差的收敛速度快不少，下面验证收敛阶。
1.4验证收敛阶：
本实验的实际误差主要由离散误差和计算过程中的舍入误差组成，这里离散误差起主导作用，故理论上实际误差的收敛阶应该与离散误差的收敛阶相同。
下面利用如下公式来计算实际的收敛阶，并与理论分析所得出的离散误差的收敛阶作比较。
err1err2=h1h2p
对上表格中所列的区间长度hi值，逐次利用相邻两个小区间长hi通过上述公式来计算收敛阶,并绘制成图形。
得到图形如下：
（a）对复合梯形公式：
由上面公式
（2）可知，离散误差关于h为二阶收敛，同时由上图可知实验结果的收敛阶将近为2，故与理论分析相符。
（b）对复合Simpson公式：
这里却有些奇怪，由上面公式（4）可知，离散误差理论上为4阶收敛，可实验结果却是将近6阶收敛。
下面将进一步深入探究。
探究如下：
考虑这是由于被积函数f（x）的特殊性导致，而不是由于Simpson公式离散误差真的能达到6阶收敛。
由误差的余项公式
Emf=-mh590f（4）ξ=-h4b-a180f（4）ξ,a<ξ考虑到f（x）=
1.5计算过程中舍入误差的影响
从上表格可以看出：
当区间数n从1到200时，随着h的减小，实际误差在减小。
考虑如下问题：
若不断减小h的值，即不断增加区间数n，是否实际误差会一直减小?
1.5.1理论分析：
我们知道影响实验结果的精确度的因素主要有离散误差和舍入误差，而离散误差的大小可以通过离散误差的余项来体现。
由公式
（2）和（4），可以看出当不断增加区间数，即区间长度h不断减小时，离散误差会越来越小。
但相应地由于计算精度的限制，当h不断减小时，舍入误差却会变大。
故理论上会存在一个阈值H。
当h大于H时，由于离散误差起主导作用，随着区间长度h的减小，实际误差会变小；而当h小于H时，此时舍入误差将在计算结果的精确度起主导作用，再减小h,会导致计算结果的精确度基本保持不变甚至可能会有减小。
1.5.2实验检验：
（a）对于复合梯形公式：
注意到误差收敛的速度较小，故我们首先选取区间数n=100，101,102,103…108进行分析，得到下面图形。
从图中可以看出，复合梯形公式的阈值H在区间数n为10^6到10^8之间取到。
下面对n处于区间10^6到10^8进行分析。
由上图可以看出在n取10^6（对应h=10^-6）左右h达到阈值，故可取阈值H=10^-6.
（b）对于复合Simpson公式：
注意到复合Simpson公式得收敛速度较快，取i从0到100（对应区间数m=2i）,逐次计算出对应于m的实际误差，并作图如下。
从图中可以看出，在i=42（对应m=84,h=1/168）左右，h达到阈值，故可取阈值H=1/168。
通过上述理论分析和实验验证，我们得到如下结果：
使用复合梯形公式和复合Simpson公式计算积分值时，所分成的小区间长h都存在阈值H，当h原因如下：
当h小于H时，此时舍入误差将在计算结果的精确度起主导作用，再减小h，虽然离散误差依旧会减小，但舍入误差会增大，这就导致计算结果的精确度基本保持不变甚至可能会有减小。
2.1公式分析
对Romberg求积方法:
T1,1=h12fa+fb,h1=b-a
Ti,1=12[Ti-1,1+hi-1k=12i-1f（a+k-1hi-1）],i=2,3,…,hi=b-a2i-1=12hi-1
Tm,j=4j-1Tm,j-1-Tm-1,j-14j-1-1,j=2,3,…（m≥j）
分析离散误差为：
对序列{Tm,j}，m=1,2,3…，考虑Em,j=abfxdx-Tm,j,离散误差的数量级为O（h2j），h为序列的小区间长。
（下表中Ef=π-Ti,if,这里π为MATLAB中的数，可以认为是足够精确的）
I
hi
Ti,i
E（Ti,i）
1
1
3
0.1416
2
1/2
3.133333333333333
0.0083
3
1/4
3.142117647058823
5.2499e-04
4
1/8
3.141585783761874
6.8698e-06
5
1/16
3.141592665277717
1.1688e-08
6
1/32
3.141592653638244
4.8451e-11
7
1/64
3.141592653589723
7.0166e-14
8
1/128
3.141592653589793
0
9
1/256
3.141592653589794
8.8818e-16
10
1/512
3.141592653589792
8.8818e-16
11
1/1024
3.141592653589792
1.3323e-15
12
1/2048
3.141592653589794
4.4409e-16
13
1/4096
3.141592653589792
1.3323e-15
14
1/8192
3.141592653589792
1.3323e-15
15
1/16384
3.141592653589794
4.4409e-16
2.2考虑舍入误差的影响：
2.2.1理论分析：
如同1.5小节中所作分析，此处Romberg积分方法同样存在阈值H。
当h小于H时，此时舍入误差将在计算结果的精确度起主导作用，再减小h，虽然离散误差依旧会减小，但舍入误差会增大，这就导致计算结果的精确度基本保持不变甚至可能会有减小。
2.2.2实验验证：
取i=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15，逐次计算E（Ti,i）,并作图如下：
从图中可以看出，当i=9（对应h=1/256）时，再减小h的值（即增大i值），误差基本上不变。
故可取阈值H=1/256.
3.1自适应算法实现：
算法思想：
通过递归调用实现。
（此处感觉课本提供的算法虽然巧妙，但比较冗杂，且会导致占用太多内存。
如用递归调用算法直观上更容易理解。
）
实现算法如下：
（自己感觉还算比较巧妙）
1.脚本文件zishiying.m
clear;
a=0;b=1;
disp（'误差限为：
'）,e=0.0000001
h=（b-a）/2;
f1=f（a）;
f2=f（（a+b）/2）;
f3=f（b）;
tic
S=h/3*（f1+f2+4*f3）;
t=Simpson_auto（a,b,e,S,h/2,f1,f2,f3）;
disp（'近似值为：
'）,t
disp（'误差为：
'）,abs（pi-t）
disp（'耗时为：
'）,toc
2.函数文件Simpson_auto.m
functiony=Simpson_auto（A,B,e,S,h,C1,C2,C3）%将已计算好的函数值传递下去，避免重复计算
f1=f（A+h）;
f2=f（A+3*h）;
%%利用Simpson公式分别计算左半区间和右半区间的近似值
S1=h/3*（C1+C2+4*f1）;
S2=h/3*（C2+C3+4*f2）;
%%判断是否达到所需精度，若未达到将区间分半，进行递归调用
ifabs（S-S1-S2）<10*e
y=S1+S2;
elsey=Simpson_auto（A,（A+B）/2,e/2,S1,h/
升级会员 