# Use of ParallelDo with an iterative or recursive calculation

Is it possible to use ParallelDo (or any similar parallelization method) to speed up the calculation of an iterative calculation, which takes a previous calculation result as input to new calculation? It seems to me not possible conceptually, since each round will have to wait for the first to get its input, but I wonder if there are ways around this?

For example, if f1[r] is some predefined arbitrary continuous input function and f2[r] is some output function (over the r variable that runs from 0 to some number a, with some step called step, e.g. say step=a/100), then a calculation like this example

f1=1;
Do[
Do[
f2[r2]=NIntegrate[r1 f1[r1] BesselJ[0, r1 r2] E^(-I (r1^2 + r2^2)/2),{r1,0,a}],
(r2,0,a,step)];
f1=Interpolation[Table[{r2,f2[r2]},{r2,0,a,step}]],
{kk,100}]


would be difficult to speed up by parallelization (e.g. using ParallelDo), right?

Are there any suggestions to speed this up by parallel methods or other methods that are not to too difficult for a beginner in Mathematica to understand?

The only thing that can be parallelized here is the generation of the interpolation points. For example:

Clear[integrate, interpolate];
integrate[f_, r2_?NumericQ, a_?NumericQ] := NIntegrate[
r1 f[r1] BesselJ[0, r1 r2] E^(-I (r1^2 + r2^2)/2),
{r1, 0, a}
];
interpolate[f_, a_?NumericQ, step_?NumericQ] := Interpolation[
ParallelTable[
{r2, integrate[f, r2, a]},
{r2, 0, a, step}
]
];


Iterate the procedure twice, starting with a constant function:

results = NestList[
interpolate[#, 1, 0.01] &,
1 &,
2
];


Plot the obtained functions at each iteration:

ReImPlot[Evaluate[Through[results[r2]]],
{r2, 0, 1},
PlotLegends -> Range[0, Length[results] - 1]
]


# Edit

I don't think that with NIntegrate you can speed this up any further. The free parameter r2 appears in the body of the integral, so I don't think that interpolation is going to help you too much here. You could make an interpolation over both r1 and r2 of which you then take the antiderivative over r1, but you also risk introducing numerical error that way. There's probably a reason why NIntegrate isn't faster here.

I do have another idea for speeding this up, though. We can use ParametricNDSolveValue to do the integration instead by re-phrasing the problem and considering r2 to be a parameter. This should reduce the top-level overhead involved in calling NIntegrate multiple times. Here's how to do it:

Clear[integrate, interpolate];
integrate[f_, a_?NumericQ] := ParametricNDSolveValue[
{y'[r1] == r1 f[r1] BesselJ[0, r1 r2] E^(-I (r1^2 + r2^2)/2), y[0] == 0},
y[a],
{r1, 0, a},
{{r2, 0, a}}
];
interpolate[f_, a_?NumericQ, step_?NumericQ] := With[{
int = integrate[f, a]
},
Interpolation[
ParallelTable[
{r2, int[r2]},
{r2, 0, a, step},
Method -> "CoarsestGrained"
]
]
];
results = NestList[interpolate[#, 1, 0.01]&, 1&, 20];

ReImPlot[
Evaluate[Through[results[r2]]],
{r2, 0, 1},
PlotLegends -> Range[0, Length[results] - 1], PlotRange -> All
]


This should be significantly faster than the method based on NIntegrate. You can also take a look at FunctionInterpolation to do the Interpolation step, since this might save you from having to sub-sample the integral more than you really need. For example:

int = Quiet @ FunctionInterpolation[
integrate[1&, 1][r2],
{r2, 0, 1},
InterpolationPoints -> 50
]
ReImPlot[int[r2], {r2, 0, 1}]


FunctionInterpolation cannot be parallelized, though, so it's a trade-off

• I appreciate your advice. My problem is that I find the loop I mentioned to be quite slow (especially as the functions evolve from iteration to iteration), so I suspected the culprit is the numerical integration, rather than the interpolation. Do you agree? (would your parallelizing of the interpolation points make much difference here?) How would you speed up such loop? Would you use the anti-derivative techniques you mentioned in your other post instead of integrating? (mathematica.stackexchange.com/questions/225127/…) – user135626 Nov 17 '20 at 22:42
• ... or will the primitive anti-derivative of the interpolation not help much in this case because the integration has a kernel made up of other functions (Bessel and r1)? It would be nice if a method like that using anti-derivative of interpolation (without integration) could be used here. I think it would speed it up greatly when combines with your answer here, which basically paralellizes the loops over the variable r2 but not the integration step itself. Any advice on this? – user135626 Nov 17 '20 at 23:36
• Another question: isn't your answer fundamentally the same as parallelizing the inner Do loop in my original question (i.e. the Do loop over r2 variable)? Or is your coding method faster? – user135626 Nov 18 '20 at 0:13
• @user135626 I updated the answer. I think you'll find it useful, but let me know. As for your last question: yes, it's sorta similar, but parallelizing assignments to variables is not the best way to do things because it has a lot of communication overhead. That's why ParallelTable is better here than ParallelDo. – Sjoerd Smit Nov 18 '20 at 14:13
• This is a great answer. Really appreciate it! I just love the clarity of your treatment. May I finally ask: (1) I would like to be good in Mathematica and have the opportunity to enroll in any courses to master it more: do you recommend any particular courses? (2) Am I correct in that you used the With command in your answer to solve the paramteric ODE before the r2 parameter values are assigned in ParallelTable as that would cause error in ParametricNDSolveValue? (3) But isn't ParallelDo also assigning different r2 values to different kernels, just as ParallelDo would do here? – user135626 Nov 19 '20 at 4:15