实验3 插值与数值积分.

Presentation on theme: "实验3 插值与数值积分."— Presentation transcript:

1 实验3 插值与数值积分

2 实验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

3 一、 三种插值方法 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） 这样我们就可以得到插值多项式函数了。

4 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=[ ]

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

6 为作拉格朗日插值，先作 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+’) 即可作出插值函数的图形

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

8 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)实现 。

9 三次样条函数记作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)实现。

10 对三种插值作出比较： 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')

11 对数值结果 作出比较 result = -5.0000 0.0385 0.0385 0.0385 0.0385 x y y1 y2 y3
x y y y y3 对数值结果 作出比较

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

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

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

15 subplot(1,3,1),plot(x0,y0),axis([0 5 -5 5])
]; 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([ ]) subplot(1,3,2);plot(v1,-u);axis([ ]) subplot(1,3,3);plot(v2,-u);axis([ ]) figure plot(x0,y0,v1,-u,v2,-u); axis equal

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