integral method gives an answer of zero

Hello,

I have been trying compute

where β is a fixed constant. My code consists of the following:

k = 1.380649*10^(-23);
t = 5800;

beta = (k*t)^(-1);

fun = @(x) x.^3 ./ (exp(beta .* x) - 1);

q = integral(fun,0,Inf)

This results in the integral being evaluate as 0. However, I know that this integral is infact not zero over the domain

Am I making somekind of mistake?

Answer 1

[Walter Roberson]

The values are so small relative to your constants that they might as well be zero.

With your beta being 1/(k*t) you are ending up with some large exponents in the exp() -- large enough that you are going to overflow numerically.

syms x
Q = @(v) sym(v);

k = Q(1380649)*Q(10)^(-29);
t = Q(5800);
beta = (k*t)^(-1);
fun = x.^3 ./ (exp(beta .* x) - 1)

fun = 

int(fun, 0, Inf)

ans = 

vpa(ans)

ans = 

2.6702586874891152137386272298896e-76