文档介绍：四、三次样条插值
给定区间上划分,若分段函数满足:
在各个子区间,上均为的三次多项式;
在整个区间上有直至二阶的连续导数。
则称为上依次划分的三次样条函数,简称样条函数。
具体地有分段表达式:
共有个参数,它们在内节点处满足
满足样条函数定义的函数集合称为分划上的三次样条函数空间,记为,可以证明为线性空间。
若,且进一步满足插值条件
其中为节点处的给定函数值(若被插函数已知,则用代替之),则称为以为节点的三次样条函数。
其中式(3)插值节点提供了个约束条件,加上式(2)的个,合起来共有个,欲求个待定参数的唯一解,尚缺两个条件。这两个条件一般由样条函数的边界条件提供。常用三类边界条件,他们分别与三次样条函数,构成不同边界条件的样条函数插值问题。
第二类边界条件
给定边界条件两端的一阶导数值:
这相当于样条两短处的方向给定(压铁在两端点的压力方向确定),对应的插值问题如下:
对于分划,给定节点对应的函数值,以及两端点处的一阶导数值,,求三次样条函数,使
第一类边界条件
给定边界两端的二阶导数值:
这相当于在样条两端处外加一个力矩,使梁两端点处有相应的曲率。对应的插值问题如下:
对于划分,给定节点对应的函数值,以及两端点处的二阶导数值,,求三次样条函数,使
特别地,若,这相当于样条边界上不加力矩,样条在边界处是自由的,这样的样条称为自由样条,边界条件称为“自由边界条件”。
第三类边界条件
被插函数是以为周期的周期函数时,则要求也是周期函数,此时边界条件应满足:
而且还要加上。这样得出的称为周期样条插值函数。
本部分分别编写了三次样条插值在第一类边界条件(二阶导数),第二类边界条件(一阶导数),第三类边界条件(周期函数)下求解的matlab程序,并给出了实际例子进行求解。
第一类边界条件的插值求解程序
function [f,f0] =spline1 (x,y,y_1, y_N,x0)
syms t;
f = ;
f0 = ;
if(length(x) == length(y))
n = length(x);
else
disp('x和y维数不相等');
return;
end %维数检查
for i=1:n
if(x(i)<=x0)&& (x(i+1)>=x0)
index = i;
break;
end
end %找到待求解x0值所在的区间
A = diag(2*ones(1,n)); %求解m的系数矩阵
u = zeros(n-2,1);
lamda = zeros(n-1,1);
c = zeros(n,1);
for i=2:n-1
u(i-1) = (x(i)-x(i-1))/(x(i+1)-x(i-1));
lamda(i) = (x(i+1)-x(i))/(x(i+1)-x(i-1));
c(i) = 3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+ ...
3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));
A(i, i+1) = u(i-1);
A(i, i-1) = lamda(i); %形成系数矩阵及向量C
end
c(1) = 2*y_1;
c(n) = 2*y_N;
m = followup(A,c); %用追赶法求解三对角方程组
h = x(index+1) - x(index); %求x0所在区间长度
f = y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h + ...
y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h + ...
m(index)*(t-x(index))*(x(index+1)-t)^2/h/h - ...
m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h;
%求x0所在区间的插值函数
f1=simplify(f);
f=collect(f1); %化简f并合并同类项。
f0 = subs(f,'t',x0); %x0处的插值
end
其中[f,f0] =spline1 (x,y,y_1, y_N,x0)中的x,y分别是相应的插值向量,y_1和 y_N为对应的边界条件值,x0为待求插值点。当x,y的维度n不是很大的时候,我们可以选择n-1个分别位于不同且无重合区域的区间段里的点Xi,分别插值,