% 周期関数の自然スプライン補間
% x = 0, 1, 2, ..., n-1 において
% f = f(1), f(2), f(3), ..., f(n) となる周期 n の自然スプライン関数
% すなわち，x = n, n+1, n+2, ... において f = f(1), f(2), f(3), ...

f = [ 3; 2; 4; 5; 4; 2 ];
sz = size(f); n = sz(1);

e1 = ones(n,1);
e0 = 4*ones(n,1);
A = spdiags([e1 e0 e1], -1:1, n, n); % 三重対角行列
A(1,n) = 1;
A(n,1) = 1;

fprintf("係数行列\n");
full(A)

b = zeros(n,1);
b(1) = 3*(f(2)-f(1));
b(n) = 3*(f(n)-f(n-1));
b(2:n-1) = 3*(f(3:n)-f(1:n-2));

% Ad = b --> U' U d = b
[U] = chol(A);
y = U'\b;
d = U\y;

set(0,'defaultAxesFontSize',16);
set(0,'defaultTextFontSize',16);

xdata = [];
gdata = [];   % 関数値
gddata = [];  % 一階微分
gdddata = []; % 二階微分
pfin = 3*n;
step = 0.1;
for p=1:pfin
    xi=p-1; xj=p;
    x = [xi:step:xj]';
    i = mod(p-1,n)+1;
    j = mod(p,n)+1;
    fi = f(i); fj = f(j);
    di = d(i); dj = d(j);
    g = @(x) fi*phi0(x-xi)+fj*phi1(x-xi)+di*psi0(x-xi)+dj*psi1(x-xi);
    gd = @(x) fi*phi0d(x-xi)+fj*phi1d(x-xi)+di*psi0d(x-xi)+dj*psi1d(x-xi);
    gdd = @(x) fi*phi0dd(x-xi)+fj*phi1dd(x-xi)+di*psi0dd(x-xi)+dj*psi1dd(x-xi);
    xdata = [xdata; x];
    gdata = [gdata; g(x)];
    gddata = [gddata; gd(x)];
    gdddata = [gdddata; gdd(x)];
end

plot(0:n-1,f,'o', xdata,gdata);
xlim([0,pfin]);
ylim([1,6]);
xlabel('x');
ylabel('f');
grid on;
saveas(gcf, 'natural-spline-cyclic.eps', 'eps');
fprintf("自然スプライン\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;

plot(0:n-1,d,'*', xdata,gddata);
xlim([0,pfin]);
ylim([-3,3]);
xlabel('x');
ylabel('f^{\prime}');
grid on;
saveas(gcf, 'natural-spline-cyclic-d.eps', 'eps');
fprintf("自然スプラインの一階微分\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;

plot(xdata,gdddata);
xlim([0,pfin]);
ylim([-12,12]);
xlabel('x');
ylabel('f^{\prime\prime}');
grid on;
saveas(gcf, 'natural-spline-cyclic-dd.eps', 'eps');
fprintf("自然スプラインの二階微分\n");
fprintf("一時停止：何かのキーを押してください\n");
pause;
