Numerically integrating the Planck distribution

[heafnerj]

I'm working on a project that requires me to numerically integrate the Planck spectral distribution. The object is to find the median wavelength, with exactly 50% of the radiance on either side. I'm using a standard composite Simpson rule method and I get good convergence with temps around 5000K. For temps at and above 7000K, I can't get convergence. Is there a better scheme to use for this application?

 

[Hurkyl]

How familiar are you with issue of numeric stability? Have you analyzed the error term for Simpson's rule?

 

[heafnerj]

How familiar are you with issue of numeric stability? Have you analyzed the error term for Simpson's rule?

Rudimentary. No.

 

[f95toli]

Sometimes the simplest methods are the best when you are trying to solve numerical problems. When I am integrating something that I suspect might "misbehave" I always use methods based on polynomials of order 1; e.g. the trapezoidal rule with very small intervals.
This is very inefficient but it always works; and unless you need to integrate thousands of equations you will rarely notice the difference on a modern computer.

Of course you can always try to use e.g. the composite Simpson's rule; it should be better than the "plain" Simpson's rule for your problem.