% Class 23 - Review of Numerical Differentiation
% Coded by NFR
% 4.9.2019
%
% Warm up activity
D = xlsread('Gluc_v2');
mdl = fitlm(D(:,1:2),D(:,3),'quadratic')
% From inspection, we should drop the cross term x1:x1
New_mdl = removeTerms(mdl,'x1:x2')
% Plot the reduced model
P = New_mdl.Coefficients.Estimate
func = @(x,y) P(1) + P(2)*x + P(3)*y...
    + P(4)*x.^2 + P(5)*y.^2;
xvec = linspace(5,7,100);
yvec = linspace(1,2.5,100);
[xmat, ymat] = meshgrid(xvec,yvec);
zmat = func(xmat,ymat);
contourf(xvec,yvec,zmat,'ShowText','on')
xlabel('Growth Time')
ylabel('Exposure Time')
% Third method of improving numerical derivatives
% [other methods are
% 1) reducing the step size (h)
% 2) including more points ]
% In romberg method, we use two estimates 
% and combine to make a more accurate estimate
%
% Example is 21.4 from textbook
F1 = @(x) cos(x);
h1 = pi/3;
h2 = h1/2;
x = pi/4;
D1 = (F1(x+h1)-F1(x-h1))/(2*h1)
D2 = (F1(x+h2)-F1(x-h2))/(2*h2)
Dimp = 4/3*D2-1/3*D1
% Example of Partial Derivative
% 21.38 from textbook
F2 = @(x,y) 3*x.*y+3*x-x.^3-3*y.^3;
% Evaluate dF2/dy analytically
dF2 = @(x,y) 3*x-9*y^2;
% Also do this numerically
h = 0.0001;
x = 1; y = 1;
dF2_num = (F2(x,y+h)-F2(x,y-h))/(2*h)
dF2(1,1)
% Two built in functions for derivatives of
% data collected in table
% example data
tvec = 1:10;
pvec = [0 2 5 7 10 15 20 40 45 60];
subplot(1,2,1)
plot(tvec,pvec,'o')
xlabel('Time (s)')
ylabel('Position (m)')
% diff function will subtract each point from the point after
% to use diff to find the velocity
vel_vec = diff(pvec)./diff(tvec);
% These would be velocities at the midpoints
% to plot, create a vector of midpoints
n = length(tvec);
tmvec = (tvec(1:n-1)+tvec(2:n))./2;
subplot(1,2,2)
plot(tvec,vel_2)
plot(tmvec,vel_vec,'o')
xlabel('Time (s)')
%
% You can also use the gradient function
vel_2 = gradient(pvec,1)
figure
plot(tvec,vel_2,'o')
ylabel('Velocity (m/s)')
% Use the quiver function to visualize gradients
% also another example of partial derivative
%
f = @(x,y) y-x-2*x.^2-2*x.*y-y.^2;
[xmat,ymat] = meshgrid(-2:.25:0, 1:.25:3);
zmat = f(xmat,ymat);
[fx,fy] = gradient(zmat,0.25);
contour(xmat,ymat,zmat);
hold on
quiver(xmat,ymat,-fx,-fy)
hold off




%
%{
% Practice
f = @(x,y) y-x-2*x.^2-2.*x.*y-y.^2
[x,y] = meshgrid(-2:.25:0, 1:.25:3);
z = f(x,y);
[fx,fy] = gradient(z,0.25)
[C,h] = contour(x,y,z);
clabel(C,h);
hold on
quiver(x,y,-fx,-fy)
hold off
%}