Blog

Computing derivatives numerically by taking differences

In the lesson 2 we discussed taking numerical derivatives using the finite difference method.

We considered the function y=x^2.  The task is to numerically calculate the first derivative by taking differences.

For example setting up an array for x :

x= [-3   -2  1  0  1  2  3]

with the step size

dx = (x(end)-x(1) )/(size(x,2)-1) = (3-(-3))/(7-1) = 6/6 = 1

we get the following array for the function y=x.^2

y = [9  4  1  0 1  4  9]

To compute the numerical central derivatives we take the difference between the  arrays that are shifted by two positions and divide by twice the step size. The factor 2 for the step size is because we shifted by two positions.

dydx =
( y(3:end)-y(1:end-2)  ) / (2*dx) =
(       [   1      0     1    4   9]
–    [  9      4     1     0   1]     ) / 2
=        [  -8  -4     0    4  8 ]    / 2
=        [   -4   -2    0   2  4]

which are the values of the first derivative of y for the corresponding values of x.

The array of dydx we got so far is two elements shorter than the array of y, it misses the first and the last element. To make both array the same size we add one additional element to the front and one to the back of dydx by simply extrapolation.

For example, we see that the difference between the last two elements of dydx is 4-2=2. Thus we add 2 to the last element to get to 6=4+(4-2).  Similarly we extend the beginning of dydx with -4+(-4 – (-2) ) =  -6. Thus we get

dydx = [ -6 -4 -2  0  2  4  6]

This corresponds to the expected result of dydx = 2x .

—————————————–

The task for the lesson was to compute the second derivative d2ydx for the function y=x.^3 , using the three point formula

d2ydx(i) =  ( y(i-1) – 2*y(i) + y(i+1) ) / (dx^2)

The formula is derived from first taking the one sided first differences

dydx(i+1) =  (  y(i+1)- y(i    ) )/dx
dydx(i) =      (  y(i)     – y(i-1) )/dx

and taking the difference of the differences above

d2ydx(i) =  (  dydx(i+1)-dydx(i-1) )/dx
(  (y(i+1)-y(i))-(y(i)-y(i-1))  ) / (dx^2)

Below is the starter code for the task, just missing the specification for d2ydx:

x = linspace(-3,3,61); 
dx = (x(end)-x(1))/(size(x,2)-1);
% Function y(x) definition
y = x.^3;

% First derivative
dydx = y(3:end)-y(1:end-2);
dydx = [(2*dydx(1)-dydx(2)) dydx (2*dydx(end)-dydx(end-1))]/(2*dx);

% Second derivative
d2ydx = FILL IN THE THE MISSING CALCULATION HERE
d2ydx = [(2*d2ydx(1)-d2ydx(2))   d2ydx   (2*d2ydx(end)-d2ydx(end-1))]/(dx^2);

% Plotting
plot(x,y,'b-'); hold on; plot(x,dydx,'g-'); plot(x,d2ydx,'r-'); 
xlabel('x');ylabel('y');
title('Numerical derivatives');legend('y=x^3','first deriv.','second deriv.'); 
grid on; grid minor; hold off;

Numerical derivatives for y=x^3

Initial Post — Week One

Thank you for spending your Saturday morning with your fellow students and your instructor in the applied analysis lab. I uploaded last weeks lecture slides as PDF to the OpenLab course page, including some corrections. In particular the xline function to create a vertical line did not seem to work for the older Matlab version on student computers. I provided a work around in the slides.

To recap, last week we practiced:

  • Using Matlab
    • Use Live Scrip notebooks for creating documents that contain text, code and graphs
    • Create Matlab single variables and arrays
    • Perform simple computations
    • Create a Matlab function
    • Create a graph in Matlab
  • Using complex numbers in Matlab, including:
    • Creating a complex number.
    • Cartesian and polar coordinate (Euler theorem) representation of a complex number.
    • Calculating absolute value, phase and complex conjugate of a complex number
  • RCL circuits
    • Application of Kirchhoff’s law
    • Impedance of resistor, capacitor and inductor
    • Calculate total impedance of simple circuits and resulting AC current phase shifts
    • Studying a RCL circuit simulation to find the resonance frequency

Looking forward to practice more of the above and also to continue where we left off. Please don’t be shy to raise questions if anything is unclear.

See you in class!