三次样条插值多项式matlab

《三次样条插值多项式matlab》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、三次样条插值多项式——计算物理实验作业四陈万物理学2013级13020011006l主程序：clear,clc;formatratx=[1,4,9,16,25,36,49,64];y=[1,2,3,4,5,6,7,8];f1=0.5;fn=1/16;[a,b,c,d,M,S]=spline(x,y,f1,fn);l子程序1：function[a,b,c,d,M,S]=spline(x,y,f1,fn)%三次样条插值函数%x是插值节点的横坐标%y是插值节点的纵坐标%u是插值点的横坐标%f1是左端点的一阶导数%fn是右端点的一阶导数%a是三对角矩阵对角线下边一行%b是三对

2、角矩阵对角线%c是三对角矩阵对角线上边一行%S是插值点的纵坐标n=length(x);h=zeros(1,n-1);deltay=zeros(1,n);miu=zeros(1,n-1);lamda=zeros(1,n-1);d=zeros(1,n-1);forj=1:n-1h(j)=x(j+1)-x(j);deltay(j)=y(j+1)-y(j);end%得到h矩阵forj=2:n-1sumh=h(j-1)+h(j);miu(j)=h(j-1)/sumh;lamda(j)=h(j)/sumh;d(j)=6*(deltay(j)/h(j)-(deltay(j-1)/h

3、(j-1)))/sumh;end%根据第一类边界条件，作如下规定lamda(1)=1;d(1)=6*(deltay(1)/h(1)-f1)/h(1);miu(1)=1;d(n)=6*(fn-deltay(n-1)/h(n-1))/h(n-1);%输出三对角矩阵的a,b,ca=miu;b=2*ones(1,n);c=lamda;M=chase(a,b,c,d);%调用chase函数得到Msymu;forj=1:n-1u=x(j):0.1:x(j+1);v=ones(size(u));S=(M(j)*(x(j+1)*v-u).^3/(6*h(j))+M(j+1)*(u-x

4、(j)*v).^3/(6*h(j))...+(y(j)-M(j)*h(j)^2/6)*(x(j+1)*v-u)/h(j)+(y(j+1)...-M(j+1)*h(j)^2/6)*(u-x(j)*v)/h(j));plot(u,S,'-k');holdonendplot(x,y,'-.*r');xlabel('x'),ylabel('y'),title('cubicsplineinterp');endl子程序2：functionM=chase(a,b,c,f)%追赶法求解三对角矩阵方程，Ax=f%a是对角线下边一行的元素%b是对角线元素%c是对角线上边一行的元素n=le

5、ngth(b);beta=ones(1,n-1);y=ones(1,n);M=ones(1,n);fori=(n-1):(-1):1a(i+1)=a(i);end%将a矩阵和n对应beta(1)=c(1)/b(1);fori=2:(n-1)beta(i)=c(i)/(b(i)-a(i)*beta(i-1));endy(1)=f(1)/b(1);fori=2:ny(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1));endM(n)=y(n);fori=(n-1):(-1):1M(i)=y(i)-beta(i)*M(i+1);enden

6、dl三次样条插值结果：l与拉格朗日插值作对比：分析图一知，三次样条插值结果与预期结果吻合得很好，曲线平滑连续性好，在左右短点处的小区间也吻合得很好，可以延伸到区间[a,b]外的一小段,用最邻近的小区间插值函数可以近似求得[a,b]区间外一小段范围内的函数值。分析图一和图二，三次样条插值在x(i)和x(i+1)间隔很小的地方插值效果稍差，间隔稍大时插值效果非常好。拉格朗日插值则不同，x(i)和x(i+1)间隔越大插值效果越差。