%differentiation.m

%Numerical
Differentiation

%Engineering
6, U.C. Davis

%The
derivative f'(x) may be interpreted geometrically as the

% slope
of the function f at the point x, where the slope

% is the
slope of the tangent line at that point.

%So
numerical differentiation techniques estimate the

%
derivative of a function at a point x_k by approximating

% the
slope of the tangent line at x_k using other known

% values
of the function near x_k.

%Three
ways to approximate the slope at x_k (Fig. 11.7):

%1.Backward
difference approximation

% This
calculates the slope by drawing a line between

% f(x_k)
and f(x_k_prev), where x_k_prev is the value

% of x
just previous to x_k that has a known value for f(x).

% The
slope is then simply:

% f'(x_k) = slope = (f(x_k) - f(x_k_prev))/(x_k
- x_k_prev)

% If the
samples of x are uniformly spaced, then we can

% write
x_k - x_k_prev = delta_x.

%2.
Forward difference approximation

% This
calculates the slope by drawing a line between

% f(x_k)
and f(x_k_next), where x_k_next is the next value

% of x
after x_k that has a known value for f(x).

% The
slope is then simply:

% f'(x_k) = (f(x_k_next) - f(x_k))/(x_k_next -
x_k)

% Again,
if the samples are uniformly spaced, then

%
x_k_next - x_k = delta_x.

%3.
Central difference approximation

% As can
be seen in Figure 11.7, however, the backward

% and
forward difference approximations for f'(x)

% can
lead to significantly different answers. In general,

% the
quality of the result depends on the spacing between

% the
data points (closer is better) and the scatter

% in the
data due to measurement error.

% In
order to lessen the effects of scatter, we calculate

% the
central difference approximation, which is the

%
average of the backward and forward difference approximations.

% Doing
the math (p. 232 of the notes) yields:

% f'(x_k) = slope = (f(x_k_next) -
f(x_k_prev))/2*delta_x

%Approximating
the Second Derivative

%The
second derivative may be approximated by finding the

% slope
of the first derivative:

% Using
the backward difference approximation:

% f''(x_k) = (f'(x_k) - f'(x_k_prev))/(x_k -
x_k_prev)

% Then
substitute the backward difference approximations for

%
f'(x_k) and f'(x_k_prev) into this equation (result shown

% on p.
232 of notes).

%
Similar results may of course also be derived using either

% the
forward difference or central difference approximations.

clc; clf; format compact;

%Numerical
Differentiation with Matlab

%The
diff function

%Consider
the following data for some function y = f(x):

y = [20 30 10 50 80 100]

x = [0 5 10 15 20 25]

disp(' '); pause;

%We may
calculate the difference between adjacent values of

% x and
adjacent values of y using the diff function:

dy = diff(y)

dx = diff(x)

disp(' '); pause;

%The
approximate derivative is then:

df = diff(y)./diff(x)

disp(' '); pause;

%Consider
the third element of x (10, in this case). The

%
approximate derivative at this x using the backward

%
difference is (10-30)/(10-5) = -4. Note that this is

% the
second element of the calculated df.

%Now
consider the second element of x (5, in this case).

% The
approximate derivative at this x using the forward

%
difference is (10-30)/(10-5) = -4, which is the same as

% the
backward difference for the third value of x.

%So in
general the values of df may be used for either

% the
backward difference or the forward difference, the

% only
difference being which value of x is referred to.

%Also
note that when we consider the whole set of elements

% in the
vector x (0, 5, 10, 15, 20, 25 in this case), the

%
backward difference uses the values 5 through 25 and

% the
forward difference uses the values 0 through 20.

% The
reason is simply that it is impossible to calculate

% a
backward difference for x=0, because there is no

%
previous value of x or y. And similarly for the forward

%
difference for x=25: there is no next value to use.

%So to
plot the backward difference results for the

%
derivative at this set of x values, we would write:

x_backward = x(2:end)

disp(' '); pause;

figure(2)

subplot(2,1,1),...

plot(x_backward,df),axis([0
25 -10 10]),...

xlabel('x'), ylabel('Backward diff results')

disp(' '); pause;

%To plot
the forward difference results for the derivative

% at
this set of x values:

x_forward = x(1:end-1)

disp(' '); pause;

subplot(2,1,2),...

plot(x_forward,df),
axis([0 25 -10 10]),...

xlabel('x'), ylabel('Forward diff results')

disp(' '); pause;

%Comparing
backward difference and central difference methods

x = linspace(0,pi,50);

n = length(x);

%Create
sine signal with random Gaussian error

noisydata = sin(x) +
0.025*randn(1,n);

%The
derivative of the noiseless sine signal (i.e., derivative of

%a sine
function is a cosine)

truederiv = cos(x);

%Backward
difference estimate for the noisy sine signal

noisydata_BD =
diff(noisydata)./diff(x);

clf;

subplot(2,1,1),
plot(x(2:n),truederiv(2:n),x(2:n),noisydata_BD,'ro'),...

xlabel('x'),ylabel('Derivative'),
axis([0 pi -2 2]),...

legend('True derivative', 'Backward
difference')

%Central
difference estimate

noisydata_CD =
(noisydata(3:n) - noisydata(1:n-2))./(x(3:n) - x(1:n-2));

subplot(2,1,2),
plot(x(2:n-1),truederiv(2:n-1),x(2:n-1),noisydata_CD,'ro'),...

xlabel('x'),ylabel('Derivative'),
axis([0 pi -2 2]),...

legend('True derivative', 'Central
difference')

disp(' '); pause;

%Numerical
Differentiation of a Polynomial Function

%Given
the linear acceleration of an object whose speed is

%
defined by s(t) = t^3 - 2t^2 +2 m/s over the interval t =

% 0 to 5
seconds, determine the specific acceleration at

% t =
2.5 seconds.

%Define
t and s vectors (i.e., pretend we just have some data on the

% speed
for various times t, not the actual function s(t) above)

t = 0:0.1:5;

s = t.^3 - 2*t.^2 + 2;

%The
acceleration is given by the derivative of the speed

ds = diff(s)./diff(t);

%To find
the acceleration at t = 2.5 seconds, we need to figure out which

% index
it is in the t vector. The basic equation is

% t_given
= delta_t*(k - 1)

% where
k is the index and delta_t is the time interval (0.1 in

% this
case). Solving for k: k = (t_given/delta_t) + 1. So:

k = 2.5/0.1 + 1

disp(' '); pause;

%Find
the acceleration at that t

%(Note:
It’s not immediately obvious, but with k given as above,

% we are
doing a forward difference here.)

accel = ds(k)

disp(' '); pause;

%Do it
again for a delta_t of 0.01 seconds

t = 0:0.01:5;

s = t.^3 - 2*t.^2 + 2;

ds = diff(s)./diff(t);

k = 2.5/0.01 + 1

disp(' '); pause;

accel = ds(k)

disp(' '); pause;

%And
once again for a delta_t of 0.001 seconds

t = 0:0.001:5;

s = t.^3 - 2*t.^2 + 2;

ds = diff(s)./diff(t);

k = 2.5/0.001 + 1

disp(' '); pause;

accel = ds(k)

disp(' '); pause;

%The
analytical result is 8.75 m/s^2 (calculated by

%
differentiating the polynomial and substituting 2.5 s).

% So we
see that increasing the resolution of the time

%
interval (i.e., the closer together the data points

% are) gives
better results.

%Differentiation
Error Sensitivity

%Generate
some noisy data

x = 0:0.1:1;

y = [-0.447 1.978 3.28
6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];

%Fit a
second degree curve to this data

a = polyfit(x,y,2);

%Plot
the data and curve

x_fitted =
linspace(0,1,101);

y_fitted =
polyval(a,x_fitted);

subplot(2,1,1),plot(x,y,'--o',x_fitted,y_fitted),...

xlabel('x'), ylabel('y'),...

title('Noisy data and 2nd degree
curve fit'),...

legend('Noisy data','2nd
degree curve fit',4);

%Differentiate
using the noisy data (use forward difference)

y_FD = diff(y)./diff(x);

x_FD = x(1:end-1); %Set up for the forward difference

%Differentiate
the fitted curve

a_deriv = polyder(a); %Coefficients of the derivative polynomial

%Calculate
values of the deriv of the fitted polynomial (so can plot it)

y_fitted_deriv =
polyval(a_deriv,x_fitted);

%Plot
the derivative results from the data and from the curve

subplot(2,1,2),plot(x_FD,y_FD,'--o',x_fitted,y_fitted_deriv),...

xlabel('x'), ylabel('dy/dx'),...

title('Derivative approximations'),...

legend('Noisy data','2nd degree fit',3);