数值积分上机实验报告Word格式.docx
- 文档编号:4641936
- 上传时间:2023-05-03
- 格式:DOCX
- 页数:17
- 大小:279.67KB
数值积分上机实验报告Word格式.docx
《数值积分上机实验报告Word格式.docx》由会员分享,可在线阅读,更多相关《数值积分上机实验报告Word格式.docx(17页珍藏版)》请在冰点文库上搜索。
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
从上表中可以看出:
复合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<
b.
考虑到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小于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
Ti,i
E(Ti,i)
3
3.142117647058823
5.2499e-04
3.141585783761874
6.8698e-06
5
1/16
3.141592665277717
1.1688e-08
1/32
3.141592653638244
4.8451e-11
7
1/64
3.141592653589723
7.0166e-14
1/128
9
1/256
3.141592653589794
8.8818e-16
1/512
3.141592653589792
11
1/1024
1.3323e-15
1/2048
4.4409e-16
13
1/4096
14
1/8192
15
1/16384
2.2考虑舍入误差的影响:
2.2.1理论分析:
如同1.5小节中所作分析,此处Romberg积分方法同样存在阈值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);
近似值为:
),t
误差为:
),abs(pi-t)
耗时为:
),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/2,C1,f1,C2)+...
Simpson_auto((A+B)/2,B,e/2,S2,h/2,C2,f2,C3);
end
end
碰到的问题:
写的第一个递归函数,每次执行递归函数都调用了被积函数4次。
这导致了计算精度不够,并且计算量大增。
之后改为每次递归调用,将已计算好的函数值,当作参数直接传递下去,这样每次只需计算新的节点的函数值,避免了重复调用被积函数计算旧点的函数值所带来的问题。
运行结果:
预先给定的误差限e
积分的近似值
实际误差
10-1
2.40261E-05
10-2
1.51131E-07
10-3
10-4
10-5
3.141592651224822
2.36497E-09
10-6
3.141592653552836
3.69571E-11
10-7
3.141592656602465
3.01267E-09
10-8
3.141592653153092
4.36701E-10
10-9
3.141592653452466
1.37327E-10
10-10
3.141592653570822
1.89715E-11
10-11
3.141592653589052
7.41E-13
10-12
3.141592653589790
3.55271E-15
10-13
1.33227E-15
10-14
8.88E-16
10-15
10-16
3.2计算过程中舍入误差的影响:
由于自适应Simpson算法是预先给定误差容限,然后通过算法的执行,保证最终得到的结果满足误差容限。
这里无法将误差刻画为h的函数,也无法分析是否存在阈值H。
从上面的结果表中可以看出实际误差满足预先给定的误差界。
题二:
Plank关于黑体辐射的理论推出积分
0∞x3ex-1dx
用所掌握的所有数值计算的方法计算这个积分,比较不同方法的计算效率和精度。
1.计算该无穷积分的准确值:
推导如下:
0∞x3ex-1dx=0∞x3e-x1-e-xdx=0∞x3e-xn=0∞e-nxdx=0∞n=0∞x3e-n+1xdx=n=1∞0∞x3e-nxdx=n=1∞0∞t3n4e-tdt=n=1∞1n4Γ4=π4903!
=π415≈6.4939394…
2.利用数值积分的方法进行求解:
2.1使用Gauss-Laguerre积分公式:
区间0,+∞中关于权函数Wx=e-x的n+1次直交多项式为:
Ln+1x=exdn+1xn+1e-xdxn+1
对于Laguerre多项式:
L0x=1
L1=1-x
L2=x2-4x+2
L3x=-x3+9x2-18x+6
L4x=x4-16x3+72x2-96x+24
L5x=-x5+25x4-200x3+600x2-600x+120
...
选取Ln+1(x)的n+1个根作为求积基点,便得到计算积分0∞fxe-xdx的n+1点的Gauss--Laguerre求积公式:
0∞e-xf(x)dx≈k=0nAkf(xk)
令fx=x3exex-1=x31-e-x,便可利用上式计算所要求的积分近似值。
已知前五个Gauss—Laguerre求积公式得节点和系数如下:
n
xk
Ak
2-2
(2+√2)/4
2+2
(2-√2)/4
0.415774557
0.711093010
2.294280360
0.278517734
6.289945083
0.103892565E-1
0.322547690
0.603154104
1.745761101
0.357418692
4.536620297
0.388879085E-1
9.395070912
0.539294706E-3
0.263560320
0.521755611
1.413403059
0.398666811
3.596425771
0.759424497E-1
7.085810006
0.361175868E-2
12.640800844
0.233699724E-4
0.222846604
0.458964674
1.188932102
0.417000831
2.992736326
0.113373382
5.775143569
0.103991975E-1
9.837467418
0.261017203E-3
15.982873981
0.898547906E-6
利用公式:
求得结果为:
近似值:
误差:
6.41372746951758
0.0802
6.48113017594569
0.0128
6.49453563566895
5.9623e-04
6.494313366207618
3.7396e-04
6.493941422339378
2.0201e-06
2.2采用截断法
首先,limx→0+x3ex-1=0,故x=0不是瑕点。
补充定义:
fx=0,x=0x3ex-1,x>
其次由上面推导可知:
无穷积分0∞x3ex-1dx收敛。
考虑用区间截去法,将无穷积分转化为在有限区间上的积分。
注意到:
b∞x3ex-1dx≤b∞x3x55!
dx=b∞120x2dx=120b<
ε
对于预先给定的误差容限ε,可取b值为:
b>
120ε.
这里若取误差容限为:
ε=10-7⇒b=1.2*109
但应用mathematica求解:
100∞x3ex-1dx=1.138566453220937×
10-39
50∞x3ex-1dx=1.505582131320634×
10-18
从这里可以看出上面的放缩是很不合理的,选取b值比实际所需值大了太多,造成可后续计算量过大和计算结果失真。
下面取b=50,用求积公式计算050x3ex-1dx的积分值,并与准确值π415进行比较。
(a).复合梯形求积公式
公式如下:
Tnf=h2[fa+fb+2i=1n-1fa+ih],h=b-an
Enf=-nh312f
(2)ξ=-h2b-a12f
(2)ξ,a<
b
这里a=0,b=50,f(x)=x3ex-1,且f(0)=0.
求解结果为:
(b).复合Simpson求积公式
Smf=h3[fa+fb+4i=1mfa+2i-1h+2i=1m-1f(a+2ih)]
b.
这里a=0,b=50,f(x)=x3ex-1.
(c).使用Gauss-Legendre求积公式
If=050x3ex-1dx
对上面积分,先将区间由[0,50]变换为[-1,1],令x=25*(t+1),则
If=-11254*t+13e25t+1-1dt
已知而前6个高斯—勒让德求积公式结点和系数为:
±
0.577350269
0.774596669
0.555555556
0.861136312
0.347854854
0.339981044
0.652145155
0.906179846
0.236926885
0.538469310
0.478628670
0.932469514
0.171324492
0.661209386
0.360761573
0.238619186
0.467913935
0.949107912
0.129484966
0.741531186
0.279705391
0.405845151
0.381830051
0.417959184
(d).Romberg求积公式
使用Romberg求积序列进行计算积分
可得结果如下:
6.02734E-16
6.4939394
7.2333E-06
6.49393217
0.129399609
6.36453979
4.287105281
2.20683412
7.40969437
0.91575497
6.538473388
0.04453399
6.492829841
0.00110956
6.493942337
2.9344E-06
6.493939407
4.6484E-09
6.493939402
8.5727E-12
1.3323E-14
4.4409E-15
3.1086E-14
2.8422E-14
4.1744E-14
其中Ti,i与准确值之间的误差见下图。
(e).自适应Simpson公式
在区间[0,50]上进行积分的结果:
预先设置误差限
近似值
0.1
0.121312585
6.3726268
0.01
6.49498615
0.0010467
0.001
6.494024148
8.475E-05
0.0001
6.493917421
2.198E-05
0.00001
6.493940614
1.212E-06
1E-06
6.493939661
2.584E-07
1E-07
6.493939624
2.219E-07
1E-08
6.493939403
1.145E-09
1E-09
1.337E-10
1E-10
1.077E-11
1E-11
1.585E-12
1E-12
6.59E-13
1E-13
1.776E-14
1E-14
2.665E-15
1E-15
1.776E-15
结论:
从上面分析可以看出,如果不用
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 积分 上机 实验 报告