三次样条插值多项式matlab.docxVIP

三次样条插值多项式 ——计算物理实验作业四陈万 物理学2013级 程序：clear,clc;format ratx = [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);子程序1：function [a,b,c,d,M,S]=spline(x,y,f1,fn)% 三次样条插值函数% x是插值节点的横坐标% y是插值节点的纵坐标% u是插值点的横坐标% f1是左端点的一阶导数% fn是右端点的一阶导数% a是三对角矩阵对角线下边一行% b是三对角矩阵对角线% 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); for j = 1:n-1 h(j) = x(j+1)-x(j); deltay(j) = y(j+1)-y(j);end % 得到h矩阵 for j = 2:n-1 sumh = 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(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函数得到M sym u;for j = 1:n-1 u = 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(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); hold onendplot(x,y,-.*r);xlabel(x),ylabel(y),title(cubic spline interp);end 子程序2：function M = chase(a,b,c,f)% 追赶法求解三对角矩阵方程，Ax=f% a是对角线下边一行的元素% b是对角线元素% c是对角线上边一行的元素 n = length(b);beta = ones(1,n-1);y = ones(1,n);M = ones(1,n); for i = (n-1):(-1):1 a(i+1) = a(i);end% 将a矩阵和n对应 beta(1) = c(1)/b(1);for i = 2:(n-1) beta(i) = c(i)/( b(i)-a(i)*beta(i-1) );endy(1) = f(1)/b(1);for i = 2:n y(i) = (f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1));endM(n) = y(n);for i = (n-1):(-1):1 M(i) = y(i)-beta(i)*M(i+1);end end三次样条插值结果：与拉格朗日插值作对比：分析图一知，三次样条插值结果与预期结果吻合得很好，曲线平滑连续性好，在左右短点处的小区间也吻合得很好，可以延伸到区间[a,b]外的一小段,用最邻近的小区间插值函数可以近似求得[a,b]区间外一小段范围内的函数值。分析图一和图二，三次样条插值在x(i)和x(i+1)间隔很小的地方插值效果稍差，间隔稍大时插值效果非常好。拉格朗日插值则不同，x(i)和x(i+1)间隔越大插值效果越差。