Slow evaluation of Recurrence plot data from NDSolve: Performance Tuning

Question

I am trying to use a recurrence plot to pinpoint the location when a system re-visits a previous point in its phase portrait. The system (a thin liquid film) is governed by a non-linear differential equation and is quite easily solved using NDSolve as shown:

{xMin, xMax} = {-(π/0.0677), π/0.0677}; 
k = 0.0677; 
TMax = 1650; 
m = 5; 
S = 100; 
Ga = -3^(-1); 
Quiet[uSolpbc = 
     u /. NDSolve[{D[u[t, x], 
        t] == (-S)*D[u[t, x]^3*D[u[t, x], x, x, x], x] + 
                 Ga*D[u[t, x]^3*D[u[t, x], x], x] - 
        m*D[(u[t, x]/(1 + u[t, x]))^2*
                       D[u[t, x], x], x], 
      u[0, x] == 1 - 0.01*Cos[k*x], u[t, xMin] == u[t, xMax],                  
      Derivative[0, 1]*u[t, xMin] == Derivative[0, 1]*u[t, xMax],                  
      Derivative[0, 2]*u[t, xMin] == Derivative[0, 2]*u[t, xMax],                  
      Derivative[0, 3]*u[t, xMin] == Derivative[0, 3]*u[t, xMax]}, 
     u, {t, 0, TMax}, {x, xMin, xMax}, MaxSteps -> 100000, 
     Method -> {"MethodOfLines", "Method" -> "LSODA", "TemporalVariable" -> t,"SpatialDiscretization" ->                      {"TensorProductGrid", "MinPoints" -> 800, "MaxPoints" -> 1200, 
                  "DifferenceOrder" -> 5}}][[1]]]

The film profile is plotted as:

Plot[uSolpbc[0.99 TMax, x], {x, xMin, xMax}]

I am coding the Recurrence Plot in a block:

Block[{stepSize = 2, end = TMax, tt, ττ, rd},
 rd = ParallelTable[
   UnitStep[
    0.01 - Norm[uSolpbc[tt, 0] - uSolpbc[ττ, 0], 2]], {tt, 
    0, end, stepSize}, {ττ, 0, end, stepSize}];
 MatrixPlot[rd]]

The grainy nature of the plot can be improved by choosing a smaller time step; say, stepSize = 1/20 or less. However, the current stepSize of 2 itself is quite slow (Timing reveals 0.344 seconds; AbsoluteTiming reveals 3.33 seconds).

I have a feeling that it is the nested nature of these functions (ParallelTable, UnitStep and Norm) that is taking a long time. Is there some way I can improve this? Currently, choosing a much smaller time step leads to an "out of memory" error (I have ~ 8 GB of memory that is exhausted!).

It is not Nest or Map that I am looking for, and if it is, I guess I am lost on its application in this case.

Answer 1

One can use the new function DistanceMatrix[] for the purpose; this avoids repeated computations (since the underlying matrix is symmetric).

With[{stepSize = 2, end = TMax},
     MatrixPlot[UnitStep[0.01 - DistanceMatrix[uSolpbc[Range[0, end, stepSize], 0]]]]]

and your plot is produced very quickly, without the need to invoke parallelization.

Answer 2

The plot in the question must have been obtained with stepSize = 15, not stepSize = 2. Using the latter value gives a smooth plot,

The computation takes about 78 sec on my PC. To address the specific issue in the question, the run time can be reduced by two orders of magnitude using

Block[{stepSize = 2, end = TMax, tt, rd}, 
    tSolpbc = Table[uSolpbc[tt, 0], {tt, 0, end, stepSize}]; 
    rd = ParallelTable[
        UnitStep[0.01 - Norm[tSolpbc[[nt]] - tSolpbc[[nτ]], 2]], 
        {nt, end/stepSize}, {nτ, end/stepSize}]; MatrixPlot[rd]]

which produces the same plot. Evidently, most of the run time used in the original computation was consumed by computing uSolpbc (end/stepSize)^2 times. The revised computation computes it only end/stepSize times.