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 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:

N=8; 

dt=1/N; 

s=0; 

for i=1:N-1 

     t=i*dt; 

     s=s+dt*t^3/((exp(t/(1-t))-1)*(1-t)^5); 

     endfor

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

N	I	% error
2	1.1640	–82
4	6.2823	–3.26
8	6.6911	3.04
16	6.5055	0.178
32	6.4940	0.000282
64	6.4939	0.00000235
128	6.4939	0.00000000174

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 π⁴/15. I wonder how many people could calculate 6.4939 and then say “Hey, I know that number; It’s π⁴/15”!