大学数学实验三_插值与数值积分.pdf.pdf

1 大学数学实验三 插值与数值积分 实验报告实验报告 【实验目的】 1、掌握用 matlab 计算拉格朗日、分段线性、三次样条三种插值方法，改 变节点的数目，对三种插值结果进行初步分析。 2、掌握用 matlab 及梯形公式、辛普森公式计算数值积分。 3、通过实例学习用插值和数值积分解决实际问题。 【实验内容】 1.2 对函数对函数 在在 n 个（个（511）节点上，用拉格朗日、分段线性、三次样条三种）节点上，用拉格朗日、分段线性、三次样条三种 插值方法，计算插值方法，计算 m 个（个（50100）插值点的函数值。通过数值和图形输出，将三种插值结果）插值点的函数值。通过数值和图形输出，将三种插值结果 与精确值进行比较。适当增加与精确值进行比较。适当增加 n 再作比较，由此初步分析。再作比较，由此初步分析。 本题要求在 matlab 中，用三种插值方法计算 m 个插值点的函数值，并输出结果数据和图 线，通过比较判断三种插值方法的误差。 编写实现拉格朗日插值的函数lagr(x0,y0,x) function y=lagr(x0,y0,x) n=length(x0); m=length(x); for i=1:m z=x(i); s=0; for k=1:n p=1; 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; end 2 编写实现分段线性插值的函数xianxing(x0,y0,x) function y=xianxing(x0,y0,x) n=length(x0); m=length(x); for i=1:m z=x(i); s=0; for j=2:n if (x0(j-1)=z) elseif (x0(j)=z) else k=0; end s=k*y0(j)+s; end y(i)=s; end end 源程序如下。 n=11;m=51; x=-1:2/(m-1):1; y=(1-x.2).0.5; x0=-1:2/(n-1):1; y0=(1-x0.2).0.5; y1=lagr(x0,y0,x); y2=xianxing(x0,y0,x); y3=spline(x0,y0,x); hold on; grid on; plot(x,y,g); gtext(original,fontsize,12); plot(x,y1,r); gtext(lagr,fontsize,12); plot(x,y3,b); gtext(spline,fontsize,12); plot(x,y2,k); gtext(lin(fun.),fontsize,12); legend(精确,拉格朗日插值,三次样条插值,分段线性插值); hold off; s= x y y1 y2 y3 x(1:2:51) y(1:2:51) y1(1:2:51) y2(1:2:51) y3(1:2:51) 3 经验证，利用自编的分段线性插值函数xianxing(x0,y0,x)得到的结果与用 matlab 自 带的interp1(x0,y0,x)得到的结果一致。 下图为 n=11, m=51 时的输出。 从上图可以看出，n=11, m=51 时，拉格朗日插值法的曲线最接近精确值曲线，三次样条插 值法次之，分段线性插值法误差最大。 输出的数据如下。 x y y1 y2 y3 -1.0000 0 0 0 0 -0.9200 0.3919 0.3498 0.2400 0.3131 -0.8400 0.5426 0.5368 0.4800 0.5245 -0.7600 0.6499 0.6518 0.6400 0.6595 -0.6800 0.7332 0.7344 0.7200 0.7429 -0.6000 0.8000 0.8000 0.8000 0.8000 -0.5200 0.8542 0.8539 0.8466 0.8510 -0.4400 0.8980 0.8979 0.8932 0.8968 -0.3600 0.9330 0.9330 0.9292 0.9336 -0.2800 0.9600 0.9600 0.9545 0.9607 -0.2000 0.9798 0.9798 0.9798 0.9798 -0.1200 0.9928 0.9928 0.9879 0.9926 -0.0400 0.9992 0.9992 0.9960 0.9992 0.0400 0.9992 0.9992 0.9960 0.9992 0.1200 0.9928 0.9928 0.9879 0.9926 0.2000 0.9798 0.9798 0.9798 0.9798 0.2800 0.9600 0.9600 0.9545 0.9607 4 0.3600 0.9330 0.9330 0.9292 0.9336 0.4400 0.8980 0.8979 0.8932 0.8968 0.5200 0.8542 0.8539 0.8466 0.8510 0.6000 0.8000 0.8000 0.8000 0.8000 0.6800 0.7332 0.7344 0.7200 0.7429 0.7600 0.6499 0.6518 0.6400 0.6595 0.8400 0.5426 0.5368 0.4800 0.5245 0.9200 0.3919 0.3498 0.2400 0.3131 1.0000 0 0 0 0 从数据表中也可看出， 拉格朗日插值法误差最小， 在-0.36,0.36区间上， 几乎与原函数无异， 分段线性插值法的误差最大。 上图左为 n=21, m=51 时的输出结果，上图右为 n=31, m=51 时的输出结果。比较上两图和前 面 n=11, m=51 时输出的结果，可以看出，随着 n 增大，各种插值方法的误差都会减小。但 是各种插值方法与精确值的误差仍然是拉格朗日方法最小、 分段线性插值方法最大。 从上右 图中可以看出， 此时拉格朗日插值的曲线与精确值曲线已经几乎重合， 可见对于本函数利用 拉格朗日插值法最合适。 5 4.3 对函数对函数 用梯形、辛普森和高斯三种方法计算积分。改变步长（对梯形） ，改用梯形、辛普森和高斯三种方法计算积分。改变步长（对梯形） ，改 变精度要求（对后两种方法） ，进行比较、分析。变精度要求（对后两种方法） ，进行比较、分析。 求精确值。 计算公式（形式上根据编程需要或与书上不同） ： （ 为偶数） 编写用梯形方法计算积分的程序。 n=2; x=0:(2/n):2; y=(1+x.2).0.5; z=(5.0.5)+0.5*log(2+(5.0.5); %精确值 s1=0; s2=0; for i=0:(2/n):(2-(2/n) a=i; b=(1+a.2).0.5; s1=(2/n)*b +s1; end for j=(2/n):(2/n):2 c=j; d=(1+c.2).0.5; s2=(2/n)*d +s2; end s=0.5*(s1+s2); %梯形面积 w=(s-z)/z; %相对误差 z,s,w 6 改变 n 的值得到 7 组数据： n=2 n=5 n=8 n=10 n=20 n=30 n=100 z 2.9579 2.9579 2.9579 2.9579 2.9579 2.9579 2.9579 s 3.0322 2.9698 2.9625 2.9609 2.9586 2.9582 2.9579 |w| 0.0251 0.0040 0.0016 0.0010 0.0003 0.0001 1.0080e-05 经验证，利用上述自编程序得到的结果与用 matlab 自带的trapz(x)得到的结果完全一 致。 编写用辛普森公式计算积分的程序。 n=2; %必须输入一个偶数 x=0:(2/n):2; y=(1+x.2).0.5; z=(5.0.5)+0.5*log(2+(5.0.5); %精确值 s=0; for i=0:(4/n):(2-(4/n) a=i; b1=(1+a.2).0.5; b2=(1+(a+2/n).2).0.5; b3=(1+(a+4/n).2).0.5; s=(2/3)/n)*(b1+4*b2+b3)+s; end %辛普森求面积 w=(s-z)/z; %相对误差 z,s,w 改变 n 的值得到 3 组数据： n=2 n=4 n=6 z 2.9579 2.9579 2.9579 s 2.9643 2.9580 2.9579 |w| 0.0022 2.3627e-05 1.5884e-06 从上述数据可以看出，随着 n 增大，w 越来越小，得到的积分值越来越接近精确值。比较梯 形公式和辛普森公式两种方法，我们发现辛普森公式的精度更高，因为 n=6 时相对误差已 经小于 0.0002%，而用梯形公式要到 n=100 时相对误差仍大于 0.001%。 编写用常用的高斯公式计算积分的程序。 n=1; x=0:(2/n):2; y=(1+x.2).0.5; z=(5.0.5)+0.5*log(2+(5.0.5); %精确值 s=0; for i=(2/n):(2/n):2 7 a=i; z1=0.5*(2*i-2/n)+(2/n)/2/sqrt(3); z2=0.5*(2*i-2/n)-(2/n)/2/sqrt(3); b1=(1+z1.2).0.5; b2=(1+z2.2).0.5; s=(1/n)*(b1+b2)+s; end %常用的高斯公式求面积 w=(s-z)/z; %相对误差 z,s,w 改变 n 的值得到 3 组数据： n=1 n=2 n=3 z 2.9579 2.9579 2.9579 s 2.9533 2.9578 2.9579 |w| 0.0016 1.7729e-05 1.0102e-06 可见，常用的高斯公式精度略高于辛普森公式。 用自适应辛普森公式和自适应高斯公式求积分： z=(5.0.5)+0.5*log(2+(5.0.5); s1=quad(1+x.2).0.5,0,2); %自适应辛普森公式 s2=quadl(1+x.2).0.5,0,2); %自适应高斯公式 w1=(s1-z)/z; w2=(s2-z)/z; s1,w1,s2,w2 输出结果为 2.957885697765622，-5.8567e-09，2.957886220517614，1.7087e-07 比较相对误差可看出，自适应辛普森公式精度较高。 将精度要求提高到 时，输出结果为 2.957885715089356，5.4350e-14，2.957885715089237，1.4113e-14 可见，随着精度要求的提高，输出的结果越来越接近精确值，但有时我们可能不需要如此精 确的结果，要根据实际的需要确定精度要求。 11 为了计算欧洲某国的国土面积，首先对地图作如下测量：以由西向东的方向为为了计算欧洲某国的国土面积，首先对地图作如下测量：以由西向东的方向为 x 轴，由轴，由 南到北的方向为南到北的方向为 y 轴轴， 选择方便的原点， 并将从最西边界点到最东边界点在， 选择方便的原点， 并将从最西边界点到最东边界点在 x 轴上的区间适轴上的区间适 当地划分为若干段，在每个分点的当地划分为若干段，在每个分点的 y 方向测出南北边界点的纵坐标方向测出南北边界点的纵坐标 y1和和 y2，得到下表数据，得到下表数据 （单位：（单位：mm） 。已知） 。已知 18mm 相当于相当于 40km，试计算该国国土的近似面积， 并与精确值（，试计算该国国土的近似面积， 并与精确值（41288 平方公里）比较。平方公里）比较。 8 x 7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 y1 44 45 47 50 50 38 30 30 34 y2 44 59 70 72 93 100 110 110 110 x 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 y1 36 34 41 45 46 43 37 33 28 y2 117 118 116 118 118 121 124 121 121 x 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0 y1 32 65 55 54 52 50 66 66 68 y2 121 122 116 83 81 82 86 85 68 本题可利用 matlab 自带的梯形公式 trapz(x,y)来计算国土面积。 编写的计算国土面积的程序如下。 clear all; org=7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0 44 45 47 50 50 38 30 30 34 36 34 41 45 46 43 37 33 28 32 65 55 54 52 50 66 66 68 44 59 70 72 93 100 110 110 110 117 118 116 118 118 121 124 121 121 121 122 116 83 81 82 86 85 68 hold on title(original outline,fontsize,14); plot(org(1,:),org(2,:),b,org(1,:),org(3,:),g); hold off z1=trapz(org(1,:),org(2,:); z2=trapz(org(1,:),org(3,:); z=(z2-z1)/18/18*40*40; z 输出图形为 9 输出国土面积计算结果为：4.2414e+004 与精确值相比，相对误差为 2.73% 亦可用如下仿照课堂例题编写的程序实现面积计算： clear all; x0=7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0 158.0 157.0 150.0 146.0 142.0 136.5 123.5 118.0 111.5 106.5 104.0 10