插值法——三次樣條插值

　　1.三次樣條插值函數

%%三次樣條插值
%%bc為boundary conditions(邊界條件)，當已知兩端點的一階導數值時為-1，當已知兩端的二階導數時為0，當函數為周期函數時為1
%%X為節點值，Y為函數表達式（attribute=0）或者具體值（attribute=1）
function CSI = Cubic_spline_interpolation(X,Y,precision,attribute,bc)
[m,n] = size(X);a = min(X);b = max(X);n = n-1;
X = sort(X);

%%畫已知函數或已知點圖像
if attribute == 0
    F = subs(Y,X);
    t = a:(b-a)/precision:b;
    Y_real = subs(Y,t);
    pic = figure;
    set(pic,'color','w');
    plot(X,F,'r*',t,Y_real,'b');
    grid on
    xlabel('x shaft');ylabel('y shaft');
    title('三次樣條插值');
    hold on
elseif attribute ==1
    F = Y;
    pic = figure;
    set(pic,'color','w');
    plot(X,F,'r*');
    grid on
    xlabel('x shaft');ylabel('y shaft');
    title('三次樣條插值');
    hold on
end

if bc == -1
    left_endpoint = input('輸入左端點的一階導數值：');
    right_endpoint = input('輸入右端點的一階導數值：');
    for i = 1:n
        X_one_interval{i} = [X(i),X(i+1)];%%節點構造區間
        F_one_interval{i} = [F(i),F(i+1)];%%構造節點值區間
        F_all_bad(i) = The_first_order_All_bad(F_one_interval{i},X_one_interval{i},attribute);%%節點區間的一階均差
        h(i) = X(i+1)-X(i);
    end
    %%樣條插值算法
    for i = 1:n-1
        mu(i) = h(i)/(h(i)+h(i+1));
    end
    mu(n) = 1;
    lambda(1) = 1;
    for i = 2:n
        lambda(i) = h(i)/(h(i-1)+h(i));
    end
    d(1)=6*(F_all_bad(1)-left_endpoint)/h(1);
    for i = 2:n
        d(i) = 6*(F_all_bad(i)-F_all_bad(i-1))/(h(i-1)+h(i));
    end
    d(n+1) = 6*(right_endpoint-F_all_bad(n))/h(n);
    A = zeros(n+1,n+1);
    for i = 1:n+1
        for j = 1:n+1
            if i == j
                A(i,j) = 2;
            elseif (i-j) == 1
                A(i,j) = mu(j);
            elseif (j-i) == 1
                A(i,j) = lambda(i);
            end
        end
    end
    %%解樣條算法中的線性方程組得出解
    disp('樣條函數初始系數：');
    M = vpa((A\d')',6)
    %%輸出節點區間及對應的樣條函數
    syms x;
    for i = 1:n
        CSI{1,i} = X_one_interval{i};
        CSI{2,i} = vpa((X(i+1)-x)^3*M(i)/(6*h(i))+(x-X(i))^3*M(i+1)/(6*h(i))+(F(i)-M(i)*h(i)^2/6)*(X(i+1)-x)/h(i)+(F(i+1)-M(i+1)*h(i)^2/6)*(x-X(i))/h(i),5);
    end
    %%畫樣條函數的圖像
    for i = 1:n
        t_pic{i} = X(i):(X(i+1)-X(i))/precision:X(i+1);
        T_pic{i} = subs(CSI{2,i},t_pic{i});
    end
    for i = 1:n
        t_Pic(1+(precision+1)*(i-1):(precision+1)*i) = t_pic{i};
        T_Pic(1+(precision+1)*(i-1):(precision+1)*i) = T_pic{i};
    end
    plot(t_Pic,T_Pic,'g');
    
    if attribute ==0
        legend('F:數據點','Y_real:真實圖像','T_Pic:樣條插值擬合函數');
    elseif attribute == 1
        legend('F:數據點','T_Pic:樣條插值擬合函數');
    end
elseif bc == 0
    left_endpoint = input('輸入左端點的二階導數值：');
    right_endpoint = input('輸入右端點的二階導數值：');
    for i = 1:n
        X_one_interval{i} = [X(i),X(i+1)];
        F_one_interval{i} = [F(i),F(i+1)];
        F_all_bad(i) = The_first_order_All_bad(F_one_interval{i},X_one_interval{i},attribute);
        h(i) = X(i+1)-X(i);
    end
    for i = 1:n-1
        mu(i) = h(i)/(h(i)+h(i+1));
    end
    mu(n) = 0;
    lambda(1) = 0;
    for i = 2:n
        lambda(i) = h(i)/(h(i-1)+h(i));
    end
    d(1)=2*left_endpoint;
    for i = 2:n
        d(i) = 6*(F_all_bad(i)-F_all_bad(i-1))/(h(i-1)+h(i));
    end
    d(n+1) = 2*right_endpoint;
    A = zeros(n+1,n+1);
    for i = 1:n+1
        for j = 1:n+1
            if i == j
                A(i,j) = 2;
            elseif (i-j) == 1
                A(i,j) = mu(j);
            elseif (j-i) == 1
                A(i,j) = lambda(i);
            end
        end
    end
    disp('樣條函數初始系數：');
    M = vpa((A\d')',6)
    syms x;
    for i = 1:n
        CSI{1,i} = X_one_interval{i};
        CSI{2,i} = vpa((X(i+1)-x)^3*M(i)/(6*h(i))+(x-X(i))^3*M(i+1)/(6*h(i))+(F(i)-M(i)*h(i)^2/6)*(X(i+1)-x)/h(i)+(F(i+1)-M(i+1)*h(i)^2/6)*(x-X(i))/h(i),5);
    end
    for i = 1:n
        t_pic{i} = X(i):(X(i+1)-X(i))/precision:X(i+1);
        T_pic{i} = subs(CSI{2,i},t_pic{i});
    end
    for i = 1:n
        t_Pic(1+(precision+1)*(i-1):(precision+1)*i) = t_pic{i};
        T_Pic(1+(precision+1)*(i-1):(precision+1)*i) = T_pic{i};
    end
    plot(t_Pic,T_Pic,'g');
    if attribute ==0
        legend('F:數據點','Y_real:真實圖像','T_Pic:樣條插值擬合函數');
    elseif attribute == 1
        legend('F:數據點','T_Pic:樣條插值擬合函數');
    end
end
end

　　2.一階均差函數

%%一階均差
function FAb = The_first_order_All_bad(f,X,attribute)
[m,n] = size(X);
if attribute == 0
    F = subs(f,X);
    FAb = (F(n)-F(1))/(X(n)-X(1));
elseif attribute == 1
    FAb = (f(n)-f(1))/(X(n)-X(1));
end
end

　　3.樣條插值擬合值函數

%%三次樣條插值擬合值
%%用法同三次樣條函數
function CSIV = Cubic_spline_interpolation_value(X,Y,precision,attribute,bc,x_value)
[m,n] = size(X);a = min(X);b = max(X);n = n-1;
X = sort(X);
if attribute == 0
    F = subs(Y,X);
elseif attribute ==1
    F = Y;
end
if bc == -1
    left_endpoint = input('輸入左端點的一階導數值：');
    right_endpoint = input('輸入右端點的一階導數值：');
    for i = 1:n
        X_one_interval{i} = [X(i),X(i+1)];
        F_one_interval{i} = [F(i),F(i+1)];
        F_all_bad(i) = The_first_order_All_bad(F_one_interval{i},X_one_interval{i},attribute);
        h(i) = X(i+1)-X(i);
    end
    for i = 1:n-1
        mu(i) = h(i)/(h(i)+h(i+1));
    end
    mu(n) = 1;
    lambda(1) = 1;
    for i = 2:n
        lambda(i) = h(i)/(h(i-1)+h(i));
    end
    d(1)=6*(F_all_bad(1)-left_endpoint)/h(1);
    for i = 2:n
        d(i) = 6*(F_all_bad(i)-F_all_bad(i-1))/(h(i-1)+h(i));
    end
    d(n+1) = 6*(right_endpoint-F_all_bad(n))/h(n);
    A = zeros(n+1,n+1);
    for i = 1:n+1
        for j = 1:n+1
            if i == j
                A(i,j) = 2;
            elseif (i-j) == 1
                A(i,j) = mu(j);
            elseif (j-i) == 1
                A(i,j) = lambda(i);
            end
        end
    end
    M = (A\d')';
    syms x;
    for i = 1:n
        S{1,i} = X_one_interval{i};
        S{2,i} = vpa((X(i+1)-x)^3*M(i)/(6*h(i))+(x-X(i))^3*M(i+1)/(6*h(i))+(F(i)-M(i)*h(i)^2/6)*(X(i+1)-x)/h(i)+(F(i+1)-M(i+1)*h(i)^2/6)*(x-X(i))/h(i),5);
    end
    for i =1:n
        if x_value >= X(i) && x_value <= X(i+1);
            s = i;
        end
    end
    CSIV{1,1} = '擬合值';
    CSIV{2,1} = vpa(subs(S{2,s},x_value),4);
    CSIV{1,2} = '實際值';
    CSIV{2,2} = vpa(subs(Y,x_value),4);
    CSIV{1,3} = '誤差';
    CSIV{2,3} = abs(CSIV{2,1}-CSIV{2,2});
    
elseif bc == 0
    left_endpoint = input('輸入左端點的二階導數值：');
    right_endpoint = input('輸入右端點的二階導數值：');
    for i = 1:n
        X_one_interval{i} = [X(i),X(i+1)];
        F_one_interval{i} = [F(i),F(i+1)];
        F_all_bad(i) = The_first_order_All_bad(F_one_interval{i},X_one_interval{i},attribute);
        h(i) = X(i+1)-X(i);
    end
    for i = 1:n-1
        mu(i) = h(i)/(h(i)+h(i+1));
    end
    mu(n) = 0;
    lambda(1) = 0;
    for i = 2:n
        lambda(i) = h(i)/(h(i-1)+h(i));
    end
    d(1)=2*left_endpoint;
    for i = 2:n
        d(i) = 6*(F_all_bad(i)-F_all_bad(i-1))/(h(i-1)+h(i));
    end
    d(n+1) = 2*right_endpoint;
    A = zeros(n+1,n+1);
    for i = 1:n+1
        for j = 1:n+1
            if i == j
                A(i,j) = 2;
            elseif (i-j) == 1
                A(i,j) = mu(j);
            elseif (j-i) == 1
                A(i,j) = lambda(i);
            end
        end
    end
    M = (A\d')';
    syms x;
    for i = 1:n
        S{1,i} = X_one_interval{i};
        S{2,i} = vpa((X(i+1)-x)^3*M(i)/(6*h(i))+(x-X(i))^3*M(i+1)/(6*h(i))+(F(i)-M(i)*h(i)^2/6)*(X(i+1)-x)/h(i)+(F(i+1)-M(i+1)*h(i)^2/6)*(x-X(i))/h(i),5);
    end
    CSIV{1,1} = '擬合值';
    CSIV{2,1} = vpa(subs(S{2,s},x_value),4);
end

　　4.例子

clear all
clc
syms x;
Y=sin(x)/(1+x^2)+exp(-x^2);
X=-5:1:5;
precision=500;
attribute=0;
bc=-1;

%%三次樣條插值
Cubic_spline_interpolation(X,Y,precision,attribute,bc)

　　結果為

輸入左端點的一階導數值：subs(diff(Y),-5)
輸入右端點的一階導數值：subs(diff(Y),5)
樣條函數初始系數：
M =
[ -0.0869325, 0.0232929, -0.00623918, 0.00166378, -0.000415945, 1.32899e-12, 0.000415945, -0.00166378, 0.00623918, -0.0232929, 0.0869325]
S =
  2×10 cell 數組
  列 1 至 4
    {1×2 double}    {1×2 double}    {1×2 double}    {1×2 double}
    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }
  列 5 至 8
    {1×2 double}    {1×2 double}    {1×2 double}    {1×2 double}
    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }    {1×1 sym   }
  列 9 至 10
    {1×2 double}    {1×2 double}
    {1×1 sym   }    {1×1 sym   }
>> S{2,:}
ans =
0.014489*(x + 4.0)^3 - 0.010735*x + 0.0038822*(x + 5.0)^3 - 0.0023031
ans =
- 0.053584*x - 0.0010399*(x + 4.0)^3 - 0.0038822*(x + 3.0)^3 - 0.1737
ans =
0.0010399*(x + 2.0)^3 - 0.15087*x + 0.0002773*(x + 3.0)^3 - 0.46557
ans =
0.11103*x - 0.0002773*(x + 1.0)^3 - 0.000069324*(x + 2.0)^3 + 0.058248
ans =
1.0528*x + 2.215e-13*(x + 1.0)^3 + 0.000069324*x^3 + 1.0
ans =
0.000069324*x^3 - 2.215e-13*(x - 1.0)^3 - 0.21145*x + 1.0
ans =
1.3766 - 0.0002773*(x - 1.0)^3 - 0.000069324*(x - 2.0)^3 - 0.58809*x
ans =
0.0010399*(x - 2.0)^3 - 0.18726*x + 0.0002773*(x - 3.0)^3 + 0.57497
ans =
0.17469 - 0.0010399*(x - 4.0)^3 - 0.0038822*(x - 3.0)^3 - 0.053831*x
ans =
0.014489*(x - 4.0)^3 - 0.010735*x + 0.0038822*(x - 5.0)^3 + 0.0023042
>> 

　　函數圖像為