# Solving the following system of integrodifferential equations: speed up of the code

Consider a system of equations $$\begin{cases}\frac{\partial f(E,t)}{\partial t}-H[T_{\gamma},f(E,t)]E\frac{\partial f(E,t)}{\partial E} -I[f(E,t),E,T_{\gamma}(t)]=0, \\ \frac{\partial T_{\gamma}}{\partial t} + \frac{3H[T_{\gamma},f(E,t)](4\rho_{\gamma}/3 + \rho_{e}+ p_{e})+\int \limits_{0}^{\infty} dE\ I[f(E,t),E,T_{\gamma}(t)]}{\partial \rho_{\gamma}(T_{\gamma})/\partial T_{\gamma} + \partial \rho_{e}(T_{\gamma})/\partial T_{\gamma}}=0\end{cases}$$ Here: $$I[f(E,t),E,T_{\gamma}(t)] \equiv\int \limits_{\Gamma(E)} dE' G_{1}[E',E,f(E,t),f(E',t),T_{\gamma}(t)]+ G_{2}(f(E,t),T_{\gamma}(t),E)$$ is some integral including $$f(E)$$, $$f(E')$$ and $$f(E)f(E')$$ terms, with $$\Gamma(E)$$ being the domain of integration. $$H,\rho_{e},p_{e},\rho_{\gamma}$$ are known functions; in particular, $$H = b \sqrt{\frac{\pi^{2}}{30}g(T_{\gamma}) + 6\int \limits_{0}^{\infty}f(E)\frac{E^{3}}{2\pi^{2}}dE}, \quad b = \text{const},$$ while $$\rho_{\gamma},p_{e},\rho_{e}$$ do not include $$f$$.

Formally, $$E$$ varies in the limits $$0, but below I restrict myself by the domain $$0.

The implementation of the solution of this system in Mathematica for a particular $$I$$, given below, follows an approach from this question:

(*Definitions*)
ME = 0.5;
etaBPlanck = 6.1*10^-10;
ng[T_] = 2*Zeta[3]/Pi^2 T^3;
mPl = 1.2*10^22;
hbar = 6.582119*10^-25;
sToGeVminusOne = 1/hbar;
sToMeVminusOne = sToGeVminusOne*10^-3;
rhoeTemp[T_] :=
4/(2*Pi^2)NIntegrate[(Ee^2 Sqrt[Ee^2 - ME^2])/(
Exp[Ee/T] + 1), {Ee, ME, Infinity}]
peTemp[T_] :=
4/(6*Pi^2)NIntegrate[(Ee^2 - ME^2)^(3/2)/(Exp[Ee/T] + 1), {Ee, ME, Infinity}]
rhoe[T_] =
Quiet[10^Interpolation[
Join[Table[{10^T, Log10[etaBPlanck*ng[10^T]*ME]}, {T,
Log10[10^-8], Log10[0.0049], 0.1}],
Table[{10^T,
Log10[etaBPlanck*ng[10^T]*ME + rhoeTemp[10^T]]}, {T,
Log10[0.005], Log10[19.99], 0.01}],
Table[{T, Log10[7/8*4*Pi^2/30*T^4]}, {T, 20, 300, 0.5}]],
InterpolationOrder -> 1][T]]
pe[T_] = Quiet[
10^Interpolation[
Join[Table[{10^T, -90}, {T, Log10[10^-8], Log10[0.0049], 0.1}],
Table[{10^T, Log10[10^-99 + peTemp[10^T]]}, {T, Log10[0.005],
Log10[19.99], 0.01}],
Table[{T, Log10[1/3 7/8*4*Pi^2/30*T^4]}, {T, 20, 300, 0.5}]],
InterpolationOrder -> 1][T]]
DrhoeDT[T_] = D[rhoe[T], T];
rhog[T_] = 2*Pi^2/30 T^4;
DrhogDT[T_] = D[rhog[T], T];

(*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}]

(*Integral I - discretized*)
imax = 100;
DEvalue = 0.5;
Emin = 0.01;
EStep[k_, DE_] = Emin + DE*k;
fnThermalDiscrete[k_, Tn_] = Exp[-EStep[k, DEvalue]/Tn];
dS11dE2discrete[i_, j_] =
dS11dE2[E1, E2] /. {fn[E1] -> fn[i][t],
fn[E2] -> fn[j][t]} /. {E1 -> EStep[i, DEvalue],
E2 -> EStep[j, DEvalue]};
S12discrete[i_] =
S12[E1, T] /. {fn[E1] -> fn[i][t], fn[E2] -> fn[j][t]} /. {E1 ->
EStep[i, DEvalue], E2 -> EStep[j, DEvalue]};
S21discrete[i_] =
S21[E1, T] /. {fn[E1] -> fn[i][t], fn[E2] -> fn[j][t]} /. {E1 ->
EStep[i, DEvalue], E2 -> EStep[j, DEvalue]};
dS221dE2discrete[i_, j_] =
dS221dE2[E1, E2, T] /. {fn[E1] -> fn[i][t],
fn[E2] -> fn[j][t]} /. {E1 -> EStep[i, DEvalue],
E2 -> EStep[j, DEvalue]};
dS222dE2discrete[i_, j_] =
dS222dE2[E1, E2, T] /. {fn[E1] -> fn[i][t],
fn[E2] -> fn[j][t]} /. {E1 -> EStep[i, DEvalue],
E2 -> EStep[j, DEvalue]};
CollisionIntegrali[
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]}

(*A derivative dfn/dE - discretized*)
fnEnDerivative[i_] =
If[i != imax, (fn[i + 1][t] - fn[i][t])/DEvalue, -fn[imax][t]/
DEvalue];
(*A function H*)
Hdynamical =
sToMeVminusOne/mPl Sqrt[
8*Pi/3 (rhog[Tg[t]] + rhoe[Tg[t]] +
6*Sum[fn[i][t] 4*Pi (EStep[i, DEvalue])^3/(2*Pi)^3, {i, 0, imax,
1}])];

(*System of equations and its solution*)
EquationsWithHubbleTable =
Join[{Tg'[t] +
1/(DrhoeDT[Tg[t]] +
DrhogDT[Tg[t]]) (4*Hdynamical*rhog[Tg[t]] +
3*Hdynamical*(rhoe[Tg[t]] + pe[Tg[t]]) +
DEvalue*Sum[
CollisionIntegrali[i] 4*Pi EStep[i, DEvalue]^3/(
2*Pi^3), {i, 0, imax, 1}]) == 0},
Table[fn[i]'[t] - Hdynamical*EStep[i, DEvalue]*fnEnDerivative[i] -
CollisionIntegrali[i] == 0, {i, 0, imax, 1}]];
H[T_] = 1.66 Sqrt[10.75] T^2/mPl sToMeVminusOne;
t0 = 0.05;
Tt[t_] = Quiet[T /. Solve[H[T] == 1/(2 t), T][[2]]];
InitialConditionTable =
Join[{Tg[t0] == Tt[t0]},
Table[fn[i][t0] == 1/(Exp[EStep[i, DEvalue]/(Tt[t0])] + 1), {i, 0,
imax, 1}]];
functionsTable = Join[{Tg}, Table[fn[i], {i, 0, imax, 1}]];
tmax = 4;
solWithDynamicalHubble = NDSolve[{EquationsWithHubbleTable, InitialConditionTable}, functionsTable, {t, t0, tmax}, Method -> {"EquationSimplification" -> "Solve"}][[1]];


My question is the following: is it possible to speed up the code above? In my machine, it requires 400 seconds or so to evaluate solWithDynamicalHubble for imax = 100 (apart from computing the table of equations), whereas my goal is to solve for imax = 20000.

P.S. My machine is 16 Gb RAM, i7-10875H. When running Mathematica, only ~15% of the CPU is used.

• What do you suppose to get from this computation? Commented Oct 29, 2020 at 18:49
• @AlexTrounev : maybe, a dataset {time, {{E1, f(E1,time)},{E2,f(E2,time)},...}} with some handled time-step. Commented Oct 29, 2020 at 20:09
• Parameter Tgin is not defined. Is it Tg[tin]`? Commented Oct 30, 2020 at 10:59
• @AlexTrounev : yes. Thanks, I have corrected the typo. Commented Oct 30, 2020 at 11:52
• For the first equation see my last code on mathematica.stackexchange.com/questions/230366/… And second equation can be solve with same method. See also my answer on mathematica.stackexchange.com/questions/230346/… Commented Oct 30, 2020 at 12:28