Analytical evaluation of Fermi-Dirac integrals

edit

Maxima (i.e. Sage) is perfectly able to find a primitive of your function,. For your second example, it gives :

sage: f(x)=x^2/(1+e^x)
sage: f(x).integrate(x)
sage: g(x)=f(x).integrate(x)
sage: g
x |--> 1/3*x^3 - x^2*log(e^x + 1) - 2*x*dilog(-e^x) + 2*polylog(3, -e^x)

This can be easily checked:

sage: g(x).diff(x).factor()
x^2/(e^x + 1)

Sage can find the limit at x=0...

sage: g(x).limit(x=0)
-3/2*zeta(3)

... but not at x=$\infty$:

sage: g(x).limit(x=oo)
limit(1/3*x^3 - x^2*log(e^x + 1) - 2*x*dilog(-e^x) + 2*polylog(3, -e^x), x, +Infinity)

Mathematica agrees with Sage for x=0...

sage: mathematica("Limit[x^3/3-x^2*Log[1+Exp[x]]-2*x*PolyLog[2,-Exp[x]]+2*PolyLo
....: g[3,-Exp[x]],x->0]")
(-3*Zeta[3])/2

And gives an answer for x=$\infty$:

sage: mathematica("Limit[x^3/3-x^2*Log[1+Exp[x]]-2*x*PolyLog[2,-Exp[x]]+2*PolyLo
....: g[3,-Exp[x]],x->Infinity]")
0

A quick numeric check (by comparison to numerical integration) lets us think that Mathamatica's answer is reasonable:

sage: numerical_integral(f(x),0,oo):
(1.8030853547394121, 5.780826801595601e-08)
sage: g(x).limit(x=0).n()
-1.80308535473939

Therefore, this is a bug. I'll check if this is already reported...

BTW, Sympy gives a wrong answer here :

sage: g(x).limit(x=oo, algorithm="sympy")
-Infinity

...but fails to find a primitive anyway:

sage: f(x).integrate(x, algorithm="sympy")
integrate(x^2/(e^x + 1), x)