实验3 插值与数值积分

实验3 插值与数值积分 插值 可以理解为，在若干已知的函数值之间插入计算一些未知的函数值. 实验3 插值与数值积分 插值 可以理解为，在若干已知的函数值之间插入计算一些未知的函数值. 根据一组已知数据(xk,yk)(节点),认为来自于g(x),寻找一个靠近g(x)的简单函数f(x),使f(xk)=yk,并用f(x)求插值点的函数值. 三种插值方法： 拉格朗日多项式插值、 分段线性插值、 三次样条插值。 数值积分根据一组已知数据(xk,yk),认为来自于f(x),近似地求一个定积分 . 四种数值积分方法: 矩形公式、 梯形公式、 辛普森 (Simpson)公式、 高斯(Gauss)公式 教材P44----P64 参考书P219---P227

一、 三种插值方法 1 拉格朗日插值 插值多项式 根据n+1个节点得到唯一的n次多项式; 设 Ln (x)是n 次多项式： 1 拉格朗日插值 插值多项式 根据n+1个节点得到唯一的n次多项式; 设 Ln (x)是n 次多项式： 系数矩阵X为向量[x0,x1,…,xn]生成的凡德蒙阵，方程组简记作 XA=Y (4) 其中A=(an , an-1,...,a1, a0)' 为多项式系数构成的列向量。方程组 有唯一的解，即根据n+1个节点可以确定唯一的 n 次插值多项式。 （1） 对于节点(xj, yj )，应有 （2） 即多项式的系数满足方程组 f 向量A可以这样求出： A=X \ Y （3） 这样我们就可以得到插值多项式函数了。

x0=[0 2 4 5 7 8 10]; y0=[-3 -1 2 5 7 3 4]; X=vander(x0); Y=y0'; A=X\Y y=polyval(A,x); plot(x,y,'r-',x0,y0,'b*') 所求插值多项式为：p=A’ P=[ 0.0014 -0.0302 0.2151 -0.5135 0.1303 1.5119 -3.0000]

li(x)是 n 次多项式，满足 拉格朗日插值多项式 为求满足 (2)的多项式Ln(x) ，作 （5） （6） 令 计算机编程容易! （7） 显然(7)是满足(2)的多项式。由方程组(4)解的唯一性，（7）式的Ln(x)与(1)式相同。(5),(7)称为拉格朗日插值多项式，用Ln(x)计算插值称为拉格朗日插值。

为作拉格朗日插值，先作 M文件定义一个函数lagr1，保 存在搜索路径内。 function y=lagr1(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 s=p*y0(k)+s; y(i)=s; 最里层循环对一个插值点x构造li(x)，外一层循环对一个插值点构造Ln(x)，最外层循环求出每个插值点 x 处的多项式值 y。 输入节点x0,y0及插值x 后，运行 y=lagr1(x0,y0,x)， plot(x,y,’r-’,x0,y0,’b+’) 即可作出插值函数的图形

还是用前面的例子： x0=[0 2 4 5 7 8 10]; y0=[-3 -1 2 5 7 3 4]; x=0:0.1:10; y=lagr1(x0,y0,x); plot(x,y,'r-',x0,y0,'b+') Lagrange插值多项式震荡的很厉害,看19世纪range给出的一个 不收敛的反例.

o y x 作分段函数 2 分段线性插值 将每两个相邻的节点用直线连起来，如此形成的一条折线 (12) 将每两个相邻的节点用直线连起来，如此形成的一条折线 就是分段线性插值函数，记作In(x)，它满足In(xj)=yj，且In(x)在 每个小区间[xj, xj+1]上是线性函数( j=0,1,…,n)。 x y o (13) 满足In(xj)=yj 的分段线性插值函数就是： l1(x) l3(x) 1 l0(x) l5(x) 分段线性插值函数In(x)有良好的收敛性，对于x∈[a, b]有 x0 x1 x2 x3 x4 x5 分段线性插值由Matlab中的库函数interp1(x0,y0,x)实现 。

三次样条函数记作S(x)，a≤x≤b，要求满足以下条件： a) 在每个小区间[xi-1, xi ](i=1,…,n)上是3次多项式； 3 三次样条插值 三次样条函数记作S(x)，a≤x≤b，要求满足以下条件： a) 在每个小区间[xi-1, xi ](i=1,…,n)上是3次多项式； b) 在a≤x≤b上二阶导数连续； (14) c) S(xi )=yi , i=0,1,…,n. 共有(3n-3)+(n+1)=4n-2个条件确定4n个待定系数，尚需要补充两 个条件，常用的是自然边界条件： S(x)={Si(x), x∈[xi-1, xi ](i=1,…,n)} Si(x)=ai x3+bi x2+ci x+di (15) 其中ai , bi , ci , di 为待定系数，共4n个。 (17) 由条件b)，应有 在不太苛刻的条件下，也有三次样条插值函数是收敛的，即 （16） 三次样条插值可以用库函数spline(x0,y0,x)实现。

对三种插值作出比较： n=11;m=21; %n为节点数，m为插值点数 x=-5:10/(m-1):5; %对x从-5到5按m分度 y=1./(1+x.^2); %定义函数y，计算出分度点值 z=0*x; %x轴 x0=-5:10/(n-1):5; %节点横坐标 y0=1./(1+x0.^2); %节点纵坐标 y1=lagr1(x0,y0,x); %拉格朗日插值 y2=interp1(x0,y0,x); %分段线性插值 y3=interp1(x0,y0,x,‘spline’); %三次样条插值另一种方式 result=[x’ y’ y1’ y2’ y3’] %输出数值结果 plot(x,z,‘k’,x,y,‘y-’,x,y1,‘b-’,x,y2,‘g:’,x,y3,‘r-’) %绘图 legend('x-Axis','y=1/(1+x^2)','Lagr','PieceLinear','Spline')

对数值结果 作出比较 result = -5.0000 0.0385 0.0385 0.0385 0.0385 x y y1 y2 y3 -5.0000 0.0385 0.0385 0.0385 0.0385 -4.5000 0.0471 1.5787 0.0486 0.0484 -4.0000 0.0588 0.0588 0.0588 0.0588 -3.5000 0.0755 -0.2262 0.0794 0.0745 -3.0000 0.1000 0.1000 0.1000 0.1000 -2.5000 0.1379 0.2538 0.1500 0.1401 -2.0000 0.2000 0.2000 0.2000 0.2000 -1.5000 0.3077 0.2353 0.3500 0.2973 -1.0000 0.5000 0.5000 0.5000 0.5000 -0.5000 0.8000 0.8434 0.7500 0.8205 0.0000 1.0000 1.0000 1.0000 1.0000 0.5000 0.8000 0.8434 0.7500 0.8205 1.0000 0.5000 0.5000 0.5000 0.5000 1.5000 0.3077 0.2353 0.3500 0.2973 2.0000 0.2000 0.2000 0.2000 0.2000 2.5000 0.1379 0.2538 0.1500 0.1401 3.0000 0.1000 0.1000 0.1000 0.1000 3.5000 0.0755 -0.2262 0.0794 0.0745 4.0000 0.0588 0.0588 0.0588 0.0588 4.5000 0.0471 1.5787 0.0486 0.0484 5.0000 0.0385 0.0385 0.0385 0.0385 x y y1 y2 y3 对数值结果 作出比较

下面是三种插值的图形比较 通过三种插值与精确值比较，节点处相等，Spline插值最好。

三种插值的比较 拉格朗日插值多项式Ln(x)近似g(x)，节点增加，Ln(x)的次 数n变大， Ln(x)的光滑性变坏，甚至出现振荡，当n→∞时并 不能保证Ln(x)处处收敛于g(x)。主要用于理论分析，实际意义 不大。 分段线性插值函数In(x)有良好的收敛性，但光滑性较差。 n 越大，分段越多，插值误差越小。数学、物理中的函数表，数 理统计中用的概率分布表等，用分段线性插值就够了。 三次样条插值函数S(x)光滑性好，也有良好的收敛性，是一 种非常适合要求得到光滑曲线的插值方法。计算量相对较大。但误差估计较困难.

二、数控机床加工零件 y 问题：右图为待加工零件外形，其上部分点坐标如书P45表3.1给出,数控机床加工时根据工艺的要求，需对已知数据进行加密，从面得到步长很小的（x，y）坐标。 x 因图形关于Y轴对称，故只须对右半部进行插值，将右半部分逆时针旋转90O时变为上半部分,此时作变换 v 5 -5 5 u u

subplot(1,3,1),plot(x0,y0),axis([0 5 -5 5]) 0.86 0.74 0.64 0.57 0.50 0.44 0.40 0.36 0.32 ... 0.29 0.26 0.24 0.20 0.15 0 -1.40 -1.96 -2.37 ... -2.71 -3.00 -3.25 -3.47 -3.67 -3.84 -4 -4.14 ... -4.27 -4.39 -4.49 -4.58 -4.66 -4.74 -4.80 ... -4.85 -4.90 -4.94 -4.96 -4.98 -4.99 -5.0]; u0=-y0;v0=x0; u=-5:0.05:5; v1=interp1(u0,v0,u); v2=spline(u0,v0,u); [v1' v2' -u'] subplot(1,3,1),plot(x0,y0),axis([0 5 -5 5]) subplot(1,3,2);plot(v1,-u);axis([0 5 -5 5]) subplot(1,3,3);plot(v2,-u);axis([0 5 -5 5]) figure plot(x0,y0,v1,-u,v2,-u); axis equal

上机作业 P64 1、学号1~20做1（1）（3）； 21~ 做1（2）（4）； 2、