Numerical Integration

% of readers think this story is Fact. Add your two cents.

A homework problem in Chapter 14 of Intermediate Physics for Medicine and Biology states

Problem 28. Integrate Eq. 14.33 over all wavelengths to obtain the Stefan-Boltzmann law, Eq. 14.34. You will need the integral

Equation 14.33 is Planck’s blackbody radiation law and Eq. 14.34 specifies that the total power emitted by a blackbody.

Suppose Russ Hobbie and I had not given you that integral. What would you do? Previously in this blog I explained how the integral can be evaluated analytically and perhaps you’re skilled enough to perform that analysis yourself. But it’s complicated, and I doubt most scientists could do it. If you couldn’t, what then?

You could integrate numerically. Your goal is to find the area under the curve shown below.

Unfortunately x ranges from zero to infinity (the plot shows the function up to only x = 10). You can’t extend x all the way to infinity in a numerical calculation, so you must either truncate the definite integral at some large value of x or use a trick.

A good trick is to make a change of variable, such as

When x equals zero, t is also zero; when x equals infinity, t is one. The integral becomes

Although this integral looks messier than the original one, it’s actually easier to evaluate because the range of t is finite: zero to one. The integrand now looks like this: 

 

The colored stars in these two plots are to guide the reader’s eye to corresponding points. The blue star at t = 1 is not shown in the first plot because it corresponds to x = ∞.

We can evaluate this integral using the trapezoid rule. We divide the range of t into N subregions, each extending over a length of Δt = 1/N. Ordinarily, we have to be careful dealing with the two endpoints at t = 0 and 1, but in this case the function we are integrating goes to zero at the endpoints and therefore contributes nothing to the sum. The approximation is shown below for N = 4, 8, and 16. 

The area of the purple rectangles approximates the area under the red curve This approximation gets better as N gets bigger. In the limit as N goes to ∞, you get the integral.

I performed the calculation using the software Octave (a free version of Matlab). The program is: