# Numerical Integration

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 theStefan-Boltzmann law, Eq. 14.34. You will need theintegral

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:

N=8;

dt=1/N;

s=0;

for i=1:N-1

t=i*dt;

s=s+dt*t³/((exp(t/(1-t))-1)*(1-t)⁵);

endfor

I found the results shown below. The error is the difference between the numerical integration and the exact result (π^4/15 = 6.4939…), divided by the exact result, and expressed as a percent difference.

These results show that you can evaluate the integral accurately without too much effort. You could even imagine doing this by hand if you didn’t have access to a computer — using, say, *N* = 16 — and getting an answer accurate to better than two parts per thousand.

For many purposes, a numerical solution such as this one is adequate. However, 6.4939… doesn’t look as pretty as π^4/15. I wonder how many people could calculate 6.4939 and then say “Hey, I know that number; It’s π^4/15”!

*Originally published at **http://hobbieroth.blogspot.com**.*