数值积分上机实验报告.docx
- 文档编号:2503124
- 上传时间:2023-05-03
- 格式:DOCX
- 页数:17
- 大小:279.67KB
数值积分上机实验报告.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/2,C1,f1,C2)+... Simpson_auto((A+B)/2,B,e/2,S2,h/2,C2,f2,C3); end end 碰到的问题: 写的第一个递归函数,每次执行递归函数都调用了被积函数4次。 这导致了计算精度不够,并且计算量大增。 之后改为每次递归调用,将已计算好的函数值,当作参数直接传递下去,这样每次只需计算新的节点的函数值,避免了重复调用被积函数计算旧点的函数值所带来的问题。 运行结果: 预先给定的误差限e 积分的近似值 实际误差 10-1 3.141568627450980 2.40261E-05 10-2 3.141592502458707 1.51131E-07 10-3 3.141592502458707 1.51131E-07 10-4 3.141592502458707 1.51131E-07 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 3.141592653589792 1.33227E-15 10-14 3.141592653589792 8.88E-16 10-15 3.141592653589793 0 10-16 3.141592653589793 0 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 1 2-2 (2+√2)/4 2+2 (2-√2)/4 2 0.415774557 0.711093010 2.294280360 0.278517734 6.289945083 0.103892565E-1 3 0.322547690 0.603154104 1.745761101 0.357418692 4.536620297 0.388879085E-1 9.395070912 0.539294706E-3 4 0.263560320 0.521755611 1.413403059 0.398666811 3.596425771 0.759424497E-1 7.085810006 0.361175868E-2 12.640800844 0.233699724E-4 5 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 利用公式: 0∞e-xf(x)dx≈k=0nAkf(xk) 求得结果为: n 近似值: 误差: 1 6.41372746951758 0.0802 2 6.48113017594569 0.0128 3 6.49453563566895 5.9623e-04 4 6.494313366207618 3.7396e-04 5 6.493941422339378 2.0201e-06 2.2采用截断法 首先,limx→0+x3ex-1=0,故x=0不是瑕点。 补充定义: fx=0,x=0x3ex-1,x>0 其次由上面推导可知: 无穷积分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<ξ 这里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)] h=b-a2m=b-an 离散误差为: Emf=-mh590f(4)ξ=-h4b-a180f(4)ξ,a<ξ 这里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个高斯—勒让德求积公式结点和系数为: n xk Ak 1 ±0.577350269 1 2 ±0.774596669 1 0 0.555555556 3 ±0.861136312 0.347854854 ±0.339981044 0.652145155 4 ±0.906179846 0.236926885 ±0.538469310 0.478628670 5 ±0.932469514 0.171324492 ±0.661209386 0.360761573 ±0.238619186 0.467913935 6 ±0.949107912 0.129484966 ±0.741531186 0.279705391 ±0.405845151 0.381830051 0 0.417959184 (d).Romberg求积公式 使用Romberg求积序列进行计算积分 If=050x3ex-1dx 可得结果如下: i hi Ti,i E(Ti,i) 1 1 6.02734E-16 6.4939394 2 1/2 7.2333E-06 6.49393217 3 1/4 0.129399609 6.36453979 4 1/8 4.287105281 2.20683412 5 1/16 7.40969437 0.91575497 6 1/32 6.538473388 0.04453399 7 1/64 6.492829841 0.00110956 8 1/128 6.493942337 2.9344E-06 9 1/256 6.493939407 4.6484E-09 10 1/512 6.493939402 8.5727E-12 11 1/1024 6.493939402 1.3323E-14 12 1/2048 6.493939402 4.4409E-15 13 1/4096 6.493939402 3.1086E-14 14 1/8192 6.493939402 2.8422E-14 15 1/16384 6.493939402 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 6.493939402 1.337E-10 1E-10 6.493939402 1.077E-11 1E-11 6.493939402 1.585E-12 1E-12 6.493939402 6.59E-13 1E-13 6.493939402 1.776E-14 1E-14 6.493939402 2.665E-15 1E-15 6.493939402 1.776E-15 结论: 从上面分析可以看出,如果不用
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 积分 上机 实验 报告