# How to Optimize Integro-Differential Equation Solving Time in Mathematica?

I'm trying to solve an integro-differential equation numerically in Mathematica that models the evolution of a physical system . However, the computation takes a very long time and doesn't yield a result. The main equation involves calculating an integral that depends on the solution itself.

The Mathematica code snippet I'm using is as follows:

y1e[T_] :=
Exp[-0.848259 + 0.685952 (1 + Tanh[(Log[T] + 1.76853)/0.174619])];
y2e[T_] :=
Exp[-2.36415 + 1.55878 (1 + Tanh[(Log[T] + 1.2981)/0.656795])];
y3e[T_] :=
Exp[-1.52151 + 1.50891 (1 + Tanh[(Log[T] + 1.64656)/1.5446])];
y4e[T_] :=
Exp[-0.792215 + 0.464922 (1 + Tanh[(Log[T] + 3.63313)/0.572627])];
ye[T_] := y1e[T] + y2e[T] + y3e[T] + y4e[T];

(*the cosmic constant*)
GF = 1.166*10^-5 (*GeV^-2*);
Ba = 10.88*10^-9(*GeV^-4*);
Mpl = 2.4*10^18(*GeV*);
Lini = 0.001;(*initial asymmetry*)
LogSpace[start_, end_, steps_] :=
Exp[Subdivide[Log[start], Log[end], steps - 1]];

initialT = 1(*GeV*);
finalT = 0.001(*GeV*);
Ti = 1(*reference T GeV*);

data = {{1.1733446880479033*^-6,
3.924573951846154}, {1.575801106756097*^-6,
3.925551384983876}, {2.116299799494522*^-6,
3.9265290615558004}, {2.842189171043508*^-6,
3.92750698162256}, {3.817058096365375*^-6,
3.928485145244796}, {5.126306390675286*^-6,
3.9294635524831656}, {6.884625946904336*^-6,
3.9304422033983455}, {9.246047898152396*^-6,
3.9314210980510236}, {0.000012417435949932553,
3.9324002365019015}, {0.000016676607916068563,
3.933379618811702}, {0.000022396672928906124,
3.9343592450411573}, {0.000030078716295840698,
3.94550039512886}, {0.00004039569523909017,
3.9659153370897418}, {0.00005425139084393372,
4.011479839380392}, {0.00007285958048453544,
4.251914278258622}, {0.00009274189139404798,
4.672634609257866}, {0.00011039695561422658,
5.193076590174616}, {0.00012793665511201275,
5.746906216145927}, {0.00014628886148662263,
6.295139137829342}, {0.00016727364785577256,
6.877333089489969}, {0.00019384982924560326,
7.548040826736746}, {0.0002338665979378173,
8.32594808859626}, {0.00029848451526735364,
9.144850067018517}, {0.00038817332806409977,
9.783338674014455}, {0.0005157558696968545,
10.249862235177805}, {0.0006926597772700497,
10.516360798861012}, {0.0009302416031247371,
10.67906453764181}, {0.001249313831380024,
10.768646915037799}, {0.0016778276138528598,
10.807100636807162}, {0.002253321328154625,
10.85768463992212}, {0.0030262089895260265,
10.864390394431076}, {0.004064196585663218,
10.8550928579246}, {0.0054582132113432155,
10.873807787817615}, {0.007330376578134934,
10.904599771224218}, {0.009844690688447049,
10.935478950051282}, {0.013221412804395764,
11.02721228482572}, {0.01775634827708484,
11.313922112519805}, {0.023846763489020602,
11.937678068614467}, {0.03202618691790199,
12.858371606490087}, {0.04301114694128756,
14.014311223452552}, {0.05776394067602663,
15.178809679263422}, {0.07757693248632162,
16.39471039108891}, {0.10007912788298898,
17.674823523413238}, {0.11913096522762875,
19.238690967483862}, {0.13440619971917453,
21.215741381186064}, {0.14372338451500524,
23.360392802572335}, {0.15164005838811956,
25.704185035424487}, {0.15576045120733759,
28.025795973106142}, {0.15999280413238237,
31.309159094593443}, {0.1643401593647731,
33.5428019239196}, {0.16971319445082872,
37.24056453278068}, {0.17455850927199254,
42.95113553120909}, {0.17930164664468345,
47.82131605521459}, {0.1889246483450124,
54.55531504185653}, {0.22670586046570043,
63.11025222146649}, {0.30446581423936964,
66.14953210869378}, {0.4088973785239546,
68.57312506673385}, {0.5491488973284621,
71.37415760938619}, {0.7375065903470914,
74.31697760223719}, {0.9904708421549656,
77.49524255279282}, {1.3302016578556475,
80.63101107448512}, {1.7864599090187747,
83.18583931560693}, {2.3992144256353862,
85.22301466905155}, {3.222143318815693,
86.63723747822668}, {4.327336254757254,
87.55732026907872}, {5.811609605440912,
87.9023547893092}, {7.804987691659681,
88.24874897678148}, {10.482093086556752,
88.72715678410984}, {14.077443785420886,
89.2739091518901}, {18.905997294170486,
90.5549752395357}, {25.390741326018023,
92.56778396421544}, {34.09974808806011,
95.32508801207295}, {45.79593816261045,
98.67211125565548}, {61.50391336550636,
102.06142950274938}, {82.59971322871769,
104.86946578859701}, {110.93135789458279,
106.84561369126395}, {148.98073713962523,
108.21924567340007}, {200.08102722187982,
109.00651599297225}, {268.70868155689107,
109.63783639370376}, {360.8755739941897,
109.78640790407033}, {453.2367777248329, 109.80753583983174}};
gs[T_] := Interpolation[data, T, InterpolationOrder -> 1];

Dgs[T_] := D[gs[T0], T0] /. T0 -> T;

p[epison_, T_] := epison*T*(gs[T]/gs[Ti])^(1/3);
episonp[epison_, T_] := epison*(gs[T]/gs[Ti])^(1/3);
Gama[epison_, T_] := ye[T]*GF^2*(gs[T]/gs[Ti])^(1/3)*epison*T^5;
stot[T_] := (2*\[Pi]^2)/45*gs[T]*T^3;
(*here and below the LveT denotes the Lve(T)*)
Va[epison_, T_, LveT_] :=
Sqrt[2]*GF*2*LveT*stot[T] - Ba*episonp[epison, T]*T^5;
Vabar[epison_, T_, LveT_] := -Sqrt[2]*GF*2*LveT*stot[T] -
Ba*episonp[epison, T]*T^5;
theta2M[sin2theta_, ms_, epison_, T_, LveT_] :=
1/4*sin2theta*((1 - (2*p[epison, T])/ms^2*
Va[epison, T, LveT])^2 + (p[epison, T]^2*Gama[epison, T]^2)/
ms^4)^(-1);
theta2Mbar[sin2theta_, ms_, epison_, T_, LveT_] :=
1/4*sin2theta*((1 - (2*p[epison, T])/ms^2*
Vabar[epison, T, LveT])^2 + (p[epison, T]^2*
Gama[epison, T]^2)/ms^4)^(-1);

fve[epison_, T_,
LveT_] := (1 +
Exp[episonp[epison, T] - (4*\[Pi]^2)/15*gs[T]*LveT])^(-1);
fvebar[epison_, T_,
LveT_] := (1 +
Exp[episonp[epison, T] + (4*\[Pi]^2)/15*gs[T]*LveT])^(-1);



and importantly, I use the iterative method:

(*The advised code*)
Off[InterpolatingFunction::dmval];
sin2theta = 10^-13;
ms = 2*10^-5;
n = 25;
sfvetest[0] =
Interpolation[
Table[{T, -(45/(4*Pi^4*gs[T]))*(-Sqrt[90/(Pi^2*gs[T])]*(Mpl/
T^3)*(1 + (T/3)*Dgs[T]/gs[T]))*
NIntegrate[
epison^2*
Gama[epison,
T]*(theta2M[sin2theta, ms, epison, T, Lini]*
fve[epison, T, Lini] -
theta2Mbar[sin2theta, ms, epison, T, Lini]*
fvebar[epison, T, Lini]), {epison, 0, 50},Method -> "LocalAdaptive"]}, {T, initialT,
finalT, -0.01}]];
Do[Lvetest[i] =
NDSolveValue[{D[Lve[T], T] == sfvetest[i - 1][T],
Lve[initialT] == Lini}, Lve, {T, initialT, finalT}(*,
Method\[Rule]*)];
sfvetest[i] =
Interpolation[
Table[{T, -(45/(4*Pi^4*gs[T]))*(-Sqrt[90/(Pi^2*gs[T])]*(Mpl/
T^3)*(1 + (T/3)*Dgs[T]/gs[T]))*
NIntegrate[
epison^2*
Gama[epison,
T]*(theta2M[sin2theta, ms, epison, T, Lvetest[i][T]]*
fve[epison, T, Lvetest[i][T]] -
theta2Mbar[sin2theta, ms, epison, T, Lvetest[i][T]]*
fvebar[epison, T, Lvetest[i][T]]), {epison, 0, 50},Method -> "LocalAdaptive"]}, {T, initialT,
finalT, -0.01}]];, {i, 1, n}];
Table[LogLogPlot[Lvetest[i][T], {T, initialT, finalT}], {i, 1, n}]
LogLogPlot[Lvetest[n][T], {T, initialT, finalT}]


This process involves nested numerical integrals and differential equations. I have tried various methods and options in NDSolveValue and NIntegrate, but I cannot seem to reduce the computation time to a reasonable duration.
And the result:

I suspect that the iterative method is unstable, so I've attempted an alternative approach using the summation method. In this method, I calculate $$\epsilon$$ as $$\frac{\epsilon_{max}}{Nbin}*i$$, summing from 1 to Nbin, I wish Nbin should be more than 2000, but in fact, it will evaluate slow ,I guess it maybe caused by the complexity of the equation. And the code follows:

episonmax = 20.;
Nbin = 500.;
episoni[i_] = episonmax/Nbin*i;

sfve[sin2theta_, ms_, T_,
LveT_] := -(45/(4*Pi^4*gs[T]))*(-Sqrt[90/(Pi^2*gs[T])]*(Mpl/
T^3)*(1 + (T/3)*Dgs[T]/gs[T]))*
ParallelSum[
episoni[i]^2*
Gama[episoni[i],
T]*(theta2M[sin2theta, ms, episoni[i], T, LveT]*
fve[episoni[i], T, LveT] -
theta2Mbar[sin2theta, ms, episoni[i], T, LveT]*
fvebar[episoni[i], T, LveT]), {i, 1, Nbin, 1}];

ffve[sin2theta_, ms_] :=
Block[{Lve, solution1},
solution1 =
NDSolve[{D[Lve[T], T] == sfve[sin2theta, ms, T, Lve[T]],
Lve[initialT] == Lini},
Lve, {T, initialT, finalT},(*AccuracyGoal->10,PrecisionGoal->10,
MaxStepSize->0.01,*)
Method -> {"EquationSimplification" -> "Residual"}];
Lve /. solution1[[1]]];


And the corresponding result:

LveT = ffve[10^-13., 2*10^-5.]; // AbsoluteTiming
LogLogPlot[LveT[T], {T, initialT, finalT}]
LveT[0.001]


Are there any solutions to improve the efficiency of the caculation when I increse the number of bins(Nbin)?

• It could be better to show this problem in Latex form as well, since your code is not so clear. Commented Apr 4 at 18:57
• Thank you for the suggestion; I have now included the LaTeX representation of the target equation. I apologize for any confusion caused. Commented Apr 4 at 21:39
• To improve computation time we can exclude NDSolve as well. Also summation can be vectorize. Tell me please, what code is valid - with NIntegrate or with Sum? Since results are not same. Commented Apr 7 at 11:22
• The code with Sum is Valid, The code with NIntegrate is the iterated method I used before, but I abondon it due to its instability. Commented Apr 7 at 11:28
• Are all definitions for functions and data are same for code with Sum as for iterative method? Commented Apr 7 at 11:37

To stop messages from NIntegrate use option Method -> "LocalAdaptive". To stop messages from NDSolve use Off as follows

Off[InterpolatingFunction::dmval];
sin2theta = 10^-13;
ms = 2*10^-5;
n = 25;
sfvetest[0] =
Interpolation[
Table[{T, -(45/(4*Pi^4*gs[T]))*(-Sqrt[90/(Pi^2*gs[T])]*(Mpl/
T^3)*(1 + (T/3)*Dgs[T]/gs[T]))*
NIntegrate[
epison^2*
Gama[epison,
T]*(theta2M[sin2theta, ms, epison, T, Lini]*
fve[epison, T, Lini] -
theta2Mbar[sin2theta, ms, epison, T, Lini]*
fvebar[epison, T, Lini]), {epison, 0, 50},
Method -> "LocalAdaptive"]}, {T, initialT, finalT, -0.01}]];
Do[Lvetest[i] =
NDSolveValue[{D[Lve[T], T] == sfvetest[i - 1][T],
Lve[initialT] == Lini}, Lve, {T, initialT, finalT}];
sfvetest[i] =
Interpolation[
Table[{T, -(45/(4*Pi^4*gs[T]))*(-Sqrt[90/(Pi^2*gs[T])]*(Mpl/
T^3)*(1 + (T/3)*Dgs[T]/gs[T]))*
NIntegrate[
epison^2*
Gama[epison,
T]*(theta2M[sin2theta, ms, epison, T, Lvetest[i][T]]*
fve[epison, T, Lvetest[i][T]] -
theta2Mbar[sin2theta, ms, epison, T, Lvetest[i][T]]*
fvebar[epison, T, Lvetest[i][T]]), {epison, 0, 50},
finalT, -0.01}]];, {i, 1, n}];


To accelerate code of about 53 times (on my computer) we can use Gauss quadrature instead of NIntegrate as follows

Needs["NumericalDifferentialEquationAnalysis"](*Run this first!*)


Code

episonmax = 20.;
Nbin = 500;
sfve[sin2theta_, ms_, T_, LveT_] := -
2]] . ((45/(4*Pi^4*gs[T]))*(-Sqrt[90/(Pi^2*gs[T])]*(Mpl/
T^3)*(1 + (T/3)*Dgs[T]/gs[T]))*#^2*
Gama[#, T]*(theta2M[sin2theta, ms, #, T, LveT]*
fve[#, T, LveT] -
theta2Mbar[sin2theta, ms, #, T, LveT]*fvebar[#, T, LveT]) & /@
ffve[sin2theta_, ms_] :=
Block[{Lve, solution1},
solution1 =
NDSolve[{D[Lve[T], T] == sfve[sin2theta, ms, T, Lve[T]],
Lve[initialT] == Lini}, Lve, {T, initialT, finalT},
Method -> {"EquationSimplification" -> "Residual"}];
Lve /. solution1[[1]]];
LveT = ffve[10^-13., 2*10^-5.]; // AbsoluteTiming


Visualization

LogLogPlot[LveT[T], {T, initialT, finalT},
PlotStyle -> {Dashed, Red}, PlotRange -> All, Frame -> True]


Update 1. At large Nbin>2000 it could be more effective trapezoidal rule implemented in a form

episonmax = 20.;
Nbin = 3000; h = episonmax/(Nbin - 1);
episoni = Range[0, 20, h];

trapIntegrate[data_List, h_] :=
h*(Plus @@ # - (First[#] + Last[#])/2) &[data]

sfve[sin2theta_, ms_, T_,
LveT_] := -trapIntegrate[((45/(4*Pi^4*gs[T]))*(-Sqrt[
90/(Pi^2*gs[T])]*(Mpl/T^3)*(1 + (T/3)*Dgs[T]/gs[T]))*#^2*
Gama[#, T]*(theta2M[sin2theta, ms, #, T, LveT]*
fve[#, T, LveT] -
theta2Mbar[sin2theta, ms, #, T, LveT]*fvebar[#, T, LveT]) & /@
episoni), h];
ffve[sin2theta_, ms_] :=
Block[{Lve, solution1},
solution1 =
NDSolve[{D[Lve[T], T] == sfve[sin2theta, ms, T, Lve[T]],
Lve[initialT] == Lini}, Lve, {T, initialT, finalT},
Method -> {"EquationSimplification" -> "Residual"}];
Lve /. solution1[[1]]];

LveT1 = ffve[10^-13., 2*10^-5.]; // AbsoluteTiming
`

It takes about 280 s on my laptop.

• Thank you very much for the comprehensive response. I've managed to resolve the previous issue, but I'm encountering difficulties with the complexity and singularities in using NIntegrate, which I lack experience in handling. The issue is evident as a small peak in my most recent figure. I attempted to use the Monte Carlo method for numerical integration, but unfortunately, it did not yield the desired results. Commented Apr 5 at 14:34
• @RoySmith In the new redaction your code is good. To remove messages see update to my answer. Commented Apr 5 at 20:00
• Thank you once again for your patient and helpful response. Your method has significantly aided me, although I've observed some instability when introducing a range of sin2theta and ms values. Therefore, I experimented with an alternative summation approach. Unfortunately, this method is quite slow, even after attempting to use ParallelSum. I'm beginning to wonder if this might be due to computational limitations within the Mathematica application or perhaps other reasons. Commented Apr 7 at 9:42
• @RoySmith Actually limitations are same as for C. You code could be optimize to reduce computation time from 146 s to 2 s. See update to my answer. Commented Apr 8 at 11:26
• Thank you very much for your attention. I have successfully obtained the approximate result using your summation method. Although there were instances when the process took nearly 5 hours for Nbin=2500, I consider this duration acceptable. Once again, thank you for your assistance; I am truly grateful. Commented Apr 10 at 7:21