Integration in 1D
We often calculate integrals in physics (electromagnetism, thermodynamics, quantum mechanics, etc.). Without computers we evaluate integrals analytically. There are some programs that follow the same steps that you follow to derive the analytical solution. The Wolfram web page has a 1-D integrator: http://integrals.wolfram.com/index.jsp. On their web page you can click on compute for a variety of examples, or enter your own function. The result is a function of the input parameters and variables. Maple and Mathematica both have analytical integrators and are otherwise similar to Matlab, allowing you to evaluate the resulting functions and plot the results.
Matlab’s trapz, quad, and quadl
Some functions are too difficult to integrate analytically and for these we need to use the computer to integrate numerically. A numerical integral goes back to the basic principles of calculus. The integral is the area under a function and is given by in the limit that Matlab computes integrals numerically using the trapz function. The idea is to evaluateat many x values and then sum over the areas made from small trapezoids:
% Make the function
clear all; close all;
dx = 0.1; % dx step size
x = 0:dx:20; % array of x values [0:20]
y = x.^2; % y is a function of x
% Integrate
format long % print more digits
I1 = trapz(x,y) % trapz
I2 = dx*trapz(y) % trapz again.
The integrals, I1 and I2 should have printed in the Matlab window. To obtain precise results it is important that be a small number. Making dx too small, however, will make enormous arrays that can waste Matlab’s memory. This may also cause the calculation to take a very long time to finish. Try the above example with dx=1.0 and dx=0.01. When you calculate an integral, try the calculation with several values of dx to see that your result hasn’t changed in the significant digits you want to report.
There are several other more sophisticated integration functions in Matlab: quad, and quadl. For functions with sharp curvatures, these routines will adapt the size of the interval to produce better estimates of the integral. Users can also specify the accuracy of the integral by passing it a tolerance. The added complication is that the user has to provide a reusable function as the first argument in the function calls.
xSquared = @(x) x.^2;
tol = 1e-8;
I3 = quad(xSquared, 0, 20, tol)
I4 = quadl(xSquared, 0, 20, tol)
Here I’ve created a “handle” to an “anonymous function” by writing xSquared = @(x) x.^2. The @ (x) means that this function takes the argument x and the equation that follows is the code that will be run each time @(x) is called. We could also define the function in the call to quad.
I3 = quad(@(x) x.^2, 0, 20, tol)
In theory we can call any of our reusable functions. In the last assignment, we wrote a 1D function to turn thermistor data into a continuous function. Even though there is no reason to integrate this function, let’s do it.
Itherms = quad(@thermistor, 1, 3, tol)
The problem I usually have, however, is that my functions have additional parameters. Remember the myGauss function? It needs sigma as well as x. To integrate myGauss(x,sigma) I would do the following:
Igauss= quad(@(x) myGauss(x,0.2), -2, 2)
This will call the function myGauss.m (which I created previously) with sigma=0.2. Technically this is called wrapping my function in an anonymous function.
Trapezoid Rule and Simpson’s Rule
Although “trapz” and “quad” work nicely, take a look at how easy they are to write for yourself.
Trapezoid Rule:
If we take a bunch of trapezoids, the area of each is: .
The integral is then just the sum, .
If performed directly, this means that each needs to be evaluated twice, once for the trapezoid to its left and again for the trapezoid to its right. We can eliminate this double computation by rearranging the terms in the sum:
.
Here’s how to implement the trapezoid sum in Matlab. We’ll first make an array of coefficients that multiply f(x): [½ 1 1 1 1 1 … ½]. Then we sum.
coef = ones(size(x));
coef(1) = 1/2;
coef(end) = 1/2;
prob = dx*sum(coef.*xSquared(x))
Simpson’s Rule:
Simpson’s Rule is similar to the Trapezoid Rule, except it considers the curvature of the function. Instead of using trapezoids to estimate the area under the curve, it uses 3 adjacent points to determine a quadratic curve and integrates the area under this analytically. The result is a sum:
Except for the first and last elements, all the even elements have a coefficient of 4/3 and all the odd elements have a coefficient of 2/3.
coef = 1/3*ones(size(x));
coef(2:2:end-1) = 4/3; % Even coefficients
coef(3:2:end-1) = 2/3; % Odd coefficients
prob = dx*sum(coef.*xSquared(x))
Electric field of a rod with constant charge density:
Here is the calculation of the electric field as done previously with a spreadsheet. Remember, it was a charged rod (charge density) on the x-axis between a and b and we want the electric field at a point in the plane, (. First we'll do the constant charge density case. This time, we can use “pointField.m”, to calculate the field of all the chunks of charge along the x-axis. Instead of summing directly, however, we’ll use trapz to perform the integration.
% define constants
lambda = 1.0e-6; % charge density is 1.0 microC/m
x0 = 3.0; % position of point (x0, y0) in meters
y0 = 3.0;
dx = 0.1;
xi = 1:dx:2; % limits of the rod. a=1 to b=2
% electric field at P from a chunk of charge at xi
[E_xi,E_yi] = pointField(x0-xi, y0, lambda*dx);
% integrate along charges on x-axis.
Ex = trapz(E_xi)
Ey = trapz(E_yi)
In the spreadsheet we got: Ex = 349.49 N/C and Ey = 714.43 N/C. Here we get Ex=349.26 N/C and Ey=714.42 N/C. Which do you think is more accurate?
Assignment: M7_Integration
1. It is tempting to blindly do a calculation and trust that computers never make mistakes. Explain why the following example gives imprecise results. Then write code that runs in less than 5 seconds and integrates correctly using trapz to produce a number that is less than 10-10.
% This is very imprecise:
theta = [0:0.3:pi]
Iwrong = trapz(theta,cos(theta))
% the answer should be exactly 0.0
2. In Modern Physics or Chemistry you needed to do integrals to find the probability of finding a particle in a certain region. Consider a particle in a one-dimensional box of width L. is the probability of finding the particle between a and b where is the wavefunction. (You may need to look in a modern physics or chemistry text to get the wave functions for a particle in a box. Try Google if you don’t have a text book handy or if you are under 30.) What is the probability of finding the particle between L/3 and L/2 for the ground state and for the first excited state? To see if these make sense make a plot of vs. x for L=1. Perform the integral using both trapz and quad with an accuracy of 6 sig figs. How many bins did you need to use for trapz?
3. . Redo the variable charge density case done previously on the spreadsheet, with and a =1 m and b = 3 m, look at the point . Be sure you get results similar to what you got before. Devise a test to determine the accuracy of the Matlab calculation. Which is more accurate, Matlab’s calculation or Excel’s? Can you suggest possible reasons for the difference(s)?
4. The CIA was pleased with your analysis of the LabPro data from Baghdad. They now want to know if you can tell them about any impulses that occurred during the Baghdad experiments with information recently obtained from a secret agent that the cart in the tests had a mass of 1000 kg. Impulse can be calculated as . First you need to verify the numerical accuracy of the Impulse-Momentum theorem. Then you should report your best estimate of the impulse for cia_file1.txt from 1.15 s to 1.45 s and again from 1.95 s to 2.55 s. (Note: Smoothing has a dramatic affect on the accuracy of integrating data. Please vary the smoothing window to find the best result.) Can you draw any conclusions about the impulse that would interest the CIA?
Page 4