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³/((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.