Homec4science
Diffusion A study of Hermite spline interpolation (master) Browse

R3913/Code/derivativeComparison.m
041b8992f8b7master

Last Change
Raw File

derivativeComparison.m
View Options

close all;
clear;
clc;
a = 0; % Lower bound of the interval. Keep it to 0
b = 1; % Upper bound of the interval. Keep it to 1
rangeFunctions = [1]; % Indices of functions we we want to load (see getFunctions.m)
Npts = 10 * 2.^(1:1:11) + 1; % Number of sampling points in [a, b]
% Npts = 13;
assert(sum(mod(Npts, 4)) == length(Npts));
L_b = 128; % Number of points on each side of the interval [a, b]
drawFun = 1; % Algo plots fun/deriv approx. if it is set to 1
drawError = 1; % Algo plots error graphs at the end
precision = 18; % Precision of numbers. Useful when 0 is -1e-18
M = 1001;
numFunctions = length(rangeFunctions);
[Functions, dFs, funNames] = getFunctions(rangeFunctions);
for k=1:numFunctions
ErrorF = zeros(2, length(Npts));
ErrorDF = zeros(4, length(Npts));
disp(['k = ', num2str(k)]);
for l=1:length(Npts)
disp(['Npts = ', num2str(Npts(l))]);
I = linspace(a, b, M);
f = @(x) Functions{k}(x);
df = @(x) dFs{k}(x);
h = (b - a) / (Npts(l)-1);
extendedX = linspace(a - L_b * h, b + L_b * h, Npts(l) + 2 * L_b);
X = round(shrink(extendedX, L_b), precision);
Samples = f(extendedX);
forward = forward_differences(extendedX, Samples, L_b);
forward_cont = lerp(forward, I, Npts(l));
centered = centered_differences(extendedX, Samples, L_b);
centered_cont = lerp(centered, I, Npts(l));
% Cubic B-Spline
c_k = computeCubicSplineCoefficients(Samples, L_b);
my_cubic = computeCubicInterpolation(c_k, X, h, Npts(l), L_b);
my_cubic_cont = computeCubicInterpolation(c_k, I, h, Npts(l), L_b);
my_cubic_der = computeCubicDerivativeInterpolation(c_k, X, h, Npts(l), L_b);
my_cubic_cont_der = computeCubicDerivativeInterpolation(c_k, I, h, Npts(l), L_b);
% Hermite
H = 2 * h;
hermNpts = ceil(0.5 * Npts(l));
hermX = X(1:2:end);
herm_c = f(hermX);
herm_d = df(hermX) * H;
my_herm = computeHermiteInterpolation(herm_c, herm_d, X, H, hermNpts);
my_herm_cont = computeHermiteInterpolation(herm_c, herm_d, I, H, hermNpts);
my_herm_der = computeHermiteDerivativeInterpolation(herm_c, herm_d, X, H, hermNpts);
my_herm_cont_der = computeHermiteDerivativeInterpolation(herm_c, herm_d, I, H, hermNpts);
exact_f_cont = f(I);
exact_derivative_cont = df(I);
nf2 = norm(exact_f_cont)^2;
ndf2 = norm(exact_derivative_cont)^2;
ErrorF(1, l) = sqrt(norm(exact_f_cont - my_cubic_cont)^2 / nf2);
ErrorF(2, l) = sqrt(norm(exact_f_cont - my_herm_cont)^2 / nf2);
ErrorDF(1, l) = sqrt(norm(exact_derivative_cont - forward_cont)^2 / ndf2);
ErrorDF(2, l) = sqrt(norm(exact_derivative_cont - centered_cont)^2 / ndf2);
ErrorDF(3, l) = sqrt(norm(exact_derivative_cont - my_cubic_cont_der)^2 / ndf2);
ErrorDF(4, l) = sqrt(norm(exact_derivative_cont - my_herm_cont_der)^2 / ndf2);
if(drawFun == 1)
figure;
plot(I, f(I), 'k');
hold on;
plot(X, my_cubic, '--ob', X, my_herm, '--or');
% plot(I, my_cubic_cont, '--ob', I, my_herm_cont, '--or');
legend('f', 'B-spline', 'Hermite', 'location', 'southwest');
xlabel('x'); grid;
ylabel('function approximation');
title(['Approximation of f, Npts = ', num2str(Npts(l)), ...
', L_b = ', num2str(L_b)]);
figure;
plot(I, df(I), 'k');
hold on;
plot(X, my_cubic_der, '--ob', X, my_herm_der, '--or');
% plot(I, my_cubic_cont_der, '--ob', I, my_herm_cont_der, '--or');
legend('f''', 'B-spline', 'Hermite', 'location', 'southwest');
xlabel('x'); grid;
ylabel('derivative approximation');
title(['Approximation of f'', Npts = ', num2str(Npts(l)), ...
', L_b = ', num2str(L_b)]);
end
end
if(drawError == 1)
figure;
loglog(Npts, ErrorF(1, :), '--xb', Npts, ErrorF(2, :), '--xr',...
Npts, 0.5 * Npts.^(-2), '--k', Npts, 1 * Npts.^(-3), '-.k');
legend('cubic Spline', 'cubic hermite', 'O(Npts^{-2})', 'O(Npts^{-3})');
title([funNames{k}, ', L_b = ', num2str(L_b)]);
xlabel('Npts');
ylabel('Error on f');
grid;
figure;
loglog(Npts, ErrorDF(1, :), '--xg', Npts, ErrorDF(2, :), '--xc',...
Npts, ErrorDF(3, :), '--xb', Npts, ErrorDF(4, :), '--xr',...
Npts, 0.6e1 * Npts.^(-1), '--k', Npts, 0.5 * Npts.^(-2), '-.k');
legend('forward diff', 'centered diff', 'cubic Spline', ...
'cubic hermite', 'O(Npts^{-2})', 'O(Npts^{-3})');
title([funNames{k}, ', L_b = ', num2str(L_b)]);
xlabel('Npts');
ylabel('Error on derivative');
grid;
end
end
c4science · Help