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

Question

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!

Answer 1

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... =)