Per determinare la capacità termica di un solido usando un metodo basato sulle frequenze di vibrazione del cristallo (vedi) si deve calcolare ∫ [0, k] x4·ex/(ex−1)2 dx.  Calcola in modo approssimato (usando del software) il valore di questo integrale quando k = 1.

Prima di fare calcoli mi conviene tracciare il grafico della funzione, prolungandola per continuità in 0 (per x → 0 x4·ex/(ex−1)2 → 0 in quanto ex−1 tende a 0 come x).
Dal grafico capisco che l'area compresa tra l'asse x, il grafico e la retta x = 1 è leggermente inferiore a quella di un triangolo di base 1 e altezza 0.8, ossia a 0.4.
Ora facciamo i calcoli. Se usiamo R (vedi). 		   

t = function(x) x^4*exp(x)/(exp(x)-1)^2
integral(t,0,1)
# 0.317244
# ovvero, se volessi avere più cifre:
more( integral(t,0,1) )
# 0.317244045234427
#
# Potrei usare anche il comando per calcolare l'area dei poligoni.
# Devo prolungare t dandole il valore 0 in 0
t = function(x) ifelse(x==0,0,x^4*exp(x)/(exp(x)-1)^2)
n=1e3; a=0; b=1; x=seq(a,b,len=n); y=t(x); x=c(x,b,a); y=c(y,0,0); areaPol(x,y)
# 0.3172442
n=1e4; a=0; b=1; x=seq(a,b,len=n); y=t(x); x=c(x,b,a); y=c(y,0,0); areaPol(x,y)
# 0.317244
# ...

Potrei usare più semplicemente uno script online (vedi integ, avendo preso pow(x,4)*exp(x)/pow(exp(x)-1,2) come "F(x)"):

0.3172440452344334  if a=0 b=1 n=1e7 [8.198997036856781e-14]
0.3172440452343514  if a=0 b=1 n=1e6 [6.966816012976551e-12]
0.3172440452273846  if a=0 b=1 n=1e5 [6.972915578273842e-10]
0.31724404453009303 if a=0 b=1 n=1e4 [0.31724404453009303]
0.31724404523443

    Ovvero posso usare WolframAlpha. Vedi qui.
integrate x^4*exp(x)/(exp(x)-1)^2 for x from 0 to 1

    Potrei fare anche il grafico:
plot x^4*exp(x)/(exp(x)-1)^2, 0 < x < 1