第 2 页 共 9 页
function s=csfit(x,y,dx0,dxn) %编写三次样条插值函数文件 n=length(x)-1; h=diff(x); d=diff(y)./h; a=h(2:n-1); b=2*(h(1:n-1)+h(2:n)); c=h(2:n); u=6*diff(d); b(1)=b(1)-h(1)/2; u(1)=u(1)-3*(d(1)); b(n-1)=b(n-1)-h(n)/2; u(n-1)=u(n-1)-3*(-d(n)); for k=2:n-1 temp=a(k-1)/b(k-1); b(k)=b(k)-temp*c(k-1); u(k)=u(k)-temp*u(k-1); end m(n)=u(n-1)/b(n-1); for k=n-2:-1:1 m(k+1)=(u(k)-c(k)*m(k+2))/b(k); end m(1)=3*(d(1)-dx0)/h(1)-m(2)/2; m(n+1)=3*(dxn-d(n))/h(n)-m(n)/2; for k=0:n-1 s(k+1,1)=(m(k+2)-m(k+1))/(6*h(k+1)); s(k+1,2)=m(k+1)/2; s(k+1,3)=d(k+1)-h(k+1)*(2*m(k+1)+m(k+2))/6; s(k+1,4)=y(k+1); End %主函数文件 clear all clc x=-5:0.01:5; z=5./(1+x.^2); plot(x,z,'gx') %作原函数图像
function [C,D]=newpoly(x,y) n=length(x); D=zeros(n,n); D(:,1)=y'; for j=2:n for k=j:n
%编写牛顿插值函数文件
D(k,j)=(D(k,j-1)-D(k-1,j-1))./(x(k)-x(k-j+1)); end end C=D(n,n); for k=(n-1):-1:1 C=conv(C,poly(x(k))); m=length(C); C(m)=C(m)+D(k,k); end
的图形比较。 步骤：打开 matlab 软件，编写如下代码： function [C,L]=lagran(x,y) %编写 lagran 插值函数文件 w=length(x); n=w-1; L=zeros(w,w); for k=1:n+1 V=1; for j=1:n+1 if k~=j V=conv(V,poly(x(j))/(x(k)-x(j))); end end L(k,:)=V; end C=y*L
二、
实验环境及相关情况（包含使用软件、实验设备、主要仪器及材料等）
装有 matlab 软件的计算机一台
第 1 页 共 9 页
三、
实验内容及步骤（包含简要的实验步骤流程）
5 在每个结点 x k 的值， 做出插值函数的图形并与 y f ( x) 1 x2
实验内容： 1.将区间[-5， 5]作 10 等分， 计算函数 y
第 3 页 共 9 页
x1=-5:1:5; y=5./(1+x1.^2); [C,L]=lagran(x1,y); xx=-5:0.1:5; yy=polyval(C,xx); hold on plot(xx,yy,'k*',x1,y,'o')
%描绘 lagran 函数插值图像以及插值点
[C,D]=newpoly(x1,y) x2=-5:0.01:5; y2=polyval(C,x2); plot(x2,y2,'r:') %作牛顿插值图像 x0=-5:0.05:5; y1=interp1(x1,y,x0,'linear');%求分段线性插值函数在 x0 上的值 plot(x0,y1,'-'); grid on x=-5:1:5; y=5./(1+x.^2); dx0=0.0739645; dxn=-0.0739645; S=csfit(x,y,dx0,dxn) x1=-5:0.01:-4;y1=polyval(S(1,:),x1-x(1)); x2=-4:0.01:-3;y2=polyval(S(2,:),x2-x(2)); x3=-3:0.01:-2;y3=polyval(S(3,:),x3-x(3)); x4=-2:0.01:-1;y4=polyval(S(4,:),x4-x(4)); x5=-1:0.01:-0;y5=polyval(S(5,:),x5-x(5)); x6=0:0.01:1;y6=polyval(S(6,:),x6-x(6)); x7=1:0.01:2;y7=polyval(S(7,:),x7-x(7)); x8=2:0.01:3;y8=polyval(S(8,:),x8-x(8)); x9=3:0.01:4;y9=polyval(S(9,:),x9-x(9)); x10=4:0.01:5;y10=polyval(S(10,:),x10-x(10)); plot(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5,x6,y6,x7,y7,x8,y8,x9,y9,x10,y10,x,y,'.') %作三次样条插值图像 grid on 2.估计某地居民的用水速度和每天的总用水量. function [a,b]=csd(X1,Y1) %最小二乘法 xmean=mean(X1) ymean=mean(Y1) sumx2=(X1-xmean)*(X1-xmean)'; sumxy=(Y1-ymean)*(X1-xmean)'; a=sumxy/sumx2 b=ymean-a*xmean
左图是居民用水的水位变化函数， 根据点的排列容易知道：它的每一小 段呈线性分布，而对应的直线就是通 过最小二乘法拟合得得到的
第 8 页 共 9 页
上图就是经过改善的飞机的外形轮廓
五、
实验总结（包括心得体会、问题回答及实验改进意见，可附页）
1.通过以上比较，三次样条插值的效果是最好的，而 lagran 和牛顿插值是不稳定的 2.影响误差的大小往往与区间内取得结点数有关 3.对于飞机外形的轮廓，进过二乘法拟合得到的还是不光滑的，需要细分节点，或再次对尾部进行拟合， 这样得到的线条就会变得光滑 4.居民的水位呈一个周期函数，水位变化总是从一个最高到最低值
3 9
39 78 156 280 520]; 28 35 36 30 0];
x1=[3 39 78 156 280 520 ]; y1=[9 28 35 36 30 0]; x2=[520 280 156 78 39 3]; y2=[0 -30 -36 -35 -28 -9]; xx = 0:0.001:520; yy = spline(x1,y1,xx); %三次样条 yy1 = spline(x2,y2,xx); %三次样条 y2=[-30 0 30];x2=[280 520 280]; %进一步对尾端进行三次样条使其变光滑 x4=-30:0.001:30; x3 = spline(y2,x2,x4); plot(xx,yy,xx,yy1,X,Y,'r',x3,x4,'k') %做出图像
第 5 页 共 9 页
3.已知某个直升飞机旋转记忆外形轮廓线 12 个点的坐标：
X Y 520 0 280 -30 156 -36 78 -35 39 -28 3 -9 0 0 3 9 39 28 78 35 156 36 280 30 520 0
采用三次样条技术，画出飞机外形轮廓线。 clc clear all X=[520 280 156 78 39 3 0 Y=[0 -30 -36 -35 -28 -9 0
课程名实验 （综合性实验） 及实验名称 姓 名 学 号 系 班 别 级
实验地点 指导教师 一、
孙丽英
实验日期 同组其他成员
实验时数 成 绩
实验目的及要求
实验目的： 1.观察拉格朗日插值的龙格（Runge）现象，探索避免此现象发生的方法，比较不同方法的插值效果。 2.学会用最小二乘法求拟合数据的多项式，并应用算法于实际问题。 3.体会三次样条技术在实现图像光滑度方面的优越性。 实验要求： 1.（1）对函数作拉格朗日插值。在 MATLAB 中用内部函数 plot 利用插值点绘制函数的图形。 （2）对函数作 Newton 插值。在 MATLAB 中用内部函数 plot 利用插值点绘制函数的图形。 （3）对函数作分段线性插值。在 MATLAB 中用内部函数 plot 利用插值点绘制函数的图形。 （4）对函数作三次样条插值。在 MATLAB 中用内部函数 plot 利用插值点绘制函数的图形。 （5）在 MATLAB 中用内部函数 ezplot 直接绘制函数的图形，并与以上方法做出的插值函数的图形进行 比较（自编程序，用不同颜色、不同结点符号将（1）-(5)的结果画在一张图上） 。 2.(1) 用最小二乘法求以上数据的拟合多项式 f (t ) ，并做出 f (t ) 的图形。 （2）根据题目要求，估计一天的总用水量. 3.（1）用 MATLAB 的内部函数 plot 直接画出上述数据点 ( X , Y ) 的图形。 （2）利用 MATLAB 软件求出以上数据的三次样条插值多项式，并画出图形。 （3）自编程序将以上结果画在一张图上，比较其差异，给出你的结论。
第 4 页 共 9 页
clc clear all %对每段时间内水位 f(t)进行最小二乘法拟合 X1=[0 3316 6635 10619 13937 17921 21240 25202 28543 32284]/3600; Y1=[3175 3110 3054 2994 2947 2892 2850 2795 2752 2697]*0.01; [a,b]=csd(X1,Y1); t=0:0.5:10; y1=a*t+b; t1=[35932 39332]/3600; s1=[0 0]; grid on X2=[39435 43318 46636 49953 53936 57254 60574 64554 68535 71854 75021]/3600; Y2=[3550 3445 3350 3260 3167 3087 3012 2927 2842 2767 2697]*0.05; t2=[79254 82649]/3600; s2=[0 0]; [a,b]=csd(X2,Y2); d2=11:0.5:21; y2=a*d2+b; X3=[85968 89953 93270]/3600; Y3=[3475 3397 3340]*0.05; d3=24:0.5:26; [a,b]=csd(X3,Y3); y3=a*d3+b; plot(X1,Y1,'.',X2,Y2,'.',X3,Y3,'.',t,y1,d2,y2,d3,y3,t1,s1,t2,s2)