数学实验——插值与数值积分.docx

数学实验报告实验3 插值与数值积分实验3 插值与数值积分分1 黄浩 2011011743一、 实验目的1) 掌握用matlab计算Lagrange、分段线性、分段三次和三次样条插值这四种插值方法，并通过改变节点的数目，分析插值结果2) 掌握用matlab及梯形公式、辛普森公式计算数值积分3) 通过实例学习用插值和数值积分解决实际问题二、 实验内容1.（1）题：问题叙述：考虑函数：fx=11+x2 x-5,5,取不同的节点数目,分别用Lagrange 插值、分段线性和三次插值，以及三次样条插值近似f(x)，分析插值效果方案实施：由于在matlab中，没有预设的lagrange插值法函数，因此我首先编制了Lagrange（程序见四.1）和分段线性的插值函数（程序见四.2），并进行了y=x2函数的插值检验（程序见四.3），以检验编写的函数是否准确，如下图所示：上图中，红、绿两线分别是自设Lagrange和分段线性的插值图，蓝色十字点是matlab自带的分段线性函数的插值图，蓝绿完全吻合、红绿基本吻合，证明自设函数是正确的。然后进行本题目要求的插值操作：以（x0,y0）为节点，节点在区间内均匀分布且将节点数分别调为6、11、21；x为插值点，绘图时使用101个插值点（始终大于节点数），制表时使用的插值点数目依节点数调节；y为真值，y1、y2、y3、y4分别为使用Lagrange、分段线性、分段三次和三次样条插值的方法所得的插值结果。1) 节点数为6时：（即节点为-5,-3,3,5）xyy1y2y3y401.00000.56730.50000.50000.56850.50000.80000.55010.50000.50000.55131.00000.50000.50000.50000.50000.50001.50000.30770.42130.40000.44250.41672.00000.20000.32120.30000.31330.31312.50000.13790.20970.20000.17750.20293.00000.10000.10000.10000.10000.10003.50000.07550.00780.08460.07540.01834.00000.0588-0.04810.06920.0559-0.02854.50000.0471-0.04600.05380.0431-0.02635.00000.03850.03850.03850.03850.0385图像为： 由上图可见，当节点数为6时，插值效果皆不理想，在各节点之间的插值都有很大的误差，在x=0附近误差最大2) 节点数为11时：（即节点为-5,-4,4,5）xyy1y2y3y401.00001.00001.00001.00001.00000.50000.80000.84340.75000.79690.82051.00000.50000.50000.50000.50000.50001.50000.30770.23530.35000.32190.29732.00000.20000.20000.20000.20000.20002.50000.13790.25380.15000.13850.14013.00000.10000.10000.10000.10000.10003.50000.0755-0.22620.07940.07550.07454.00000.05880.05880.05880.05880.05884.50000.04711.57870.04860.04650.04845.00000.03850.03850.03850.03850.0385图像为： 由上图可见，当节点数为11时，分段三次与三次样条插值的效果较好，其次是分段线性插值，而Lagrange插值在-2,2区间内精度较高，但在|x|3时，就会发生Runge振荡3) 节点数为21时（程序见四.4）：（即节点为-5,-4.5,4.5,5）xyy1y2y3y401.00001.00001.00001.00001.00000.10000.99010.99040.96000.98690.98910.20000.96150.96260.92000.95260.95940.30000.91740.91890.88000.90500.91520.40000.86210.86320.84000.85150.86060.50000.80000.80000.80000.80000.80000.60000.73530.73370.74000.74560.73720.70000.67110.66820.68000.68230.67420.80000.60980.60650.62000.61630.61270.90000.55250.55030.56000.55350.55411.00000.50000.50000.50000.50000.50001.10000.45250.45500.46150.45440.45171.20000.40980.41430.42310.41180.40891.30000.37170.37650.38460.37280.37111.40000.33780.34110.34620.33790.33751.50000.30770.30770.30770.30770.30771.60000.28090.27700.28620.28130.28101.70000.25710.25010.26460.25750.25721.80000.23580.22810.24310.23600.23591.90000.21690.21150.22150.21690.21702.00000.20000.20000.20000.20000.20002.10000.18480.19180.18760.18500.18482.20000.17120.18430.17520.17130.17122.30000.15900.17420.16280.15900.15892.40000.14790.15900.15030.14790.14792.50000.13790.13790.13790.13790.13792.60000.12890.11300.13030.12890.12892.70000.12060.08960.12280.12070.12062.80000.11310.07500.11520.11310.11312.90000.10630.07700.10760.10630.10633.00000.10000.10000.10000.10000.10003.10000.09430.14080.09510.09430.09423.20000.08900.18580.09020.08900.08903.30000.08410.21010.08530.08410.08413.40000.07960.18230.08040.07960.07963.50000.07550.07550.07550.07550.07553.60000.0716-0.11430.07210.07160.07163.70000.0681-0.34590.06880.06810.06813.80000.0648-0.51360.06550.06480.06483.90000.0617-0.44590.06220.06170.06174.00000.05880.05880.05880.05880.05884.10000.05611.13530.05650.05620.05614.20000.05362.67440.05410.05370.05364.30000.05134.06910.05180.05130.05134.40000.04913.94510.04940.04910.04914.50000.04710.04710.04710.04710.04714.60000.0451-10.33450.04530.04510.04514.70000.0433-28.66260.04360.04330.04334.80000.0416-50.86440.04190.04150.04164.90000.0400-58.23810.04020.03990.04005.00000.03850.03850.03850.03850.0385图像为：由上图可知，当节点数为21时，分段三次和三次样条插值的精度都很高，分段线性稍次，然而Lagrange插值在x=0附近的区间内精度很好，在|x|3.5时，却发生了更为剧烈的Runge振荡（见上表蓝色粗字，在x=4.9时，插值竟达到了-58.2381）为定量比较各种算法的误差大小，在节点数为21的情况下，在不发生Runge振荡的1.1,2.0的区间内，对上表进行数据处理，结果如下：（其中y1，y2，y3，y4分别为使用Lagrange、分段线性、分段三次和三次样条插值的方法所得结果的相对误差的绝对值）xy1y2y3y41.10000.00550.01430.01540.00591.20000.01100.02120.02670.00701.30000.01290.02150.03070.00461.40000.00980.01500.02400.00121.50000.00000.00000.00000.00001.60000.01390.03320.01710.00111.70000.02720.05800.02680.00121.80000.03270.06580.02920.00041.90000.02490.04730.02080.00052.00000.00000.00000.00000.0000SUM0.13790.27620.19070.0218由上表看出，在无Runge振荡的情况下，使用不同插值方法，误差从大到小排列为：分段线性 分段三次 Lagrange 三次样条插值。得出结论：由上述的比较，我们可以总结出以下几点结论：a) 一般而言，节点数增加时，插值精度普遍提高b) 在不发生Runge振荡时，分段三次、三次样条插值和Lagrange插值的精度近似，都有很好的插值效果，而分段线性插值次于上述三者，具体的误差大小为：分段线性 分段三次 Lagrange 三次样条插值c) 当节点数增加时，Lagrange插值法的Runge振荡会更加剧烈，在区间端点附近会发生更大的误差，但在区间内部插值效果较好综合上述三点，还可以对该四种插值方法做一个整体评价：a) Lagrange插值法：在节点数较少时，插值精度差，当节点数增加时，又有可能发生Runge振荡，造成更大的误差，这是它作为高次插值多项式的内在缺陷。因此，在我们无法评估节点数与所选区间是否合适，也无法判断插值的突跃是否为Runge振荡时，它的精度无法得到保证，可信性较低。b) 分段线性插值法：其精度随着节点数的增加而单调地增加，但在同种情况下，其精度低于分段三次与三次样条插值。因此在计算机计算时，使用后二者更有利于保证精度，但是在广泛的工程领域中，分段线性因其计算的简便性，使得工程师可以通过查表进行插值，因而比后二者有着更广泛的应用。c) 分段三次和三次样条插值：精度随着节点数的增加而单调地增加，且明显高于分段线性插值，从精度上来说是更优化的选择，但是计算量也较大。1.（2）题问题叙述：对于上一问的函数fx=11+x2 x-5,5，选取非均匀插值节点：xk=5cos2k-12n,k=1,2,n,尝试进行插值，并分析插值效果方案实施：选节点数为20，即n=20，非均匀选取；插值点21个，公差为0.5，节点与插值点互不相同；采用与上一问相似的编程语言（程序见四.5），插值结果如下：xyy1y2y3y401.00000.96240.86660.86660.94680.50000.80000.81260.80500.84820.81701.00000.50000.51190.51900.51490.51561.50000.30770.29630.33020.31150.29812.00000.20000.20280.20380.20030.20132.50000.13790.14050.14180.13790.13863.00000.10000.09640.10270.10000.09963.50000.07550.07830.07660.07550.07564.00000.05880.05670.05930.05880.05884.50000.04710.04870.04720.04710.04715.00000.03850.0370NaN0.03850.0385图像为：（其中“+”为节点）由上图可见，在使用这种非均匀插值节点的情况下，没有Runge振荡发生。和上一问类似，也对上表进行处理，以定量比较误差大小：xy1y2y3y40.00000.03760.09950.00000.09250.50000.01570.00940.05370.03681.00000.02380.01390.00790.00141.50000.03700.11440.05660.04302.00000.01400.00490.01720.00502.50000.01890.00930.02750.00513.00000.03600.06540.02630.00403.50000.03710.02170.01440.00134.00000.03570.04590.00840.00004.50000.03400.03080.00210.0000SUM0.28980.41510.21410.1891从上表可以看出，在此情况下，误差从大到小排列为：分段线性 Lagrange 分段三次 三次样条插值得出结论：由该图我们可以得到如下几点信息：a) 在这种情况下，插值精度从大到小为：分段线性 Lagrange 分段三次 三次样条插值。该结果与上一问的不同之处，在于Lagrange和分段三次的位次发生了对调，这不难理解，因为上一问是在无Runge振荡的情况下窄区间内排序的，而本次是在全插值区间进行了比较，两者并不矛盾。b) 当使用这种非均匀插值方法时，在整个插值区间内，都没有发生Runge振荡，也没有任何振荡的痕迹，是对Lagrange插值的优化。2.题问题叙述：用数值积分方法计算函数x=-x12e-t22dt在一些节点上的近似值，自行编制一张x的取值表。根据此表编写一个函数，对任意xR，给出(x)的近似值，同时输出近似值的误差上界。将你的(x)近似值与Matlab中normpdf函数的计算结果进行对比。方案实施：为了训练matlab的编程能力，我首先编制了均匀分段的梯形求积函数（程序见四.6）和均匀分段的辛普森求积函数（程序见四.7），为了探究二者是否能准确计算积分，我借助于matlab的预设函数trapz()与quad()，通过计算x=0x12e-t22dt，进行了简易的检验（程序见四.8），结果如下：xy1y2y3y40.01000.00400.00400.00400.00400.02000.00800.00800.00800.00800.03000.01200.01200.01200.01200.04000.01600.01600.01600.01600.05000.01990.01990.01990.01990.06000.02390.02390.02390.02390.07000.02790.02790.02790.02790.08000.03190.03190.03190.03190.09000.03590.03590.03590.03590.10000.03980.03980.03980.0398上表中的y2,y4,y1,y3分别为自编的梯形、辛普森求积和matlab预设函数trapz()、quad()的积分值。显而易见，自编的两个求积公式都是正确的。 然后进行题目要求的积分操作：1) 编制标准正态分布函数表（程序见四.9）：u0.000.010.020.030.040.050.060.070.080.090.00.50000.50400.50800.51200.51600.51990.52390.52790.53190.53590.10.53980.54380.54780.55170.55570.55960.56360.56750.57140.57530.20.57930.58320.58710.59100.59480.59870.60260.60640.61030.61410.30.61790.62170.62550.62930.63310.63680.64060.64430.64800.65170.40.65540.65910.66280.66640.67000.67360.67720.68080.68440.68790.50.69150.69500.69850.70190.70540.70880.71230.71570.71900.72240.60.72570.72910.73240.73570.73890.74220.74540.74860.75170.75490.70.75800.76110.76420.76730.77040.77340.77640.77940.78230.78520.80.78810.79100.79390.79670.79950.80230.80510.80780.81060.81330.90.81590.81860.82120.82380.82640.82890.83150.83400.83650.83891.00.84130.84380.84610.84850.85080.85310.85540.85770.85990.86211.10.86430.86650.86860.87080.87290.87490.87700.87900.88100.88301.20.88490.88690.88880.89070.89250.89440.89620.89800.89970.90151.30.90320.90490.90660.90820.90990.91150.91310.91470.91620.91771.40.91920.92070.92220.92360.92510.92650.92790.92920.93060.93191.50.93320.93450.93570.93700.93820.93940.94060.94180.94290.94411.60.94520.94630.94740.94840.94950.95050.95150.95250.95350.95451.70.95540.95640.95730.95820.95910.95990.96080.96160.96250.96331.80.96410.96490.96560.96640.96710.96780.96860.96930.96990.97061.90.97130.97190.97260.97320.97380.97440.97500.97560.97610.97672.00.97720.97780.97830.97880.97930.97980.98030.98080.98120.98172.10.98210.98260.98300.98340.98380.98420.98460.98500.98540.98572.20.98610.98640.98680.98710.98750.98780.98810.98840.98870.98902.30.98930.98960.98980.99010.99040.99060.99090.99110.99130.99162.40.99180.99200.99220.99250.99270.99290.99310.99320.99340.99362.50.99380.99400.99410.99430.99450.99460.99480.99490.99510.99522.60.99530.99550.99560.99570.99590.99600.99610.99620.99630.99642.70.99650.99660.99670.99680.99690.99700.99710.99720.99730.99742.80.99740.99750.99760.99770.99770.99780.99790.99790.99800.99812.90.99810.99820.99820.99830.99840.99840.99850.99850.99860.99863.00.99860.99870.99870.99880.99880.99890.99890.99890.99900.99903.10.99900.99910.99910.99910.99920.99920.99920.99920.99930.99933.20.99930.99930.99940.99940.99940.99940.99940.99950.99950.99953.30.99950.99950.99950.99960.99960.99960.99960.99960.99960.99973.40.99970.99970.99970.99970.99970.99970.99970.99970.99970.99983.50.99980.99980.99980.99980.99980.99980.99980.99980.99980.99983.60.99980.99980.99990.99990.99990.99990.99990.99990.99990.99993.70.99990.99990.99990.99990.99990.99990.99990.99990.99990.99993.80.99990.99990.99990.99990.99990.99990.99990.99990.99990.99993.91.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000经过与概率论与数理统计(第一版，茆诗松主编)一书中表格的比对，发现上表的近似值是完全正确的。2) 编写所求的matlab函数由于自适应求积函数quad()会自动调节节点以减小误差，因而无法通过公式计算输出误差，因此作者采用了梯形求积公式trap(x,y)来进行计算。又考虑到正态分布关于y轴对称，因此可通过计算x点与零点之间的定积分来转化为所求的广义积分，于是无穷限一端可能引入的端误差被消除。在课本中给出了梯形求积公式的误差上界公式：I-Tnh212M2b-a，其中M2=maxfx,x(a,b)，但对于数值积分，仍然存在一个问题，即步长与端点之间“缝隙”的误差：例如，如果使用了0.001的步长，但若求x=1.0009时的广义积分，则trapz()求出的数值积分便未包括区间1.000,1.0009（如下图所示）,这会引入另一种误差。因此，作者认为，该情况下会有三种误差，一是机器精度误差，二是梯形面积带来的误差，三是区间“缝隙”引入的误差（如上图）。从直观上考虑，在步长为0.001时，若x小数位数大于三位，则第三种误差会大一些，我通过实验事实证明了上述猜测（函数见四.10，检验程序见四.11）：xy1（积分值）y2（二、三类误差）y3（第三类误差）10.8413447259043150.00000003324519001.10.8643339190838340.00000003656970901.220.8887675432819440.00000004055913201.3330.9087340798726740.00000004431583801.44440.9256304954052940.0000562903226320.0000562423032801.555550.9400269499055090.0000655182825530.000065466567998在上表中，x是积分右限，y1是数值积分值，y2是第二、三类误差的和，y3是第三类误差。结果证明上述直观的想法是正确的，当小数位数小于等于三位时，由于无“区间缝隙”，仅有第二类误差，数量级在10-8左右；而当小数位数大于三位时，在尾部端点产生缝隙，该区间未算入数值积分中，因而产生了第三类误差，数量级较大，在10-5左右，此时总误差仅取决于第三类误差。而因为matlab都是双精度计算，第一类误差在10-16左右，远小于后两种误差，不予考虑。于是，在使用trapz函数的基础上，为了消除第三类误差，我又修改了函数，使输出的积分结果为上表所示的（y1+y3），通过单向校正消除误差；同时考虑到当x8时的积分值已经十分趋近0和1，为避免由于x很大造成的大数计算，节省计算时间，我又增添了两个elseif语句，将这种情况下的积分值直接输出为0和1，并将x=8时的误差当做此时的误差上界。经过以上两步修改，即得到最终题目所要求的函数（程序见四.12），并与matlab中的预设函数normcdf做了对比（程序见四.13）,试验结果如下表所示：xy1（自设函数）y2（误差上限）y3（normcdf）10.8413447259043150.0000000332451900.8413447460685431.10.8643339190838340.0000000365697090.8643339390536171.220.8887675432819440.0000000405591320.8887675625521651.3330.9087340798726740.0000000443158380.9087340980995581.44440.9256867377085740.0000000480193520.9256867546330211.555550.9400924164735070.0000000517145550.9400924319044251.6666660.9522095675829620.0000000554086280.95220958140905600.50000000000000000.500000000000000-1900.0000003324519000.00000000000000030.9986501008604080.0000000997355700.99865010196837080.9999999999999990.0000002659615200.999999999999999将上表再做整理，以验证误差上限的确小于y3与y1的差值：y3-y1（实际误差）y2（输出的误差上限）0.00000002020.00000003320.00000002000.00000003660.00000001930.00000004060.00000001820.00000004430.00000001690.00000004800.00000001540.00000005170.00000001380.00000005540.00000000000.00000000000.00000000000.00000033250.00000000110.00000009970.00000000000.0000002660 由上表可以看出，自设函数输出的积分值与预设函数normcdf的值十分接近，而且所给的误差上限也大于实际的误差，因此，编制的函数是符合题目要求的。三、 实验总结本次“插值与数值积分”实验是第一次数学实验，在实验中遇到了两个比较大的问题：第一个是matlab的使用，一开始对于矩阵、基本运算和常用函数都不太熟悉，所以本次实验的前期工作耗费了较长的时间；第二个是在第二题中的误差分析，这一块需要仔细地分析，否则很容易漏掉第三类误差的校正，而且误差公式也很复杂，涉及到求导运算，许多问题都是边发现、边思考、边查找、边解决的。当然，这次实验也让我熟悉了matlab的操作，初步掌握了运用数学软件进行数值分析、解决实际问题的方法，更领会了矩阵在matlab中的重要作用。今后应当在进行数学实验的过程中，发挥更多的主观能动性，积极探索更多的、潜在的问题，把问题挖得更深更透，以求更好地提高自己的科研素养与能力。四、 程序清单1. Lagrange插值法函数：function y=lagr(x0,y0,x)n=length(x0); m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0(j); end end s=p*y0(k)+s; end y(i)=s;end2. 分段线性插值法函数：function y=cutline(x0,y0,x)n=length(x0);m=length(x);for i=1:m for j=1:(n-1) if (x(i)=x0(j)&(x(i)=x0(j+1) l1=(x(i)-x0(j+1)/(x0(j)-x0(j+1); l2=(x(i)-x0(j)/(x0(j+1)-x0(j); y(i)=l1*y0(j)+l2*y0(j+1); end endend3. 自编函数的可信性检验：clear；x0=-1:0.2:1;x=-1:0.1:1;y0=x0.2;y1=lagr(x0,y0,x);y2=cutline(x0,y0,x);y3=interp1(x0,y0,x);plot(x,y1,r,x,y2,g,x,y3,+,x0,y0,k-)4. 插值主程序（共有三次插值，节点数分别是6、11、21，此处仅附节点数为21的相应程序，另外两个情况类似）：clear;x0=-5:0.5:5;y0=1./(1+x0.2);x=-5:0.05:5;y=1./(1+x.2);y1=lagr(x0,y0,x);y2=interp1(x0,y0,x); y3=interp1(x0,y0,x,cubic);y4=interp1(x0,y0,x,spline);for i=1:51 xx(i)=x(99+2*i); yy(i)=y(99+2*i); yy1(i)=y1(99+2*i); yy2(i)=y2(99+2*i); yy3(i)=y3(99+2*i); yy4(i)=y4(99+2*i);endxx;yy;yy1;yy2;yy3;yy4z=0*x;plot(x,z,k,x,y,k-,x,y1,r,x,y2,m,x,y3,g,x,y4,b)axis(-5 5 -1.5 4.05)title(节点数为21时，四种插值方法的对比)gtext(黑虚线原函数)gtext(红Lagrange)gtext(紫分段线性)gtext(绿分段三次)gtext(蓝三次样条插值)pause5. 非均匀节点时的插值主程序：clear;n=20;k=1:n;l=k.*2-1;x0=5*cos(l*pi/2/n);y0=1./(1+x0.2);x=-5:0.1:5;y=1./(1+x.2);y1=lagr(x0,y0,x);y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,cubic);y4=interp1(x0,y0,x,spline);for i=1:11 xx(i)=x(46+5*i); yy(i)=y(46+5*i); yy1(i)=y1(46+5*i); yy2(i)=y2(46+5*i); yy3(i)=y3(46+5*i); yy4(i)=y4(46+5*i);endxx;yy;yy1;yy2;yy3;yy4z=0*x;plot(x0,y0,+,x,z,k,x,y,k-,x,y1,r,x,y2,m,x,y3,g,x,y4,b)axis(-5 5 -0.05 1.05)title(非均匀节点时，四种插值方法的对比)gtext(黑虚线原函数)gtext(红Lagrange)gtext(紫分段线性)gtext(绿分段三次)gtext(蓝三次样条插值)pause6. 均匀分段的梯形求积公式：function y=fuhetixing(a,b,y0)m=length(y0);h=(b-a)/(m-1);y=h*(sum(y0)-0.5*(y0(1)+y0(m);end7. 均匀分段的辛普森求积公式：function y=fuhexinpusen(a,b,y0)m=length(y0);k=(m-1)/2;h=(b-a)/(m-1);s1=0;s2=0;for i=1:k s1=s1+y0(2*i);endfor i=1:(k-1) s2=s2+y0(2*i+1);endy=h/3*(y0(1)+y0(m)+4*s1+2*s2);end8. 检验原创积分程序与matlab自带的预设程序的吻合性：x=0:0.001:0.1;y=1/sqrt(2*pi)*exp(-0.5*x.2);for i=1:10 l=10*i+1; y1(i)=trapz(x(1:l),y(1:l); y2(i)=fuhetixing(0,x(l),y(1:l); y3(i)=quad(1/sqrt(2*pi)*exp(-0.5*x.2),0,x(l); y4(i)=fuhexinpusen(0,x(l),y(1:l);endx0=0.01:0.01:0.1;x0;y1;y2;y3;y49. 标准正态分布函数表（输出成常见格式，并保留4位小数）：clear;x=0:0.001:4;for i=1:399 l=10*i+1; y(i+1)=quad(1/sqrt(2*pi)*exp(-0.5*x.2),0,x(l);endy(1)=0;y=y+0.5;y=reshape(y,10,40);y10. 输出(x)、二三类误差和的函数：function y1,y2,y3=anothernorm(x0)if x0=