Chapter 10
Derivatives and Integrals
Matlab won't give you formulas for the derivatives and
integrals of functions like a symbolic
math program. But if you have closely spaced arrays filled with values of x and
f (x)
Matlab will quickly give you numerical approximations to the derivative f'(x)
and the
integral
10.1 Derivatives
In your first calculus class the fol lowing formula for the derivative was given:
To do a derivative numerically we use the following slightly different, but
numerically more
accurate, form:
It's more accurate because it's centered about the value of x where we want the
derivative
to be evaluated.
To see the importance of centering, consider Fig. 10.1. In
this figure we are trying to
find the slope of the tangent line at x = 0.4. The usual calculus-book formula
uses the data
points at x = 0.4 and x = 0.5, giving tangent line a. It should be obvious that
using the
"centered" pair of points x = 0.3 and x = 0.5 to obtain tangent line b is a much
better
approximation.
As an example of what a good job centering does, try differentiating sin x this way:
dfdx=(sin(1+1e-5)-sin(1-1e-5))/2e-5
% take the ratio between the numerical derivative and
the
% exact answer cos(1) to see how well this does
format long e
dfdx/cos(1)
Figure 10.1 The centered derivative approximation works best.
You can also take the second derivative numerically using
the formula
For example,
d2fdx2=(sin(1+1e-4)-2*sin(1)+sin(1-1e-4))/1e-8
% take the ratio between the numerical derivative and
the
% exact answer -sin(1) to see how well this does
format long e
d2fdx2/(-sin(1))
You may be wondering how to choose the step size h. This
is a little complicated, take a
course on numerical analysis and you can see how it's done. But until you do,
here's a rough
rule of thumb. If f(x) changes significantly over an interval in x of about L,
approximate
the first derivative of f(x) using h = 10^{-5}L, to approximate the second
derivative use
h = 10^{-4}L.
If you want to differentiate a function defined by arrays x and f, then the step
size is
already de termined , you just have to live with the accuracy obtained by using h
= ∆x,
where ∆x is the spacing between points in the x array. Notice that the data must
be evenly
spaced for the example we are going to give you to work.
The idea is to approximate the derivative at x = x_{j} in
the array by using the function
values f_{j+1} and f_{j-1} like this
This works fine for an N element array at all points from x_{2} to x_{N-1}, but it
doesn't work
at the endpoints because you can't reach beyond the ends of the array to find
the needed
values of f. So we use this formula for x_{2} through x_{N-1}, then use linear
extrapolation to
find the derivatives at the endpoints, like this
Example 10.1a (ch10ex1a.m)
% Example 10.1a (Physics 330)
clear, close all,
dx=1/1000,
x=0:dx:4,
N=length(x),
f=sin(x),
% Do the derivative at the interior points all
at once using
% the colon command
dfdx(2:N-1)=(f(3:N)-f(1:N-2))/(2*dx),
% linearly extrapolate to the end points (see
the next section)
dfdx(1)=2*dfdx(2)-dfdx(3),
dfdx(N)=2*dfdx(N-1)-dfdx(N-2),
% now plot both the approximate derivative and
the exact
% derivative cos(x) to see how well we did
plot(x,dfdx,'r-',x,cos(x),'b-')
% also plot the difference between the
approximate and exact
figure
plot(x,dfdx-cos(x),'b-')
title('Difference between approximate and exact derivatives') |
The second derivative would be done the same way. Matlab
has its own routines for
doing derivatives, look in online help for diff and gradient.
10.2 Definite Integrals
There are many ways to do definite integrals numerically, and the more
accurate these
methods are the more complicated they become. But for everyday use the midpoint
method
usually works just fine, and it's very easy to code. The idea of the midpoint
method is to
approximate the integral by subdividing the
interval [a, b] into N subintervals of
width h = (b-a)/N and then evaluating f(x) at the center of each subinterval. We
replace
f(x)dx by f(x_{j})h and sum over all the subintervals to obtain an approximate
integral. This
Figure 10.2 The midpoint rule works OK if the function is nearly
straight across each interval.
method is shown in Fig. 10.2. Notice that this method
should work pretty well over subintervals
like [1.0,1.5] where f(x) is nearly straight, but probably is lousy over
subintervals
like [0.5,1.0] where the function curves .
Example 10.2a (ch10ex2a.m)
% Example 10.2a (Physics 330)
close all,
N=1000,
a=0,
b=5,
dx=(b-a)/N,
x=.5*dx:dx:b-.5*dx, % build an x array of centered values
f=cos(x), % load the function
% do the approximate integral
s=sum(f)*dx
% compare with the exact answer, which is sin(5)
err=s-sin(5) |
If you need to do a definite integral in Matlab, this is
an easy way to do it. And to see
how accurate your answer is, do it with 1000 points, then 2000, then 4000, etc.,
and watch
which decimal points are changing as you go to higher accuracy. (A function that
does this
for you is given in Chapter 12.)
Figure 10.3 Fitting parabolas (Simpson's Rule) works better.
And if you need to find the indefinite integral, in
Chapter 12 there is a piece of code
that will take arrays of (x, f(x)) and calculate the indefinite integral
function
10.3 Matlab Integrators
Matlab also has its own routines for doing definite and indefinite integrals
using both data
arrays and M- les, look in online help for the commands trapz (works on arrays,
not
functions), cumtrapz, quad, quadl, and dblquad, or in Mastering Matlab 6,
Chapter 23.
As an example, the Matlab routine quad approximates the function with parabolas
(as
shown in Fig. 10.3) instead of the rectangles of the midpoint method. This
parabolic
method is called Simpson's Rule. As you can see from Fig. 10.3, parabolas do a
much
better job, making quad a standard Matlab choice for integration. As an example
of what
these routines can do, here is the way you would use quad to integrate cos xe^{-x}
from 0 to
2:
Note that these integration routines need function M- les,
as f zero did . These will be
discussed more fully in Chapter 12, so for now just use fint.m, given below, as
a template
for doing integrals with quad, quadl, etc.
Example 10.3a (fint.m)
% Example 10.3a (Physics 330)
% define the function to be integrated in fint.m
function f=fint(x)
%*********************************************************
% Warning: Matlab will do this integral with arrays of x,
% so use .*, ./, .^, etc. in this function. If you forget
% to use the .-form you will encounter the error:
%
% Inner matrix dimensions must agree.
%
%*********************************************************
f=cos(x).*exp(-x), |
E
xample 10.3b (ch10ex3b.m)
% Example 10.3a (Physics 330)
% once fint.m is stored in your current
directory
% you can use the following commands to integrate.
% simple integral, medium accuracy
quad(@fint,0,2)
% integrate with specified relative accuracy
quad(@fint,0,2,1e-8)
% integrate with specified relative accuracy
% using quadl (notice that quadl is faster--
% always try it first)
quadl(@fint,0,2,1e-8) |
Or you can use an inline function like this to avoid
making another le:
Example 10.3c (ch10ex3c.m)
% Example 10.3b (Physics 330)
f=inline('exp(-x).*cos(x)','x')
quadl(f,0,2,1e-8) |
And if parabolas are good, why not use cubics , quartics,
etc. to do an even better job?
Well, you can, and Matlab does. The quadl integration command used in the
example
above uses higher order polynomials to do the integration and is the best Matlab
integrator
to use.
Matlab also has a command dblquad that does double
integrals. Here's how to use it.
Example 10.3d (f2int.m)
% Example 10.3d (Physics 330)
% First define the integrand as a function of x and y
function f=f2int(x,y)
f=cos(x*y), |
Example 10.3e (ch10ex3e.m)
% Example 10.3e (Physics 330)
%*********************************************************
% This is how to obtain the double integral over
% the xy rectangle (0,2)X(0,2). It runs lots
% faster if you use the 'quadl' option, as shown below
%*********************************************************
dblquad(@f2int,0,2,0,2,1e-10,'quadl') |