Call Now: (800) 537-1660  
The Algebra Buster
The Algebra Buster


August 31st









August 31st

INTRODUCTION TO MATLAB

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-5L, to approximate the second derivative use
h = 10-4L.

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 = xj in the array by using the function
values fj+1 and fj-1 like this

This works fine for an N element array at all points from x2 to xN-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 x2 through xN-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(xj)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')

Prev Next
 
Home    Why Algebra Buster?    Guarantee    Testimonials    Ordering    FAQ    About Us
What's new?    Resources    Animated demo    Algebra lessons    Bibliography of     textbooks
 

Copyright © 2009, algebra-online.com. All rights reserved.