# Discretization of the following integral: an optimal way

Consider an integral

$$I[f_{n},T]=\int \limits_{0}^{\infty} dE_{n}F(f_{n}(E_{n}),E_{n},T)+\int\limits_{0}^{\infty} dE_{n}\int_{\mathcal{F}(E_{n})}dE_{n}^{'}G(f_{n}(E_{n}),f_{n}(E_{n}^{'}),E_{n},E_{n}^{'},T),$$ where $$T>0$$ is some parameter, while $$f_{n}$$ is the distribution function. The integral has the property that for $$f_{n} = e^{-E_{n}/T_{n}}$$, it 1) is negative if $$T_{n}>T$$, 2) vanishes at $$T_{n} = T$$, 3) positive if $$T.

The integral $$I$$ enters some partial integrodifferential equation on $$f_{n}(E_{n},t)$$ which I want to solve as a set of ODEs with the help of the discretization of $$E$$ domain (see this question). Therefore, I need to discretize this integral.

The integral in continuous and discrete forms, $$f_{n} = e^{-E_{n}/T_{n}}$$, is given below (the final functions are CollisionIntegralIntegratedBoltzmannContinuous[T, Tn], CollisionIntegralIntegratedBoltzmannDiscretized[T, Tn]). To discretize it, I have just used $$E_{n} \to E_{\text{min}}+i\Delta E$$, where $$0\leqslant i \leqslant imax$$, and $$\int dE_{n} \to \Delta E \sum_{i}$$.

(*Some parameters*)
hbar = 6.582119*10^-25;
sToGeVminusOne = 1/hbar;
sToMeVminusOne = sToGeVminusOne*10^-3;
(*Discretization of En values*)
imax = 100;
DEvalue = 0.25;
Emin = 0.01;
Estep[k_] = Emin + DEvalue*k;
fnBoltzmannDiscrete[k_, Tn_] = Exp[-Estep[k]/Tn];
fnThermalBoltzmann[E\[Nu]_, Tn_] = Exp[-E\[Nu]/Tn];

(*Integral I-continuous*)
dS11dE2[E1_, E2_] = -0.0034976546222801287 E1 E2^3*fn[E1] fn[E2];
S12[E1_, T_] = 0.02098592773368077 Exp[-E1/T] E1 T^4;
S21[E1_, T_] = -0.08394371093472308 E1 T^4 fn[E1];
dS221dE2[E1_, E2_, T_] =
0.03147889160052116/(E1^2) Exp[-E1/
T] T^2 (-E1^2 (E2^2 + 2 E2 T - 2 (-1 + Exp[E2/T]) T^2) -
2 E1 T (E2^2 + 2 E2 T - 2 (-1 + Exp[E2/T]) T^2) +
2 T^2 ((-1 + Exp[E2/T]) E2^2 - 2 (1 + Exp[E2/T]) E2 T +
4 (-1 + Exp[E2/T]) T^2)) fn[E2];
dS222dE2[E1_, E2_, T_] =
0.03147889160052116/(E1^2) Exp[-E1/
T] T^2 (2 (-1 + Exp[E1/T]) T^2 (E2^2 + 2 E2 T + 4 T^2) -
E1^2 (E2^2 + 2 E2 T - 2 (-1 + Exp[E1/T]) T^2) -
2 E1 T (E2^2 + 2 E2 T + 2 (1 + Exp[E1/T]) T^2)) fn[E2];
CollisionIntegralContinuous[E1_, T_] :=
S12[E1, T] + S21[E1, T] +
NIntegrate[dS11dE2[E1, E2], {E2, 0, Infinity}] +
NIntegrate[dS221dE2[E1, E2, T], {E2, 0, E1}] +
NIntegrate[dS222dE2[E1, E2, T], {E2, E1, Infinity}]
CollisionIntegralIntegratedBoltzmannContinuous[T_, Tn_] :=
NIntegrate[((S12[E1, T] + S21[E1, T]) /. {fn[E1] ->
fnThermalBoltzmann[E1, Tn]}) E1^3/(2*Pi^2), {E1, 0,
Infinity}] +
NIntegrate[(dS11dE2[E1,
E2] /. {fn[E1] -> fnThermalBoltzmann[E1, Tn],
fn[E2] -> fnThermalBoltzmann[E2, Tn]}) E1^3/(2*Pi^2), {E1, 0,
Infinity}, {E2, 0, Infinity}] +
NIntegrate[(dS221dE2[E1, E2,
T] /. {fn[E1] -> fnThermalBoltzmann[E1, Tn],
fn[E2] -> fnThermalBoltzmann[E2, Tn]}) E1^3/(2*Pi^2), {E1, 0,
Infinity}, {E2, 0, E1}] +
NIntegrate[(dS222dE2[E1, E2,
T] /. {fn[E1] -> fnThermalBoltzmann[E1, Tn],
fn[E2] -> fnThermalBoltzmann[E2, Tn]}) E1^3/(2*Pi^2), {E1, 0,
Infinity}, {E2, E1, Infinity}]

(*Integral I-discretized*)
dS11dE2discrete[i_, j_] =
dS11dE2[E1, E2] /. {fn[E1] -> fn[i][t],
fn[E2] -> fn[j][t]} /. {E1 -> Estep[i], E2 -> Estep[j]};
S12discrete[i_] =
S12[E1, T] /. {fn[E1] -> fn[i][t], fn[E2] -> fn[j][t]} /. {E1 ->
Estep[i], E2 -> Estep[j]};
S21discrete[i_] =
S21[E1, T] /. {fn[E1] -> fn[i][t], fn[E2] -> fn[j][t]} /. {E1 ->
Estep[i], E2 -> Estep[j]};
dS221dE2discrete[i_, j_] =
dS221dE2[E1, E2, T] /. {fn[E1] -> fn[i][t],
fn[E2] -> fn[j][t]} /. {E1 -> Estep[i], E2 -> Estep[j]};
dS222dE2discrete[i_, j_] =
dS222dE2[E1, E2, T] /. {fn[E1] -> fn[i][t],
fn[E2] -> fn[j][t]} /. {E1 -> Estep[i], E2 -> Estep[j]};
CollisionIntegralDiscrete[
i_] := (S12discrete[i] + S21discrete[i] +
DEvalue*Sum[dS11dE2discrete[i, j], {j, 0, imax, 1}] +
DEvalue*Sum[dS221dE2discrete[i, j], {j, 0, i, 1}] +
If[i == imax, 0,
DEvalue*Sum[
dS222dE2discrete[i, j], {j, i + 1, imax, 1}]]) /. {T -> Tg[t]}
CollisionIntegralIntegratedBoltzmannDiscrete[Tg_,
Tn_] := (DEvalue*
Sum[((S12discrete[i] + S21discrete[i]) /. {fn[i][t] ->
fnBoltzmannDiscrete[i, Tn]}) (Estep[i])^3/(2*Pi^2), {i, 0,
imax, 1}] +
DEvalue^2*
Sum[(dS11dE2discrete[i,
j] /. {fn[i][t] -> fnBoltzmannDiscrete[i, Tn],
fn[j][t] -> fnBoltzmannDiscrete[j, Tn]}) (Estep[i])^3/(
2*Pi^2), {i, 0, imax, 1}, {j, 0, imax, 1}] +
DEvalue^2*
Sum[(dS221dE2discrete[i,
j] /. {fn[i][t] -> fnBoltzmannDiscrete[i, Tn],
fn[j][t] -> fnBoltzmannDiscrete[j, Tn]}) (Estep[i])^3/(
2*Pi^2), {i, 0, imax, 1}, {j, 0, i, 1}] +
DEvalue^2*
Sum[If[i == imax, 0,
Sum[(dS222dE2discrete[i,
j] /. {fn[i][t] -> fnBoltzmannDiscrete[i, Tn],
fn[j][t] -> fnBoltzmannDiscrete[j, Tn]}) (Estep[i])^3/(
2*Pi^2), {j, i + 1, imax, 1}]], {i, 0, imax, 1}]) /. {T -> Tg}


Comparing the values of CollisionIntegralIntegratedBoltzmannContinuous[T, Tn], CollisionIntegralIntegratedBoltzmannDiscretized[T, Tn] at $$T = T_{n} = 2$$, $$T = T_{n} = 3$$ for DEvalue = 0.25 and imax = 100, I find

Timing[CollisionIntegralIntegratedBoltzmannDiscrete[3, 3]]
CollisionIntegralIntegratedBoltzmannContinuous[3, 3]


{0.941387,-43.7787}

-9.28833*10^-6

This is a very poor agreement. It may be improved by increasing imax and decreasing DEvalue. Namely, for DEvalue = 0.05 and imax = 1000`, I get

{98.1421, -0.182558}

The agreement is somewhat better now. However, this improvement is accompanied by a significant increase in timing. This problem suggests that the way of the discretization that is used in my code is poor.

Could you please tell me how to discretize more efficiently, assuming that it may be used for "arbitrary" $$f_{n}$$ (in the sense that this integral $$I$$ is evaluated for dynamically changed $$f_{n}$$)?