数学试验——利用MALTAB计算插值与数值积分

1、实验一插值与数值积分一、差值问题描述选择一些函数，在n个节点上用拉格朗日、分段线性、三次样条三种差值方法，计算m个插值点的函数值，通过函数和图形输出，将三种插值结果与精确值进行比较，适当增加n再作比较，由此作初步分析。方法与公式1、拉格朗日插值法上，工）-£.窜/（工）若4则，（”）=九/）=JKf）j（$一餐）-工%u）Y%一，）2、分段线性插值法3、三次样条插值法1/17S=e(x),X曰茗_"卜，=1,F?321)5.tx)=(7f.r+i.r+ctx+rf.(i=l,M)2)SG)f(f=O,L公3)S(i)£CmJ4)S"*q)=5*3)=)结

2、果与分析1、函数(1)5节点，100插值点5由图可以看出，三种差值方法中，拉格朗日插值法所得结果与真实值最接近，其次是三次样条插值法，最后是分段线性插值法。difl=0.0043;dif2=0.0196;dif3=0.0058。dif表示每种插值方法得到的插值点与实际值的差方均。(即每个插值点与实际值之差的平方的和的平均值)。由此页可以看出，拉格朗日插值法所得结果与真实值最接近，其次是三次样条插值法，最后是分段线性插值法。(2)7节点，100插值点2/17可以发现，通过增加节点，每一种插值法的结果都比原来更加接近真实值。difl=0.0015;dif2=0.0082;dif3=0.0027;(

3、3)10节点，100插值点结果进一步逼近真实值。difl=0.00062;dif2=0.0035;dif3=0.0012;(4)简要小结a.增加节点可以有效地提高插值结果的准确度。b.在本函数下，拉人格朗日插值法的效果与三次样条插值法的效果相近，两者都好于分段线性插值法。(5)对比3/17皿领igf/.atlaoS0分F殳然性插值法对比了自己编写的分段线性插值和matlab自带的分段线性插值函数。发现两者基本一致，区别无法从图中读出。difl=dif2=0.0121。误差相等。2、其它函数为了更好地对比各种插值法的优劣，仅用一种函数作测试是不够的，这里将书中提供的另外三个函数训练，结果如下：6

4、节点，100插值点此种情况与上一函数所示结果相近：拉格朗日插值法与三次样条插值法效果相近，两者优于分段线性插值法。（实验过程中发现，若节点为5,则拉格朗日插值法与三次样条插值法所得结果一致）。6节点，100插值点4/1710节点，100插值点时观察到了比较明显的Runge现象。6节点，100插值点5/1710节点，100插值点12Iiiiiii-I=交除值kvlEerrX!与1.3、总结由以上各个对比可以得出以下结论：a.拉格朗日插值法的结果与三次样条插值法的结果在许多情况下比较相像；b.除非节点取得足够多，分段线性插值法的精度才能有所保障；c.对一些有峰值的函数，是否取峰值点为节点对结果的影

5、响非常大；d.拉格朗日插值法在个别函数中会出现Runge现象，但并不是所有类似钟形的函数中都会出现，最后一个函数说明了这一点；e.在拉格朗日插值法可能出现Runge现象的情况中，三次样条插值法和分段线性插值法的效果要好于拉格朗日插值法。6/17程序清单1、拉格朗日插值法functiony=lagr(x0,y0,x)n=length(x0);m=length(x);fori=1:mz=x(i);s=0;fork=1:np=1;forj=1:nifj=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s;end2、分段线性插值法functio

6、ny=mylnterp(x0,y0,x)n=length(x0);m=length(x);fori=1:mforj=1:nifx(i)<x0(j)breakendendif(j=1)y(i)=0;elseifx0(j)=x0(j-1)y(i)=y0(j);elsey(i)=(x(i)-x0(j-1)/(x0(j)-x0(j-1)*(y0(j)-y0(j-1)+y0(j-1);endend3、“自启动”脚本x0=-1:2/5:1;y0=sqrt(1-x0.A2);7/17x=-1:0,02:1;y=sqrt(1-x,A2);y1=larg(x0,y0,x);y2=mylnterp(x0,y

7、0,x);y3=interp1(x0,y0,x,'spline');y4=interp1(x0,y0,x);dif1=0;dif2=0;dif3=0;dif4=0;fori=1:length(x)dif1=dif1+(y(i)-y1(i)A2;dif2=dif2+(y(i)-y2(i)2;dif3=dif3+(y(i)-y3(i)A2;dif4=dif4+(y(i)-y4(i)A2;enddif1=dif1/length(x);dif2=dif2/length(x);dif3=dif3/length(x);dif4=dif4/length(x);plot(x,y,'k&

8、#39;,x,y1,'r',x,y2,'b',x,y3,'g');legend('e?e?,1arg(x)','mylnterp(x)','ey'?ui?'axis(-1,1,0,1);xlabel('±?a?x');ylabel('o6y?心)；gridontext(-0.2,0.65,'ey?2?ii-?-6?eii?e?iiii?)+e?'8/17二、数值积分问题描述选择一些函数用梯形、辛普森和Gauss-Lobatto三种方法计算积分

9、。改变步长，改变精度要求，进行比较、分析方法与公式1、梯形公式梯形公式Tn=史，+丁5>+4)hl-丸T十f非等距分割梯形公式。二z八,十1Si一招)2、辛普森公式hM-1*12j_=鼠=三通+&+426阳+2士力J。=丁(5)3JM3JM2"?3、Gauss-Lobatto公式Gauss-Lobatto求积公式G"4/+£/"(乙)+4/结果与分析1、(1)解析解(2)步长为0.2,精度为时，各种积分法的结果9/172.960®673769&4531e+0002.957885697765622e+000DeltaS=L0

10、08038234920812-003-5.856741763251813c-0092.957884774745239H-0002.957886220517614e+000-3.17910848021656&G-007L7087489759832L3e-DQ7S1为自己编辑的梯形求积公式结果，s2为自己编写的Simpson求积公式结果，S3为Matlab自适应Simpson求积公式结果，S4为Matlab自适应Gauss-Lobatto求积公式结果。后两者设置的绝对误差为默认的。deltaS的各项是各种求积公示的结果与解析解的相对误差，即，s为解析解结果的近似值。从相对误差可以看出，Si

11、mpson积分公式的精度要比梯形积分公式的精度高不少。然而这里不能将前两者与Matlab自带的自适应积分公式作比较，因为两者步长不同，后者按照所要求的积分精度设置步长。不过对比后两者可以看到此种情况下，自适应Simpson求积公式结果比自适应Gauss-Lobatto求积公式结果更接近真实值。(3)改变步长及精度要求步长为0.1,精度要求为S二2.Q586310S597Q594e+0002957885715742326eW00DeltaS-2.519944859925511巴Y042.208101967794245e-01Q2.957885655644615e+0002,95788571508

12、9237e+000-2.00969S333572532e-0081.411291366429520c-014从结果可以看出，各种积分法的精度都有所提高，其中自适应积分法的精度提高得最多。步长为0.05,精度要求10/172958072055018870e+0002.957战57口365297E0002.957885714859374e+0002.957SS5715089237e400CDeltaS-6. 29976"726S15904e-005-1.25S9"30S144SSS3e-009-7.76978954918S019e-0111.411291366429520e-0

13、14前两种与直接步长相关的积分公式得到的结果的精度进一步提升。由此思考,可以自己编写“以精度为参量”的梯形积分法和Simpson积分法，即多次调用现有积分程序，判断精度是否符合要求，若不符则进一步缩小步长，知道满足要求。对比步长取0.2,比较自己编写的梯形积分公式以及Matlab自带的梯形积分公式。S1为自己所编公式的结果，S2为Matlab自带公式的结果。S=2.960S67376934531e+0002.%CS673769S4531e+000DeltaS-LOOS038234920812C-0031.QQS038234920312e-m两者结果一致。步长取0.05,再次比较2.95goi2

14、055018870+0002.958072055018871e+000DcltaS=6.29976772681&904e-0056.299767726846933e-005可以看出两者结果的最末一位不同。总体上讲，两者基本一致。(5)小结a.对于“以步长为参数”的积分公式，缩小步长可以有效地提高计算精度，同等步长情况下，对于本题函数，Simpson积分公式的结果比梯形积分公式结果精确。11/17b.对于“以精度为参数”的积分公式，自适应Simpson积分法和自适应Gauss-lobatto积分法哪种更优仅通过本函数不好予以评价。2、其它函数一解析解为ln(2)步长为0.1,精度要求为S

15、:6.931502306889305e-0016-$31471861471861e-0014.400405960962491e-0068.060684567212437e-00a6.937714031754280«'0016.93147199S629700e-001DcltaS=9005628717675必加Y042.7S4837805816601e-fl08步长为0.05,精度要求为6.933033817926941e-0016.J31471306S57081c-001DcllaS二2,253507438674981e-0041.814374322076464e-0106.

16、931473746651162c-0016.M1471861471861e-(M)12.800345675208945e-0078.060684567212437e-009解析解一步长为0.1,精度要求为S=-3.453120332461184e+001I973464947962066e4001DeltaS-1.6131195C4431303e-001l-32723217832404e-OOS-2.982425999698896H001-2.9734649084757710013.013686549992539e-003-7.23478922061489'Tl212/17步长为0.05

17、,精度要求为S=-3.003568072046070e4-001-2.97346490825647拆+001DellaS-1.012393435777952e-t)02-8.098442781228790e4)11-2.973497799783576c+001-2.973464908497234e4001L106160230733880e-005-1.4弱32263gl322)014对于这个函数，自己编写的两个公式与Matlab自带公式的差距较大。猜想原因是此函数在积分区域中变化率较大，因此须将步长缩减到一定程度才能得到较精确的值步长为0.2,精度为解析解s=(erf()-erf()/2S=9

18、.537803406031293eH)019.5449S748576686e-001DeltaS=53608527405044r-004£15547857506643氓-0089.544977324882294e-0019.544997361023964e-001-2.099126208813142e-006-i.工止步长为0.1,精度为S=9.543197962745974«-001文544997361713754e-OOIDeltaS=-L88517421弟036935004一口962如464753937Y:19.544996148317538e-001g.544997

19、3&1023964e-001-L270628250444723e-007-1.304817740556470e-0123、总结13/17a.对于相同的步长，Simpson积分公式通常比梯形积分公式更加准确（若被积函数本身为分段函数则另当别论）；b.梯形积分公式以及Simpson积分公式的好处在于可以将x,y作为参量求得最终结果，而自适应Simpson和自适应Gauss-lobatto函数则需要已知函数的解析式。程序清单1、自编梯形求积公式functions=myTrapz(x,y)n=length(x);s=0;fori=1:n-1s=s+(y(i)+y(i+1)*(x(i+1)-x(

20、i)/2;end2、自编Simpso球积公式functions=mySimpson(x,y)n=length(x);s=0;fork=1:2:n-2s=s+(x(k+2)-x(k)3*(y(k)+4*y(k+1)+y(k+2);end3、“自启动”脚本x=0:0.05:2;y=sqrt(1+x.A2);s=sqrt(5)+0.5*log(2+sqrt(5);si=myTrapz(x,y);s2=mySimpson(x,y);s3=quad('sqrt(1+x.A2),0,2,10A-8);s4=quadl('sqrt(1+x.A2),0,2,10A-8);deltaSI=(s1

21、-s)/s;deltaS2=(s2-s)/s;deltaS3=(s3-s)/s;deltaS4=(s4-s)/s;S=s1,s2,s3,s4DeltaS=deltaS1,deltaS2,deltaS3,deltaS414/17三、插值与数值积分问题描述已知一个国家地图的测量图，数据如书中P66所给的表格，试根据地图的比例及测量数据计算该国国土的近似面积。方法与公式为提高精度，在进行数值积分前需要先插值。考虑到收敛性的保障，这里选用分段线性插值和三次样条插值。由于本题没有具体的函数，因此无法使用自适应Simpson以及自适应Guass-Labatto方法进行估算，本题将使用自己编写的梯形求积公式

22、、Simpson求积公式以及Matlab自带的梯形求积公式进行计算，并进行简要对比。结果与分析1、各种公式的结果(1)分段线性插值+三种数值积分S=4.241358024647034e+0044.241356712293975e+0044.241358024682084e+004deltaS=1.125580246470337e+0031.125567122939749e+0031.125580246820842e+003第一项为自己编写的梯形求积公式所得结果，第二项为自己编写的Simpson求积公式所得结果，第三项为Matlab自带的梯形求积公式所得结果。(2)三次样条插值+三种数值积分15

23、/17S=4.246835283669813e+0044.246833971264325e+0044.246835283669541e+004deltaS=1.180352836698126e+0031.180339712643247e+0031.180352836695405e+003第一项为自己编写的梯形求积公式所得结果，第二项为自己编写的Simpson求积公式所得结果，第三项为Matlab自带的梯形求积公式所得结果。2、总结a.选择不同的插值方法做数值积分，将得到不同的积分值，且相互之间差别相对较大(4.241358和4.246835相差约0.1%);b.选择相同的插值方法不同的数值积分方法，当积分步长足够小时，所得各结果相差较小(4.246833971264325e+004和4.246835283669541e+004相差约3)程序清单x00=7.0,10.5,13.0,17.5,34,40.5,44.5,48.0,56,61,68.5,76.5,80.5,91,96,101,104,106.5,111.5,118,123.5,136.5,142,146,150,157,158;x01=158,157,150,146,142,136.5,123.5,118,111.5,106.5,104,101,96,91,80.5,76.5,68.5,61,56,4