%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);