%% Finite Differences
% Consider the function f(x)=\cos(x^2). Then, f'(x)=-2x\sin(x^2).
% Sample f(t_i)at N+1 points.

clear
N=10;
t = linspace(0,1,N+1);   % nodes
f = @(x) cos(x.^2);
dfdx = @(x) -2*x.*sin(x.^2);

%
%% Forward Difference Approximation
% f'(x)\approx \frac{f(x+h)-f(x)}h

h=1/N;

forward = @(x)(f(x+h)-f(x))/h;

x0=t(N/2+1);
exact_value=dfdx(x0)
approx_value=forward(x0)
err=abs(exact_value-approx_value)

%

t2 = linspace(0,1,50);  % For plotting exact derivatives
plot(t2,dfdx(t2),t,forward(t),'ro')
xlabel('x')
ylabel('-2xsin(x^2)')
title('Forward Difference')
legend('Exact','Approx')

%
%% Backward Difference Approximation
% f'(x)\approx \frac{f(x)-f(x-h)}h

backward = @(x)(f(x)-f(x-h))/h;

x0=t(N/2+1);
exact_value=dfdx(x0)
approx_value=backward(x0)
err=abs(exact_value-approx_value)

%

plot(t2,dfdx(t2),t,backward(t),'ro')
xlabel('x')
ylabel('-2xsin(x^2)')
title('Backward Difference')
legend('Exact','Approx')

%
%% Centered Difference Approximation
% f'(x)\approx \frac{f(x+h)-f(x-h)}{2h}

centered = @(x)(f(x+h)-f(x-h))/(2*h);

x0=t(N/2+1);
exact_value=dfdx(x0)
approx_value=centered(x0)
err=abs(exact_value-approx_value)

%

plot(t2,dfdx(t2),t,centered(t),'ro')
xlabel('x')
ylabel('-2xsin(x^2)')
title('Centered Difference')
legend('Exact','Approx')

%
%% Second Derivative
% f''(x)\approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}

dfdx2 = @(x) -2*sin(x.^2)-4*x.^2.*cos(x.^2);
ddf = @(x)(f(x+h)-2*f(x)+f(x-h))/(h^2);

x0=t(N/2+1);
exact_value=dfdx2(x0)
approx_value=ddf(x0)
err=abs(exact_value-approx_value)

%

plot(t2,dfdx2(t2),t,ddf(t),'ro')
xlabel('x')
ylabel('f''''(x)')
title('Second Derivative')
legend('Exact','Approx')

%
%% Fornberg Weights
% Estimate the derivative of \cos(x^2) at x=0.5 using five nodes.
%

t = [ 0.35 0.5 0.57 0.6 0.75 ]';   % nodes
f = @(x) cos(x.^2);
dfdx = @(x) -2*x.*sin(x.^2);
exact_value = dfdx(0.5)

%
% Shift the nodes so that the point of estimation for the derivative is at x=0.
%

w = fdweights(t-0.5,1)
fd_value = w'*f(t)

%
% We can reproduce the weights in the finite difference tables by using equally spaced nodes with h=1. For example, here are two one-sided formulas.

format rat         % Produce rational output
fdweights(0:2,1)   % a[0] f(x) + a[1] f(x+h) + a[2] f(x+2h)
fdweights(-3:0,1)  % a[-3] f(x-3h) + a[-2] f(x-2h) + a[-1] f(x-h) + a[0] f(x)

fdweights(-1:1,2)  % Second derivative weights

%

function w = fdweights(t,m)
%FDWEIGHTS   Fornberg's algorithm for finite difference weights.
% Input:
%   t    nodes (vector, length r+1)
%   m    order of derivative sought at x=0 (integer scalar)
% Output:
%   w    weights for the approximation to the jth derivative (vector)

% This is a compact implementation, not an efficient one. 

r = length(t)-1;
w = zeros(size(t));
for k = 0:r
  w(k+1) = weight(t,m,r,k);
end
end

function c = weight(t,m,r,k)
% Implement a recursion for the weights.
% Input:
%   t   nodes (vector)
%   m   order of derivative sought 
%   r   number of nodes to use from t (<= length(t))
%   k   index of node whose weight is found 
% Output:
%   c   finite difference weight 

if (m<0) || (m>r)        % undefined coeffs must be zero
  c = 0;    
elseif (m==0) && (r==0)  % base case of one-point interpolation
  c = 1;   
else                     % generic recursion 
  if k<r
    c = (t(r+1)*weight(t,m,r-1,k) - ...
        m*weight(t,m-1,r-1,k))/(t(r+1)-t(k+1));
  else
    beta = prod(t(r)-t(1:r-1)) / prod(t(r+1)-t(1:r));
    c = beta*(m*weight(t,m-1,r-1,r-1) - t(r)*weight(t,m,r-1,r-1));
  end
end
end