第二章 插值法_课堂演示实验.doc

湖北民族学院理学院数值分析课程教案第二章 插值法_课堂演示实验问题: 分别以函数和为例,在区间上,取 为节点,用本章所学的Lagrange、Newton、Hermite插值法以及分段线性、分段三次Hermite插值法、Spline插值法作插值计算.具体要求如下：1取，通过在同一坐标系中作出被插函数与插值函数的图形的方法，来观察插值函数的图象与被插函数的位置关系，并给出观察到的的结论.2通过计算时的插值误差,结合前面得到的直观结论,试对各种插值方法的应用作出述评.解答：一、各种插值方法的算法公式及MATLAB通用程序设已知函数表为且当时.1Lagrange插值法算法公式:, 其中余项为 , 有关.通用程序1(此程序用得较多):function yy=lagr1(x,y,xx)%用途：拉格朗日插值法求插值点xx（可以是多个）处的插值yy%格式：yy=lagr(x,y,xx), x是节点向量，y是节点对应的函数值向量，% xx是插值点（可以是多个），yy返回插值结果m=length(x);n=length(y);if m=n, error(向量x与y的长度必须一致);ends=0;for i=1:n t=ones(1,length(xx); for j=1:n if j=i t=t.*(xx-x(j)/(x(i)-x(j); end end s=s+t*y(i);endyy=s;通用程序2:function L ,C, l ,L1= lagr2(X,Y)%输入的量:n+1个节点(xi,yi)的横坐标向量X，纵坐标向量Y；%输出的量:n次拉格朗日插值多项式L及其系数向量C，基函数l及其系数矩阵L1m=length(X); L=ones(m,m);for k=1: m V=1; for i=1:m if k=i V=conv(V,poly(X(i)/(X(k)-X(i); endendL1(k,:)=V; l(k,:)=poly2sym (V);endC=Y*L1;L=Y*l;l=vpa(l,4);L=vpa(L,4);通用程序3:function y,R=lagr3(X,Y,x,M)%输入的量：X 是n+1个节点的横坐标向量，Y是纵坐标向量, x是以向量形式输入的m个插值点，M是被插函数在a,b区间上的n+1阶导数的最大值.%输出的量：y为m个插值构成的向量，R是误差限.n=length(X); m=length(x);for i=1:m z=x(i);s=0.0; for k=1:n p=1.0; q1=1.0; c1=1.0;for j=1:n if j=kp=p*(z-X(j)/(X(k)-X(j); end q1=abs(q1*(z-X(j);c1=c1*j; end s=p*Y(k)+s; end y(i)=s;endR=M*q1/c1;2Newton插值法算法公式:余项为其中 有关.通用程序1(此程序用得较多):function s=newton1(x,y,x0,nn)%Newton插值,x与y为已知的插值点及其函数值%x0为需要求的插值点向是,s返回插对应插值%nn为newton插值多项式的次数,即nn次newton插值多项式nx=length(x);ny=length(y);if nx=ny warning(向量x与y的长度应该相同) return;endm=length(x0);%按照公式,对需要求的插值点x0的每个元素进行计算for i=1:m t=0.0; %j=1; yy=y; kk=1; %求各级均差 while(kk=nn) kk=kk+1; for k=kk:nx yy(k)=(yy(k)-yy(kk-1)/(x(k)-x(kk-1); end end %求插值结果 t=yy(1); for k=2:nx u=1.0; jj=1; while(jj=X(i)&x0=X(i+1) k=i;break; end end a1=x0-X(k+1); a2=x0-X(k); a3=X(k)-X(k+1); y0=(a1/a3)2*(1-2*a2/a3)*Y(k)+(-a2/a3)2*(1+2*a1/a3)*Y(k+1)+(a1/a3)2*a2*DY(k)+(-a2/a3)2*a1*DY(k+1);通用程序2:function f,f0 = fd_hermite2(x,y,y_1,x0)syms t;f = 0.0;f0 = 0.0;if(length(x) = length(y) if(length(y) = length(y_1) n = length(x); else disp(y和y的导数的维数不相等！); return; endelse disp(x和y的维数不相等！); return;end %维数检查for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在区间h = x(index+1) - x(index); %x0所在区间长度fl = y(index)*(1+2*(t-x(index)/h)*(t-x(index+1)2/h/h + . y(index+1)*(1-2*(t-x(index+1)/h)*(t-x(index)2/h/h;fr = y_1(index)*(t-x(index)*(t-x(index+1)2/h/h + . y_1(index+1)*(t-x(index+1)*(t-x(index)2/h/h; f = fl + fr; %x0所在区间的插值函数f0 = subs(f,t,x0); %x0处的插值6Spline插值法定义 函数，且在每个小区间上是三次多项式，其中是给定节点，则称是节点上的三次样条函数。若在节点上给定函数值，并成立 ， 则称为三次样条插值函数。求三次样条函数的三弯矩法: 三次样条插值函数可以有多种表达方法，有时用二阶导数值 表示使用更方便。在力学上解释为细梁在截面处的弯矩，并且得到的弯矩与相邻两个弯矩有关，故称三弯矩方程。算法公式:在区间上是三次多项式:只要加上三种边界条件中的任一种就可导出求三弯矩的方程组,求出后代入上式即求得三次样条函数. 为便于编程，将上式其改写为：*对第一类边界条件：第一步：根据已知数据按下列公式计算各参数值，并取第二步：用追赶法解下列方程组求第三步：求分段三次样条函数的系数，第四步：计算某插值点x0处的值或写出s(x)的表达式.通用程序function y0,C=spline1(X,Y,s0,sN,x0) %X,Y是已知插值点坐标 %s0,sN是两端点的一次导数值 %x0是插值点(单个点) %y0是三次样条函数在x0处的值 %C是分段三次样条函数的系数 N=length(X); C=zeros(4,N-1); h=zeros(1,N-1); mu=zeros(1,N-1); lmt=zeros(1,N-1); d=zeros(1,N); %d表示右端函数值 h=X(1,2:N)-X(1,1:N-1); mu(1,N-1)=1; lmt(1,1)=1; mu(1,1:N-2)=h(1,1:N-2)./(h(1,1:N-2)+h(1,2:N-1); lmt(1,2:N-1)=h(1,2:N-1)./(h(1,1:N-2)+h(1,2:N-1); d(1,1)=6*(Y(1,2)-Y(1,1)/h(1,1)-s0)/h(1,1); d(1,N)=6*(sN-(Y(1,N)-Y(1,N-1)/h(1,N-1)/h(1,N-1); d(1,2:N-1)=6*(Y(1,3:N)-Y(1,2:N-1)./h(1,2:N-1)-(Y(1,2:N-1)-Y(1,1:N-2)./h(1,1:N-2)./(h(1,1:N-2)+h(1,2:N-1); %追赶法解三对角方程组 bit=zeros(1,N-1); bit(1,1)=lmt(1,1)/2; for i=2:N-1 bit(1,i)=lmt(1,i)/(2-mu(1,i-1)*bit(1,i-1); end y=zeros(1,N); y(1,1)=d(1,1)/2; for i=2:N y(1,i)=(d(1,i)-mu(1,i-1)*y(1,i-1)/(2-mu(1,i-1)*bit(1,i-1); end x=zeros(1,N); x(1,N)=y(1,N); for i=N-1:-1:1 x(1,i)=y(1,i)-bit(1,i)*x(1,i+1); end v=zeros(1,N-1); v(1,1:N-1)=(Y(1,2:N)-Y(1,1:N-1)./h(1,1:N-1); C(1,:)=Y(1,1:N-1); C(2,:)=v-h.*(2*x(1,1:N-1)+x(1,2:N)/6; C(3,:)=x(1,1:N-1)/2; C(4,:)=(x(1,2:N)-x(1,1:N-1)./(6*h); if nargin=X(1,j)& x0=X(1,j+1) omg=x0-X(1,j); y0=(C(4,j)*omg+C(3,j)*omg+C(2,j)*omg+C(1,j); end end end*对第二类边界条件：只需将第一类边界条件中的；即可(注:若称为自然边界条件)通用程序function y0,C=spline2(X,Y,s0,sN,x0) %X,Y是已知插值点坐标 %s0,sN是两端点的二阶导数值 %x0是插值点(单个点) %y0是三次样条函数在x0处的值 %C是分段三次样条函数的系数 N=length(X); C=zeros(4,N-1); h=zeros(1,N-1); mu=zeros(1,N-1); lmt=zeros(1,N-1); d=zeros(1,N); %d表示右端函数值 h=X(1,2:N)-X(1,1:N-1); mu(1,N-1)=0; lmt(1,1)=0; mu(1,1:N-2)=h(1,1:N-2)./(h(1,1:N-2)+h(1,2:N-1); lmt(1,2:N-1)=h(1,2:N-1)./(h(1,1:N-2)+h(1,2:N-1); d(1,1)=2*s0;d(1,N)=2*sN; d(1,2:N-1)=6*(Y(1,3:N)-Y(1,2:N-1)./h(1,2:N-1)-(Y(1,2:N-1)-Y(1,1:N-2)./h(1,1:N-2)./(h(1,1:N-2)+h(1,2:N-1); %追赶法解三对角方程组 bit=zeros(1,N-1); bit(1,1)=lmt(1,1)/2; for i=2:N-1 bit(1,i)=lmt(1,i)/(2-mu(1,i-1)*bit(1,i-1); end y=zeros(1,N); y(1,1)=d(1,1)/2; for i=2:N y(1,i)=(d(1,i)-mu(1,i-1)*y(1,i-1)/(2-mu(1,i-1)*bit(1,i-1); end x=zeros(1,N); x(1,N)=y(1,N); for i=N-1:-1:1 x(1,i)=y(1,i)-bit(1,i)*x(1,i+1); end v=zeros(1,N-1); v(1,1:N-1)=(Y(1,2:N)-Y(1,1:N-1)./h(1,1:N-1); C(1,:)=Y(1,1:N-1); C(2,:)=v-h.*(2*x(1,1:N-1)+x(1,2:N)/6; C(3,:)=x(1,1:N-1)/2; C(4,:)=(x(1,2:N)-x(1,1:N-1)./(6*h); if nargin=X(1,j)& x0=X(1,j+1) omg=x0-X(1,j); y0=(C(4,j)*omg+C(3,j)*omg+C(2,j)*omg+C(1,j); end end end*对第三类边界条件：，第一步：根据已知数据按下列公式计算各参数值，第二步：用LU分解法或直接在MATLAB中解下列方程组求第三步：求分段三次样条函数的系数（注：），第四步：计算某插值点x0处的值或写出s(x)的表达式.通用程序function f,f0 = spline3 (x,y,x0)%x,y是已知插值节点坐标(向量形式)%x0是插值点(单点)%f是插值点x0所在区间的插值函数%f0是插值点x0处的插值syms t;f = 0.0;f0 = 0.0;if(length(x) = length(y) n = length(x);else disp(x和y的维数不相等！); return;end %维数检查for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在区间for i=1:n if(x(i)=x0) index = i; break; endend %找到x0所在区间A = diag(2*ones(1,n-1); %求解m的系数矩阵h0 = x(2)-x(1);h1 = x(3)-x(2);he = x(n)-x(n-1);A(1,2) = h0/(h0+h1);A(1,n-1) = 1 - A(1,2);A(n-1,1) = he/(h0+he);A(n-1,n-2) = 1 - A(n-1,1);c = zeros(n-1,1);c(1) = 3* A(1,n-1)*(y(2)-y(1)/h0 + 3* A(1,2)*(y(3)-y(2)/h1;for i=2:n-2 u = (x(i)-x(i-1)/(x(i+1)-x(i-1); lamda = (x(i+1)-x(i)/(x(i+1)-x(i-1); c(i) = 3*lamda*(y(i)-y(i-1)/(x(i)-x(i-1)+ . 3*u*(y(i+1)-y(i)/(x(i+1)-x(i); A(i, i+1) = u; A(i, i-1) = lamda; %形成系数矩阵及向量cendc(n-1) = 3*( he*(y(2)-y(1)/h0+h0*( y(n)-y(n-1)/he)/(h0+he);m = zeros(n,1);m(2:n),Q,R = qrxq(A,c); %用qr分解法法求解方程组m(1) = m(n);h = x(index+1) - x(index); %x0所在区间长度f = y(index)*(2*(t-x(index)+h)*(t-x(index+1)2/h/h/h + . y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index)2/h/h/h + . m(index)*(t-x(index)*(x(index+1)-t)2/h/h - . m(index+1)*(x(index+1)-t)*(t-x(index)2/h/h; %x0所在区间的插值函数f0 = subs(f,t,x0); %x0处的插值二、对作插值计算Lagrange插值法%eg1_lagr.mclear;clf;xx=linspace(-5,5,50); y=sin(xx); %作被插函数的图象disp(n x=4.5处的插值 绝对误差的绝对值)x1=linspace(-5,5,3); y1=sin(x1); yy1=lagr1(x1,y1,xx); %2次插值chazhi_y45_2=lagr1(x1,y1,4.5);gd_wucha_limit_y45_2=abs(chazhi_y45_2-sin(4.5);disp(sprintf(%d %15.4f %15.4f,2,chazhi_y45_2,gd_wucha_limit_y45_2) x2=linspace(-5,5,5); y2=sin(x2); yy2=lagr1(x2,y2,xx); %4次插值chazhi_y45_4=lagr1(x2,y2,4.5);gd_wucha_limit_y45_4=abs(chazhi_y45_4-sin(4.5);disp(sprintf(%d %15.4f %15.4f,4,chazhi_y45_4,gd_wucha_limit_y45_4) x3=linspace(-5,5,9); y3=sin(x3); yy3=lagr1(x3,y3,xx); %8次插值chazhi_y45_8=lagr1(x3,y3,4.5);gd_wucha_limit_y45_8=abs(chazhi_y45_8-sin(4.5);disp(sprintf(%d %15.4f %15.4f,8,chazhi_y45_8,gd_wucha_limit_y45_8) plot(xx,y,m-);hold on,pause,plot(x1,y1,rs,xx,yy1,r-); hold on,pause,plot(x2,y2,b*,xx,yy2,b-); hold on,pause,plot(x3,y3,ko,xx,yy3,k-); hold on eg1_lagrn x=4.5处的插值 绝对误差的绝对值2 -0.8630 0.11454 -0.3715 0.60608 -0.9267 0.0508Newton插值法%eg1_newton.mclear;clf;xx=linspace(-5,5,50); y=sin(xx); %作被插函数的图象disp(n x=4.5处的插值 绝对误差的绝对值)x1=linspace(-5,5,3); y1=sin(x1); yy1=newton1(x1,y1,xx,2); %2次插值chazhi_y45_2=newton1(x1,y1,4.5,2);gd_wucha_limit_y45_2=abs(chazhi_y45_2-sin(4.5);disp(sprintf(%d %15.4f %15.4f,2,chazhi_y45_2,gd_wucha_limit_y45_2) x2=linspace(-5,5,5); y2=sin(x2); yy2=newton1(x2,y2,xx,4); %4次插值chazhi_y45_4=newton1(x2,y2,4.5,4);gd_wucha_limit_y45_4=abs(chazhi_y45_4-sin(4.5);disp(sprintf(%d %15.4f %15.4f,4,chazhi_y45_4,gd_wucha_limit_y45_4) x3=linspace(-5,5,9); y3=sin(x3); yy3=newton1(x3,y3,xx,8); %8次插值chazhi_y45_8=newton1(x3,y3,4.5,8);gd_wucha_limit_y45_8=abs(chazhi_y45_8-sin(4.5);disp(sprintf(%d %15.4f %15.4f,8,chazhi_y45_8,gd_wucha_limit_y45_8) plot(xx,y,m-);hold on,pause,plot(x1,y1,rs,xx,yy1,r-); hold on,pause,plot(x2,y2,b*,xx,yy2,b-); hold on,pause,plot(x3,y3,ko,xx,yy3,k-); hold on eg1_newtonn x=4.5处的插值 绝对误差的绝对值2 -0.8630 0.11454 -0.3715 0.60608 -0.9267 0.0508Hermite插值法%eg1_hermite.mclear;clf;xx=linspace(-5,5,50); y=sin(xx); %作被插函数的图象disp(n x=4.5处的插值 绝对误差的绝对值)x1=linspace(-5,5,3); y1=sin(x1); y1_1=cos(x1); yy1=hermite1(x1,y1,y1_1,xx); %5次插值chazhi_y45_2=hermite1(x1,y1,y1_1,4.5);gd_wucha_limit_y45_2=abs(chazhi_y45_2-sin(4.5);disp(sprintf(%d %15.4f %15.4f,2,chazhi_y45_2,gd_wucha_limit_y45_2) x2=linspace(-5,5,5); y2=sin(x2); y2_1=cos(x2); yy2=hermite1(x2,y2,y2_1,xx); %9次插值chazhi_y45_4=hermite1(x2,y2,y2_1,4.5);gd_wucha_limit_y45_4=abs(chazhi_y45_4-sin(4.5);disp(sprintf(%d %15.4f %15.4f,4,chazhi_y45_4,gd_wucha_limit_y45_4) x3=linspace(-5,5,9); y3=sin(x3); y3_1=cos(x3); yy3=hermite1(x3,y3,y3_1,xx); %17次插值chazhi_y45_8=hermite1(x3,y3,y3_1,4.5);gd_wucha_limit_y45_8=abs(chazhi_y45_8-sin(4.5);disp(sprintf(%d %15.4f %15.4f,8,chazhi_y45_8,gd_wucha_limit_y45_8) plot(xx,y,m-);hold on,pause,plot(x1,y1,rs,xx,yy1,r-); hold on,pause,plot(x2,y2,b*,xx,yy2,b-); hold on,pause,plot(x3,y3,ko,xx,yy3,k-); hold on eg1_hermiten x=4.5处的插值 绝对误差的绝对值2 -0.8341 0.14354 -0.9717 0.00598 -0.9775 0.0000分段线性插值法%eg1_fd_line.mclear;clf;xx=linspace(-5,5,50); y=sin(xx); %作被插函数的图象disp(n x=4.5处的插值 绝对误差的绝对值)x1=linspace(-5,5,3); y1=sin(x1); yy1=fd_line(x1,y1,xx);chazhi_y45_2=fd_line(x1,y1,4.5);gd_wucha_limit_y45_2=abs(chazhi_y45_2-sin(4.5);disp(sprintf(%d %15.4f %15.4f,2,chazhi_y45_2,gd_wucha_limit_y45_2) x2=linspace(-5,5,5); y2=sin(x2); yy2=fd_line(x2,y2,xx);chazhi_y45_4=fd_line(x2,y2,4.5);gd_wucha_limit_y45_4=abs(chazhi_y45_4-sin(4.5);disp(sprintf(%d %15.4f %15.4f,4,chazhi_y45_4,gd_wucha_limit_y45_4) x3=linspace(-5,5,9); y3=sin(x3); yy3=fd_line(x3,y3,xx); chazhi_y45_8=fd_line(x3,y3,4.5);gd_wucha_limit_y45_8=abs(chazhi_y45_8-sin(4.5);disp(sprintf(%d %15.4f %15.4f,8,chazhi_y45_8,gd_wucha_limit_y45_8) plot(xx,y,m-);hold on,pause,plot(x1,y1,rs,xx,yy1,r-); hold on,pause,plot(x2,y2,b*,xx,yy2,b-); hold on,pause,plot(x3,y3,ko,xx,yy3,k-); hold on eg1_fd_linen x=4.5处的插值 绝对误差的绝对值2 -0.8630 0.11454 -0.6474 0.33018 -0.8040 0.1736分段三次Hermite插值法%eg1_fd_hermite1.mclear;clf;xx=linspace(-5,5,50); y=sin(xx); %作被插函数的图象disp(n x=4.5处的插值 绝对误差的绝对值)x1=linspace(-5,5,3); y1=sin(x1); y1_1=cos(x1); for i=1:length(xx) yy1(i)=fd_hermite1(x1,y1,y1_1,xx(i); endchazhi_y45_2=fd_hermite1(x1,y1,y1_1,4.5); gd_wucha_limit_y45_2=abs(chazhi_y45_2-sin(4.5);disp(sprintf(%d %15.4f %15.4f,2,chazhi_y45_2,gd_wucha_limit_y45_2) x2=linspace(-5,5,5); y2=sin(x2); y2_1=cos(x2); for i=1:length(xx) yy2(i)=fd_hermite1(x2,y2,y2_1,xx(i); endchazhi_y45_4=fd_hermite1(x2,y2,y2_1,4.5);gd_wucha_limit_y45_4=abs(chazhi_y45_4-sin(4.5);disp(sprintf(%d %15.4f %15.4f,4,chazhi_y45_4,gd_wucha_limit_y45_4) x3=linspace(-5,5,9); y3=sin(x3); y3_1=cos(x3); for i=1:length(xx) yy3(i)=fd_hermite1(x3,y3,y3_1,xx(i); end chazhi_y45_8=fd_hermite1(x3,y3,y3_1,4.5);gd_wucha_limit_y45_8=abs(chazhi_y45_8-sin(4.5);disp(sprintf(%d %15.4f %15.4f,8,chazhi_y45_8,gd_wucha_limit_y45_8) plot(xx,y,m-);hold on,pause,plot(x1,y1,rs,xx,yy1,r-); hold on,pause,plot(x2,y2,b*,xx,yy2,b-); hold on,pause,plot(x3,y3,ko,xx,yy3,k-); hold on eg1_fd_hermite1n x=4.5处的插值 绝对误差的绝对值2 -1.0020 0.02444 -0.9518 0.02578 -0.9721 0.0054Spline插值法%eg1_fd_spline2.mclear;clf;xx=linspace(-5,5,50); y=sin(xx); %作被插函数的