# 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