三次样条插值多项式matlab

三次样条插值多项式

——计算物理实验作业四

陈万物理学2013级

主程序：

clear,clc;

format rat

x = [1,4,9,16,25,36,49,64];

y = [1,2,3,4,5,6,7,8];

f1 = ;

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,c

a = 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)::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 on

end

plot(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) ); end

y(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)); end

M(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)间隔越大插值效果越差。