5

Planck's law dependent on frequency rho is as follows:

   B[T_, h_, rho_, k_, c_] := (2 h rho^3)/c^2 1/(Exp[h rho/(T k)] - 1)

As you can see, the denominator can only be zero for rho = 0. All the rest are nonzero constants. Curiously if we try to plot it

Plot[B[5800, 6.626 10^-34, rho, 1.38 10^-23, 299792458], {rho, 
  0.1 10^-6, 3 10^-6}]

It will return "Infinite expression 1/0. encountered" and the plot will be blank. In fact, substituting any value for rho appears to result in this error (out = ComplexInfinity). I have no experience with precision, but I don't see any other possible reason. Anyway to solve this?

C. E.
  • 70,533
  • 6
  • 140
  • 264
  • 1
    when doing numerical calculations, it's better to work in dimensionless units, like here: http://mathematica.stackexchange.com/a/2057/16 so as to avoid underflow/overflow problems. so if I were you I'd set $h=k_B=c=1$ and go from there – acl May 17 '12 at 22:43
  • As @acl said, try Plot[B[5, 1, rho, 1, 1], {rho, 0, 100}] ... Plank's units – Dr. belisarius May 17 '12 at 22:45
  • Also, not sure what you're doing but there is not much 5800K black-body radiation at 10^-6 Hertz. A better place to look might be around 10^15 Hz. – dws May 17 '12 at 22:51

1 Answers1

11

To understand what is going on, look at this:

N[b[5800, 6626 10^-37, rho, 138 10^-25, 299792458]]

Mathematica graphics

look at the exponent in the second term in the denominator. If you put $\rho\approx 10^{-6}$, you're trying to calculate the difference between two numbers that are extremely close to each other, namely, between $1$ and $\exp(\alpha)$ with $\alpha\approx 10^{-21}$. So you run into numerical accuracy problems with machine numbers; see here and at the end of this answer.

In general, when doing numerical calculations for physical problems, it's better to work in dimensionless units, so as to avoid underflow/overflow problems. so if I were you I'd set $h=k_B=c=1$ (this amounts to a choice of units) and go from there.

Community
  • 1
acl
  • 19,834
  • 3
  • 66
  • 91