# Speed up for NDSolveValue containing a coefficient of interpolate function

Could you give me some advice to speed up this code? I think especially NDSolveValue takes a lot of time.

For j=1, it the NDSolveValue seems to be fast. But for j=2, I use the result of previous NDSolveValue(interpolate function) as a coefficient for another NDSolveValue, which seems takes lots of time.

I am running on a server of 24 cores, but I am not sure how I can take the advantage of that to speed up.

u = 2;
c = 3000;
it = 3;
Array[h, {u, c}];
Do[Do[h[i, k] = 0, {i, u}], {k, c}];
F[x_] := Sin[x];
AA[x_] := 1/(x^2 + 1);
w[x_] := F[x] + AA[x];
v[x_] := Sum[w[x + i*2 Pi], {i, -it, it}];

Do[n = 1;
If[j == 1,
While[n <= c,
h[j, n] =
NDSolveValue[{y''[x] == Cos[y[x]], y'[0] == 0, y[0] == n},
y, {x, 0, 1},
"Method" -> {"EquationSimplification" -> {Automatic,
"TimeConstraint" -> Infinity}}, AccuracyGoal -> 3]; n++];,
While[n <= c, h[j, n] = NDSolveValue[{y''[x] == Cos[y[x]]
+ Sum[v[h[j - 1, n][x] - h[j - 1, i][x]], {i, 1, c}],
y'[0] == 0, y[0] == n}, y, {x, 0, 1},
"Method" -> {"EquationSimplification" -> {Automatic,
"TimeConstraint" -> Infinity}}, AccuracyGoal -> 3.];
n++];];, {j, u}];


I am not sure how I can accommodate parallelize for the NDSolveValue. I tried to use Compile instead of functions, but it doesn't seem to be compatible with interpolate function.

• Questions with so much code tend not to attract much attention. Can you capture the essence of your problem with a shorter example? – bbgodfrey Oct 20 '17 at 0:01
• I revised the previous code to a simpler one to address the essence only. – user16308 Oct 20 '17 at 3:15

First setting a more reasonable amount of DEs to solve for my laptop,

c = 50;


I find the OP's code takes 4.45876 seconds to run.

The dependence of h[j, n] on previous solutions limits what can be parallelized. For each j, the following parallelizes the solutions for all n. Since you don't want the subkernels calling back to the main kernel to set the value of h[j, n], I changed the loop to ParallelTable and stored the list of solutions in hh[j]. The solution h[j, n] will then correspond to hh[j][[n]]. (I don't know if the Method setting is necessary for high values of c; it was unnecessary for the following.)

It's 2.5 times faster for c = 50 on my 4/8-core i7 laptop. For each j, the new definition of hh have to be distributed to each subkernel, which probably accounts for some of the limitations in gain of speed.

Clear[hh];
Do[
If[j == 1,
hh[j] = ParallelTable[
NDSolveValue[{y''[x] == Cos[y[x]], y'[0] == 0, y[0] == n},
y, {x, 0, 1}, AccuracyGoal -> 3],
{n, c}],
hh[j] = ParallelTable[
NDSolveValue[{y''[x] ==
Cos[y[x]] + Sum[v[hh[j - 1][[n]][x] - hh[j - 1][[i]][x]], {i, 1, c}],
y'[0] == 0, y[0] == n}, y, {x, 0, 1}, AccuracyGoal -> 3.],
{n, c}]
];,
{j, u}]; // AbsoluteTiming
(*  {1.71541, Null}  *)


Spot check a value with OP's code:

h[1, 23][0.62]
hh[1][[23]][0.62]
%% - %
(*
22.8948
22.8948
0.
*)


Check another solution over the domain (the difference is zero at all sampling points):

Plot[h[2, 23][x] - hh[2][[23]][x], {x, 0, 1}]