cubic Spline Interpolation - 三次样条插值

在MATLAB上按照非结点边界条件三次样条插值，并于自带的interp1函数对比，代码如下：

%测试文件
clc
clear all; close all;

%% 测试6个点的三次样条插值
slen = 6;
x = 1 : slen;
y = rand(1,slen)*10;
xi = 1 : 0.1 : slen;
yi = interp1(x, y, xi, 'slpine');


yii = interp1Test(x, y, xi);
figure
plot(x,y,'*-', xi, yi, 'r*-', xi, yii, 'bo-');
legend('origin data','matlab interp1','cubic interolation-not a knot');
title('interpolation')

function [yi] = interp1Test(x, y, xi)
% 测试spline插值
% 输入 - x: 距离, y: 函数值 , xi 待插值点
% 输出 - yi: 待插值点的函数值

sizex = (size(x));
sizey = (size(y));
if ~(sizex(1)*sizex(2)==sizex(1) || sizex(1)*sizex(2)==sizex(2)) || ~(sizey(1)*sizey(2)==sizey(1) || sizey(1)*sizey(2)==sizey(2))
    error('the input x or y must be one dimension')
end

if length(x) ~= length(y)
    error('the lenth of x or y must be same');
end

% 输入的长度
len = length(x);  
hi = zeros(1, len - 1);
hi = x(2:end) - x(1:end-1);

mi = zeros(1, len - 1);
% 非结点边界的终止条件
A = zeros(len , len );
A(1, 1) = - hi(2); A(1,2) = hi(1)+hi(2); A(1,3) = -hi(1);
A(len , len - 2) = -hi(len-1); A(len , len - 1) = hi(len-1)+hi(len-2);  A(len , len ) = -hi(len-2);
for i = 2 : len - 1
    A(i, i-1) = hi(i-1);
    A(i, i) = 2 * (hi(i-1) + hi(i));
    A(i, i+1) = hi(i);
end

Bm = zeros(len, 0);
Bm(1,1) = 0;
Bm(len,1) = 0;
for i = 2 : len - 1
    Bm(i,1) = 6 * ((y(i+1)-y(i))/hi(i) - (y(i)-y(i-1))/hi(i-1));
end

% Ax = B
mi = A\Bm;

a = zeros(1, len-1);
b = zeros(1, len-1);
c = zeros(1, len-1);
d = zeros(1, len-1);
for i = 1 : len-1
    a(i) = y(i);
    b(i) = (y(i+1)-y(i))/hi(i) - hi(i)*mi(i)/2 - hi(i)*(mi(i+1)-mi(i))/6;
    c(i) = mi(i)/2;
    d(i) = (mi(i+1)-mi(i))/6/hi(i);
end

yi = zeros(1, length(xi));
for i = 1 : length(xi)
    xIndex = 0;
    for j = 1 : 1 : length(x)
        if xi(i) - x(j) >= 0
            xIndex = j;
        end
    end
    fprintf(2, 'xIndex = %d\n', xIndex);
    if xIndex >= 1 && xIndex < length(x)
        yi(i) = a(xIndex) + b(xIndex) * (xi(i) - x(xIndex)) + c(xIndex) * (xi(i) - x(xIndex))^2 +d(xIndex) * (xi(i) - x(xIndex))^3;
    end
    if xIndex >= length(x)
        yi(i) = a(len-1) + b(len-1) * (xi(i) - x(len-1)) + c(len-1) * (xi(i) - x(len-1))^2 +d(len-1) * (xi(i) - x(len-1))^3;
    end
     if xIndex < 1
        yi(i) = a(1) + b(1) * (xi(i) - x(1)) + c(1) * (xi(i) - x(1))^2 +d(1) * (xi(i) - x(1))^3;
    end
end


end

对比结果显示如图所示：

在interp1Test.m文件中，计算方程组Ax = B是通过A\B来计算得出的，使用了MATLAB自带的函数，也可以自行编写代码将A转化成上三角阵，求得方程组的解，如下两段代码所示。这两段代码均是求解方程组，第二段代码使用了方程组的已知信息，不考虑系数矩阵A中值为0的区域。

function [mi] = solvingEquat(A,B)
% 求解方程组
% A - [n,n]系数矩阵 B- [n,1]

if size(A,1)~=size(A,2) || size(A,1) ~= size(B,1) || size(B,2)~=1
    error('the input x or y error')
end

% 将A , B放在同一个矩阵中
C = [A, B];
len = size(A,1);

% 将C转化为上三角矩阵，且置对角线的元素值为1
for i = 1 : len-1
    % 对角线元素不能为0，如为0，移动行使对角线元素不为0
    if C(i,i) == 0
        for j = i+1 : len
            if C(j,i) ~=0
                C([j, i], :) = C([i, j], :);
            end
        end
    end 
     if C(i,i) == 0
        continue;
    end
    if i == 1
        C(1, :) =  C(1, :) / C(1, 1);
    end
    for k = i+1 : len
        C(k,:) = C(k,:) - C(i,:) * C(k,i) / C(i,i);
        if abs(C(k,i+1))>0 
            C(k,:) = C(k,:) / C(k,i+1);
        end
    end
    
end

mi = zeros(1, len);
% 直接求出mi的最后一个值
mi(len) = C(len, len + 1);
% 逆向迭代，求出mi
for i = len-1:-1:1
    for j = len : -1 : i+1
        mi(i) =  mi(i) + mi(j) * C(i, j);
    end
    mi(i) = C(i, len + 1) - mi(i);
end

end

function [mi] = solvingEquat2(A,B)
% 求解方程组
% A - [n,n]系数矩阵 B- [n,1]
% 非结点边界, 不考虑系数矩阵A中值为0的区域

if size(A,1)~=size(A,2) || size(A,1) ~= size(B,1) || size(B,2)~=1
    error('the input x or y error')
end

len = size(A,1);
if len < 4
    error('the length of A can not smaller than 4!')
end
% initial
a = zeros(1, len);
b = zeros(1, len);
c = zeros(1, len);
d = zeros(1, len);

for i =  1 : len
    if i == 1
        a(i) = A(i, i);
        b(i) = A(i, i+1);
        c(i) = A(i, i+2);
    end
    if i >= 2 && i < len
        a(i) = A(i, i-1);
        b(i) = A(i, i);
        c(i) = A(i, i+1);
    end
    if i == len
        a(i) = A(i, i-2);
        b(i) = A(i, i-1);
        c(i) = A(i, i);
    end
    d(i) = B(i);
end

% 转化为上三角矩阵
b(1) = b(1) / a(1);
c(1) = c(1) / a(1);
d(1) = d(1) / a(1);

tmp = b(2) - a(2) * b(1);
c(2) = (c(2) - a(2) * c(1)) / tmp;
d(2) =  (d(2) - a(2) * d(1))/ tmp;

for i = 3 : len-1
    tmp = b(i) - a(i) * c(i - 1);
    c(i) = (c(i)) / tmp;
    d(i) = (d(i) - d(i-1) * a(i)) / tmp;
end

% tmp = b(len) - a(len) * c(len - 2);
b(len) = b(len) - a(len) * c(len-2);
d(len) = d(len) - a(len) * d(len-2);

tmp = c(len) - b(len) * c(len-1);
c(len) = (c(len) - b(len) * c(len-1))/tmp;
d(len) = (d(len) - d(len-1) * b(len))/tmp;

% 直接求出mi的最后一个值
mi(len) = d(len);
% 逆向迭代，求出mi
for i = len -1 : -1 : 2
    mi(i) = d(i) - c(i) * mi(i+1);
end
mi(1) = (d(1) - b(1) * mi(2) - c(1) * mi(3));

end