# How to speed up (or parallelize) NDSolve?

Is there a way to speed up the following code:

m = 1;
w = (a^2 + (m/p)^2)^(1/2);
z = (b^2 - a^2)/w^2;
y[z_] = ArcTan[z^(1/2)]/z^(1/2);
R[2, 4, 0, b_?NumericQ, a_?NumericQ] = (3 + 2 z - 3 (1 + z) y[z])/(z^2 w^3);
R[2, 2, 1, b_?NumericQ, a_?NumericQ] = (-3 + (3 + z) y[z])/(z^2 w^3);

II[n_, q_, r_, b_?NumericQ, a_?NumericQ] := ((b^(2 q + 2)) (a^(r + 1)) )/(4 \[Pi]^2 (2 q)!!) NIntegrate[p^(n + 1) R[n, q, r, b, a] E^-(p^2 + m^2)^(1/2), {p, 0, \[Infinity]}];

t0 = 0.1; tend = 1;
AbsoluteTiming[
s1 = NDSolve[{D[II[2, 4, 0, b[t], a[t]], t] == a[t], D[II[2, 2, 1, b[t], a[t]], t] == b[t], b[t0] == 1, a[t0] == .1},
{b[t], a[t]}, {t, t0, tend}]
]


Output:

Is it possible to parallelize NDSolve since that's where most of the time is spent?

This is a demo code. The actual code takes much longer (~30 mins).

• The integration NIntegrate[p^(n + 1) R[n, q, r, b, a] E^-(p^2 + m^2)^(1/2), {p, 0, \[Infinity]}] might be simplified because R[n, q, r, b, a]  doesn't depend on p! For even n analytical integration is possible! May 11, 2021 at 12:04
• @UlrichNeumann R[n, q, r, b, a] depends on p through w. May 11, 2021 at 12:09
• This video is possibly relevant to you. May 11, 2021 at 12:21
• @SjoerdSmit Thanks! It seems parallelizing NDSolve is a bit of a hassle. They don't provide a working example. Is there an easier way to speed up this code? May 11, 2021 at 17:12
• @SunilJaiswal Yes, unfortunatlely solving differential equations is not easy to parallellise. That's just in the nature of the problem. If I were you, I'd focus on trying to improve the evaluation speed of the NIntegrate in your code by finding good option settings for it. I assume that there's no symbolic version you can find for this integral? May 11, 2021 at 17:35

We can reduce computational time by 25 times with extended expressions for coefficients of equations and by reducing integration time with option PrecisionGoal -> 4 as follows

Clear["Global*"]

m = 1;
w = (a^2 + (m/p)^2)^(1/2);
z = (b^2 - a^2)/w^2;
y1 = ArcTan[z^(1/2)]/z^(1/2);
(*R[2,4,0,b_,a_]:=*)R1 = (3 + 2 z - 3 (1 + z) y1)/(z^2 w^3);
(*R[2,2,1,b_,a_]:=*)R2 = (-3 + (3 + z) y1)/(z^2 w^3);
(*II[n_,q_,r_,b_?NumericQ,a_?NumericQ]*)

A1 = (((b^(2 q + 2)) (a^(r + 1)))/(4 \[Pi]^2 (2 q)!!) p^(n +
1) R1 E^-(p^2 + m^2)^(1/2)) /. {n -> 2, q -> 4, r -> 0}; A11 =
D[A1, a]; A12 = D[A1, b];

a11[x_?NumericQ, y_?NumericQ] :=
NIntegrate[
If[x != y, A11 /. {a -> x, b -> y}, (
E^-Sqrt[1 + p^2] p^5 x^10 (7 - 8 p^2 x^2))/(
26880 Sqrt[1/p^2 + x^2] (\[Pi] + p^2 \[Pi] x^2)^2)], {p, 10^-8,
40}, PrecisionGoal -> 4];

a12[x_?NumericQ, y_?NumericQ] :=
NIntegrate[
If[x != y, A12 /. {a -> x, b -> y}, (
E^-Sqrt[1 + p^2] p^7 x^10 Sqrt[1/p^2 + x^2] (35 + 32 p^2 x^2))/(
13440 \[Pi]^2 (1 + p^2 x^2)^3)], {p, 10^-8, 40},
PrecisionGoal -> 4];

A2 = (((b^(2 q + 2)) (a^(r + 1)))/(4 \[Pi]^2 (2 q)!!) p^(n +
1) R2 E^-(p^2 + m^2)^(1/2)) /. {n -> 2, q -> 2, r -> 1}; A21 =
D[A2, a]; A22 = D[A2, b];

a21[x_?NumericQ, y_?NumericQ] :=
NIntegrate[
If[x != y, A21 /. {a -> x, b -> y}, (
E^-Sqrt[1 + p^2] p^5 x^7 (14 + 5 p^2 x^2))/(
840 Sqrt[1/p^2 + x^2] (\[Pi] + p^2 \[Pi] x^2)^2)], {p, 10^-8, 40},
PrecisionGoal -> 4];
a22[x_?NumericQ, y_?NumericQ] :=
NIntegrate[
If[x != y, A22 /. {a -> x, b -> y}, (
E^-Sqrt[1 + p^2] p^7 x^7 Sqrt[1/p^2 + x^2] (7 + 5 p^2 x^2))/(
140 \[Pi]^2 (1 + p^2 x^2)^3)], {p, 10^-8, 40}, PrecisionGoal -> 4];

eqs = {a11[x[t], y[t]] x'[t] + a12[x[t], y[t]] y'[t] - x[t] == 0,
a21[x[t], y[t]] x'[t] + a22[x[t], y[t]] y'[t] - y[t] == 0};
t0 = 0.1; tend = 1;
AbsoluteTiming[
sol = NDSolve[{eqs, y[t0] == 1, x[t0] == .1}, {x[t], y[t]}, {t, t0,
tend}]]


Note, that we need exact expressions for a11,a12,a21,a22 as a limit at a->b. Computation time is

(*{14.3547, {{x[t] ->...}}}*)


Visualization

Plot[Evaluate[{x[t], y[t]} /. sol[[1]]], {t, t0, tend},
AxesLabel -> Automatic, PlotLegends -> {"a", "b"}]


• Thanks! Although the code becomes specific for this particular case, this might work for my problem. May 12, 2021 at 7:41
• @SunilJaiswal We can also compile some functions with Parallelize option before NDSolve. But NDSolve runs quickly for this kind of system with known coefficients a11, a12, a21, a22`. Time mostly spend for computation coefficients. It could be better to improve code for your original system, not for toy example. May 12, 2021 at 11:50