Numerical Differentiation
Derivatives are everywhere in physics. Velocity is the time derivative of position. Acceleration is the time derivative of velocity. Derivatives tell us about the slope of a function and the minimum or maximum of a function occurs where the derivative is zero. We use these facts all the time so let’s investigate some numerical derivatives and their applications.
Finite difference method
Consider a function like, y = x5. Matlab represents this function as an array of y values evaluated at x values.
clear all; close all;
x = 0:0.1:10;
y = x.^5;
The numerical derivative is calculated using the pre-calculus method of finding Dy and dividing by Dx. This is only a good approximation to the derivative if D is small, but it’s the best we’ve got. We are just going to calculate the slope, Dy/Dx at every point.
dy = y(2:end) - y(1:end-1); % yn+1 – yn
dx = x(2:end) - x(1:end-1);
dydx = dy./dx;
We have used a notation called "slicing" to pick just some elements of our array.
This is done by explicitly telling Matlab which elements of the array to use. The last index in an array is named "end". To get a sense of how this works try these and be sure you understand what is happening.
z=0:0.5:4
z(3:end)
z(2:end-1)
The following won’t work. Can you see why?
plot(x,dydx); % this doesn’t work
The vector x is longer by one element than the derivative vector dydx. We have to define the x values for dydx. There are two ways to pick these:
xSingle = x(1:end-1); % single-sided
xMidPt = (x(2:end) + x(1:end-1) )/2; % midpoint
plot(xSingle, dydx, xMidPt, dydx)
Notice that the xSingle is a bit to the left of xMidPt. Although these are nice ways to find the derivative, keeping track of the arrays will become a big pain. I like to use Matlab’s interpolation routine to get back to the original spacing.
dydx = interp1(xMidPt,dydx,x,'cubic');
Remember how to use interpolation?
Accuracy of Finite difference method
Let’s plot the error introduced by our numerical derivatives. The exact derivative is 5x4.
dydxExact = 5*x.^4;
Let’s plot the percent difference between the real derivative and our double-sided derivative calculated at the same x values:
percentError = 100*abs( dydx - dydxExact )./dydxExact;
semilogy(x, percentError)
Assignment: M6_Differentiation
1. It will be easier to use a few lines of Matlab code over and over if they are in a file by themselves. Write a reusable function called finiteDifference.m that is used like …
dydx = finiteDifference(x,y);
Now write a script to that takes the derivative.
x = 0:0.1:10;
y = x.^5;
dydx = finiteDifference(x,y);
plot(x,dydx);
How do you know it works correctly?
What is the accuracy with 0.001 spacing? For what spacing is the derivative accurate to 4 significant figures?
2. LabPro has an ultrasound detector that can measure the position of a cart on a track as a function of time. The CIA has just obtained two files of such data from Baghdad and important people in Washington want you to tell them what it means. The files are called cia_file1.txt and cia_file2.txt. Copy them to your working area and read the data in using Matlab’s textread function. (Use "help textread" to find out more.)
[t,x,dummy] = textread('cia_file1.txt','','commentstyle','matlab');
Plot x vs t to see what you’ve got and then find the velocity and acceleration using your finiteDifference function. Try smoothing the derivatives (help smooth) to get a clearer picture. Is the information clearer if you smooth on velocity or acceleration or both? Include a plot of x, v, and a vs t in your report to the CIA bosses with a description of what happened in Baghdad.
(Some versions of Matlab do not include the “smooth” function. If you find that its not available in your version of Matlab, please download the “smoothJLC.m” file from the course web page and try using it.)
M4_Differentiation_mjm.doc 5/5/08 Page 2 of 2