三次样条插值多项式matlab

1、三次样条插值多项式计算物理实验作业四物理学 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

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);for j = 1:n-1h(j) = x(j+1)-x(j);deltay(j) = y(j+1)-y(j);end % 得到 h 矩阵for j = 2:n-1sumh = h(j-1) + h(j);miu(j) = h(j-1) / sumh;lamda(j) = h(j) / sumh;d(j) =

3、 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 函数得到 Msym u;for j = 1:n-1u = x(j):0.1:x(j+1);v = ones(size(u

4、);S = (M(j)*(x(j+1)*v-u)A3/(6*h(j)+M(j+1)*(u-x(j)*v)A3/(6*h(j) +(y(j)-M(j)*h(j)A2/6)*(x(j+1)*v-u)/h(j)+(y(j+1).-M(j+1)*h(j)A2/6)*(u-x(j)*v)/h(j);plot(uS'-k');hold onendplot(x,y,'-.*r');xlabel('x'),ylabel('y'),title('cubic spline interp');end子程序2:function M = c

5、hase(a,b,c,f)%追赶法求解三对角矩阵方程，Ax=f% a是对角线下边一行的元素%b是对角线元素%c是对角线上边一行的元素n = length(b);欢迎下载4beta = ones(1,n-1);y = ones(1,n);M = ones(1,n);for i = (n-1):(-1):1a(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:ny(i) = (f(i)-a(i)*y

6、(i-1)/(b(i)-a(i)*beta(i-1); endM(n) = y(n);for i = (n-1):(-1):1M(i) = y(i)-beta(i)*M(i+1);endend三次样条插值结果：10203。40505070Vew Inicrrt Tools Deaktep Window Help住势赛露）目口回 与拉格朗日插值作对比：Q Figure 1-LXFi c Edit View niert oslsWindow Help，臼4 ：瑞| b W势昵* - a 西值100豺60j 4020102。304050 GO 70分析图一知，三次样条插值结果与预期结果吻合得很好，曲线平滑连续性好，在左右短点处的小区间也吻合得很好，可以延伸到区间a,b外的一小段，用最邻近的小区间插值函数可以近似求得a,b区间外