% Class 8 - Numerical Derivatives
% 9.19.2019
% Coded by NFR
%
% Problem 21.2 from the textbook
% First derivative (centered) O(h^2)
h = 0.1;
x = 2;
x_ip1 = x+h;
x_im1 = x-h;
Func1 = @(x) exp(x);
FD_Oh2 = (Func1(x_ip1)-Func1(x_im1))/(2*h)
FD_Oh4 = (-Func1(x+2*h)+8*Func1(x+h)...
    -8*Func1(x-h)+Func1(x-2*h))/(12*h)
FD_exact = exp(x)
% Real life finite difference problem
% Rate of protein expression
% Load the fluorescence data;
F_vec = readmatrix('sample Data.xlsx','Range','I3:I290')
% Convert our fluorescent data to protein 
% concentration using a calibration curve
% (also called transfer function)
h = 5; %min
CalFunc = @(x) 0.0043*x + 0.6626;
C_vec = CalFunc(F_vec); % ug/ml
t_vec = (1:length(C_vec))*h;
plot(t_vec,C_vec,'ro','MarkerSize',3,...
    'LineWidth',3)
hold on
% Because this is noisy, we have two options
%1) We can use a model (regress the data)
%    this will be done in Unit 2
% 2) Smooth the data, in this case
%   use an average window
C_vec_S = smooth(C_vec,25);
plot(t_vec,C_vec_S,'k+','MarkerSize',3,...
    'LineWidth',3)
hold off
xlabel('Time (min)')
ylabel('Protein Conc (ug/ml)')
% Now we want to solve for the first derivative
% which is the rate of protein expression
% Preallocate to store answers
np = length(C_vec);
D_vec = zeros(np,1);
for i = 1:np
    if i == 1
        % Forward FD for first point
        Slope = (C_vec_S(i+1)-C_vec_S(i))/h;
    elseif i == np
        % Backward FD for last
        Slope = (C_vec_S(np)-C_vec_S(np-1))/h;
    else
        % Centered for rest
        Slope = (C_vec_S(i+1)-C_vec_S(i-1))/(2*h);
    end
    D_vec(i,1) = Slope;
end
figure
plot(t_vec,D_vec,'ro','MarkerSize',3,...
    'LineWidth',3)
xlabel('Time (min')
ylabel('Rate of expression (ug/(mL*min))')
%
% Code for 21.6 from Text
% Numerical derivative for not equal spaced points
% Use La Grange interpolation
F2 = @(x) 2*x^4-6*x^3-12*x-8;
% Use equation given in notes
x0 = -.5;
x1 = 1;
x2 = 2;
x = 0;
Der_0 = F2(x0)*(2*x-x1-x2)/((x0-x1)*(x0-x2))...
    +F2(x1)*(2*x-x0-x2)/((x1-x0)*(x1-x2))...
    +F2(x2)*(2*x-x0-x1)/((x2-x0)*(x2-x1))
% Compare to true
D3 = @(x) 8*x^3-18*x^2-8;
Der_0_true = D3(0)
% Compare to centered finite