Without seeing your code, it is hard to know where you went wrong. I did write a few lines of code myself to see what might be going on - and found that the values converge quite nicely as you approach higher temperatures (for the lowest temperatures of 500 K and 1000 K, significant amounts of energy do indeed extend beyond 10,000 nm).
updated
There is a small problem with scaling: the SB law gives radiation from a black body in W/m$^2$ which is different than Planck's Law which gives radiance per unit area, per steradian.
Now there is a modification for the SB law: if you divide the usual equation by $\pi$ you get the result in "per steradian" (see the Wiki article).
This gives me the following code snippet (Python - not my strongest language but it's good for this kind of thing because of the scipy.integrate functions which do a lot of the hard work; the built in physical constants and plotting routines are nice too):
# integrate Planck's Law and compare results with Stefan-Boltzman
from math import pow,exp,pi
from scipy.constants import codata
from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
# get physical constants we need:
D = codata.physical_constants
h = D['Planck constant'][0]
k = D['Boltzmann constant'][0]
c = D['speed of light in vacuum'][0]
s = D['Stefan-Boltzmann constant'][0]
def planck(l, T):
# reference https://en.wikipedia.org/wiki/Planck's_law
if (l <= 0):
return 0 # no valid values <=0
p = c*h/(k*l*T)
if (p > 700):
return 1e-99 # preventing under/overflow in the exp function
else:
return (2*h*c*c)/(pow(l, 5.0) * (exp(c*h/(k*l*T))-1)) # <<< fixed
def SB(T, emissivity = 1.0):
# Stefan-Boltzmann law : returns radiance = total power per steradian
# reference https://en.wikipedia.org/wiki/Stefan-Boltzmann_law
return emissivity * s * pow(T, 4.0) / pi # <<< expression "per steradian"
# T values to use:
Tvec = np.arange(500,10500, 500)
Planck_integrals=[]
SB_integrals=[]
max_wavelength = 10000*1e-9 # 10,000 nm
for T in Tvec:
y,abserr = integrate.quad(planck, 0, max_wavelength, args=(T))
Planck_integrals.append(y)
SB_integrals.append(SB(T))
plt.figure()
plt.plot(Tvec, np.array(Planck_integrals)/np.array(SB_integrals))
plt.xlabel('Temperature (K)')
plt.ylabel('Ratio integral/SB value')
plt.title('Comparing direct integration and Stefan-Boltzmann - integrating 0-10 um')
plt.show()
And this plot:
I played around a bit with discrete integration to see if I can get the results to match (so you don't need to use existing libraries... although that's usually a really good idea). This was the result (using the same definition of functions planck and SB as before):
# stepwise integration
n = 1000
simpson = 0.
trapezoid = 0.
max_wavelength = 10000.0
min_wavelength = 0.
H = (max_wavelength - min_wavelength)/n
wavelength_scale = 1e-9
T = 2000.0
for i in range(n):
if ((i==0) | (i == n-1)):
w_trap = 1
w_simpson = 1
else:
w_trap = 2
if (i%2==0):
w_simpson = 4
else:
w_simpson = 2
wavelength = (min_wavelength + H * i) * wavelength_scale
hp = H * planck(wavelength, T) * wavelength_scale
trapezoid = trapezoid + hp * w_trap
simpson = simpson + hp * w_simpson
print 'trapezoid: ', trapezoid / 2
print 'simpson: ', simpson / 3
print 'SB: ', SB(T)
This gave the following output:
trapezoid: 284606.470099
simpson: 284606.473505
SB: 288789.72548
Clearly the two integration methods give nearly the same result; and it's very comparable with the result from the Stefan-Boltzman formula. I have to conclude the error is in your integration formula (as I pointed out in my comment), and/or in extremely coarse sampling of the function. But even when I set n=100 in the above I get comparable results...
Best Answer
The solution to your problem(1) is this : \begin{align} R(T) & =\frac c4 \int\limits_{\lambda=0}^{\lambda=\infty} \lambda^{-5}f(\lambda T)\mathrm{d}\lambda=\frac c4 \int\limits_{\lambda=0}^{\lambda=\infty} T^{4}\dfrac{f(\lambda T)}{(\lambda T)^{5}}\mathrm{d}(\lambda T)\qquad \Longrightarrow \nonumber\\ R(T) & =\frac c4 \underbrace{\left(\:\:\int\limits_{\mu=0}^{\mu=\infty} \dfrac{f(\mu)}{\mu^{5}}\mathrm{d}\mu\right)}_{A=\text{constant}}T^{4}= \underbrace{\left(\frac c4 A\right)}_{\sigma}\, T^{4}=\sigma\, T^{4} \tag{01} \end{align}
But, sincerely, trying to find this directly you are missing important facts about the "before Planck" adventure of the blackbody radiation theory. For example, you must try to find from Wien's Law (your first equation) why if you know the function $\;\rho(\lambda,T_{1})\;$ for a given temperature $\;T_{1}\;$ then you know it for any temperature $\;T\;$ or that $\;\lambda_{\rm max}\cdot T=b=\rm constant\;$ (Wien's Displacement Law), see Emilio Pisanty answer therein : Showing Wien's Displacement Law from Wien's Law.
(1) "Quantum Mechanics" B. H. Bransden-C. J. Joachain, 2nd Edition 2000, Pearson Education Limited (Problem 1.3, page 45)