# Generating and analyzing an Interpolating function from many interpolating functions: speed up

Is there a clever way to speed up this code?

    v0 = 2 10^-5;
ppi = E^(-(p^2/4)) (2/\[Pi])^(1/4);
h = 1/60; (*stepsize*)
a = -10; (*discretization range start*)
b = 10;  (*discrtization range end*)
ic = Table[f[p, 0] == ppi, {p, a, b, h}];
state = Table[f[p, t], {p, a, b, h}];
eq = Table[{I D[f[p, t], t] ==
1/4 p^2 f[p, t] -
I v0/(2 h^3) (f[p + 2 h, t] - f[p - 2 h, t] -
2 (f[p + h, t] - f[p - h, t]))}, {p, a, b,
h}] /. {f[b + h, t] -> 0, f[b + 2 h, t] -> 0, f[a - h, t] -> 0,
f[a - 2 h, t] -> 0};
sol = First@NDSolve[{eq, ic}, state, {t, 0, 2 \[Pi] 10}];
pplist[t1_] :=
pplist[t1] =
Table[{i ,
Abs[(state /. sol) /. t -> (t1*2 \[Pi])][[-a/h + i/h +
1]]^2}, {i, a, b, h}];
plot = Interpolation[pplist[10]]
Plot[plot[x], {x, -3, 2}, PlotRange -> All]


I think the biggest bottleneck is the creation of the Table from the list of interpolating functions generated by NDSolve.

plot = Interpolation[pplist[10]]

This takes like 70 seconds on my Laptop.

A tremendous speedup would be amazing because I need to run this code probably a lot of times. In end the end I have to analyze (and maximize) the distance between the maximums of the resulting function as a function of $$t$$ for a lot of different input parameters. Thank you a lot for any help and hint and whatever :)

Edit: Not sure if I should post that as an answer to my own question. However, I found out that I get a good speedup by changing this line:

pplist[t1_] :=
pplist[t1] =
Table[{i ,
Abs[(state /. sol) /. t -> (t1*2 \[Pi])][[-a/h + i/h +
1]]^2}, {i, a, b, h}];


into this line:

pplist[t1_] :=
pplist[t1] =
Table[{i ,
Abs[(state[[-a/h + i/h +
1]] /. sol) /. t -> (t1*2 \[Pi])]^2}, {i, a, b, h}];


Now the bottleneck seems to be NDSolve (which takes 13 seconds on my laptop). Of course any ideas for further speed up still very appreciated!

• Should have the tag performance-tuning, but which to remove? I'd get rid of list-manipulation or compile. It's not really about either, imo. – Michael E2 Mar 16 '20 at 13:23
• How long does the NDSolve take? – Michael E2 Mar 16 '20 at 13:33
• thanks for your comment. I got rid of list-manipulation and compile and added performance-tuning. NDSolve takes 13 seconds on my laptop, which I also added to the question. – Luke Mar 16 '20 at 13:38
• Are you aware that Mathematica is able to solve the underlying pde directly ? – Ulrich Neumann Mar 16 '20 at 13:56
• I tried that, but the results I got this way do not make any sense. I guess because NDSolve has no "method of lines" implemented for third order pde's. Please correct me if I am wrong. – Luke Mar 16 '20 at 14:05

As far as I have checked, the following should produce the same result within a couple of seconds. The most important point towards performance is to use NDSolveValue because that avoids all sorts of symbolic computations and replacements. Using a sparse matrix to set up the system is just convenient (if one has a bit of experience with that).

v0 = 2 10^-5;
h = 1/60;
a = -10;
b = 10;
plist = Range[N@a, b, h];
A = With[{c = N[I v0/(2 h^3)]},
SparseArray[{
Band[{1, 1}] -> 1./4 plist^2 + 0. I,
Band[{1, 3}] -> -1 c,
Band[{1, 2}] -> +2 c,
Band[{2, 1}] -> -2 c,
Band[{3, 1}] -> +1 c
},
{1, 1} Length[plist], 0. + 0. I
]
];
Y = NDSolveValue[{
I X'[t] == A.X[t],
X[0] == E^(-(plist^2/4)) (2/π)^(1/4)
}, X, {t, 0, 2 π 10}];
plot = Interpolation[Transpose[{plist, Abs[Y[2. Pi 10]]^2}]];
Plot[plot[x], {x, -3, 2}, PlotRange -> All]


The interpolation part is done within a 1.5 milliseconds which is quite a speed-up compared to the over 70 seconds from OP's code. Moreover, the ODE is solved approximately twice as fast this way.

By the way, if you require the solution only at few specific times, you should get it faster and more precisely with

MatrixExp[A/I (2. Pi 10), E^(-(plist^2/4)) (2/π)^(1/4)]


as you have probably learnt in your quantum mechanics class... =)

• Thank you so much! This is extremely helpful :) Do you mean that if I want to produce a 3DPlot of the solution with p and t as x and y axis I should use the NDSolveValue and if I already know which time I am interested in the matrix exponential? – Luke Mar 16 '20 at 14:18
• Yes. Exactly. =) – Henrik Schumacher Mar 16 '20 at 14:23
• Do you know of a way to further speed up the Matrix Exponential calculation by using packed arrays or other tricks? – Luke Apr 23 '20 at 10:43
• Nope. Matrix exponentials a nice to write down but very slow to compute numericallly. They can be used to compute the solution at a few given times. But if you want to compute the full time evolution, a better strategy would be to solve the PDE with some implicit time integration scheme (e.g., backward Euler or Crank-Nicolson). – Henrik Schumacher Apr 23 '20 at 11:41
• Yes, posting it in a separate post might be meaningful. But in this case, please describe first what you actually want to do. This is almost certainly an X-Y-question: You ask how to do something that does not solve your problem or at least not very efficiently. – Henrik Schumacher Apr 23 '20 at 12:20