kernel Local crashes during the course of a For loop

Question

I want to find the minimum of the function r[x,p] subject to the constraint w[x,p,t]==0 (t is a fixed parameter). To find the minimum I use the function NMinimize.

The constraint w[x,p,t]==0 is a rather complicated numerical function. It depends on the function evolve[x,p,t] which is the numerical solution of some ordinary differential equation.

I use a For loop to find the minima for different values of the constraint parameter t.

My problem is the following: At certain instants of time, the kernel Local quits without errors or warnings before the loop is finished. The error message is "The kernel Local has quit (exited) during the course of an evaluation". The time when the kernel quits is always different, sometimes it happens very early, sometimes rather late. When I evaluate NMinimize for the parameter t where the kernel crashed (and also for parameters beyond that point), it works well without any error.

First I thought the reason for the crash could be a lack of memory, but if I monitor the memory during the calculation, everything looks fine.

Here is the code of a simplified version of my problem.

(*definition of the constraint*)
h[x_, p_] := -x^2 + p^2 + x^4;

evolve[ xini_, pini_, time_] := 
 Module[{x, p, t}, 
  Hold[{x[time], p[time]} /. 
   NDSolve[{x'[t] == p[t], p'[t] == -x[t], x[0] == xini, p[0] == pini}, {x, p}, 
           {t, 0, time}, MaxSteps -> Infinity][[1]]
 ]];

w[x_?NumericQ, p_?NumericQ, t_?NumericQ] := h @@ ReleaseHold[evolve[x, p, t]] - h[x, p];

(*definition of the function which should be minimized*)
r[x_, p_] := {x - 1, p}.{{0.1, 0}, {0, 0.1}}.{x - 1, p};

fct = First[{Print[#]; 
 NMinimize[{r[x, p], w[x, p, #] == 0}, {x, p}, 
  Method -> "NelderMead", MaxIterations -> 10^3]}] &;

(*loop over different parameters t*)
dt := 0.05;
tmin := 0;
tmax := 40;
maxt := Floor[(tmax - tmin)/dt];
times = Table[tmin + dt*n, {n, 0, maxt}];

positions = Table[0, {Length[times]}];
values = Table[0, {Length[times]}];

For[n = 1, n <= Length[times], n++,
  tmp = fct@times[[n]];
  positions[[n]] = {x, p} /. (tmp[[2]]);
  values[[n]] = tmp[[1]];
];

When I replace the loop over t by

fct/@times;

I get the same error.

I am using Mathematica 9.0.1.0

Answer 1

The random crashing of the kernel seems to arise due to your definition of evolve and fct. I'll keep

h[x_, p_] := -x^2 + p^2 + x^4;

and change the definition of evolve to

evolve[xini_, pini_, time_] := 
 Module[{x, p}, 
  Hold[{x[time], p[time]} /. 
   First@NDSolve[{x'[t] == p[t], p'[t] == -x[t], x[0] == xini, p[0] == pini}, 
                 {x, p}, {t, 0, time}, MaxSteps -> Infinity]
  ]];

as the t you put into Module might occasionally collide with the t inside NDSolve, although it shouldn't. The following definitions are kept unchanged

w[x_?NumericQ, p_?NumericQ, time_?NumericQ] := h @@ ReleaseHold[evolve[x, p, time]] - h[x, p];
r[x_, p_] := {x - 1, p}.{{0.1, 0}, {0, 0.1}}.{x - 1, p};

and Set(=) instead of SetDelayed (:=) is used for

dt = 0.05;
tmin = 0;
tmax = 40;
maxt = Floor[(tmax - tmin)/dt];
times = Table[tmin + dt*n, {n, 0, maxt}];

just because it is more appropriate.
Your definition of fct looks a little strange to me. But I guess the purposse of this construction is to monitor the progress of the calculation, while having the result as the output, and you simply misused the curly braces for lists, where you should have used parentheses for grouping (see tutorial):

fct = (Print[#];
       NMinimize[{r[x, p], w[x, p, #] == 0}, {x, p}, 
                  Method -> "NelderMead", MaxIterations -> 10^3]) &;

Unfortunately the kernel still crashes occasionally when fctis mapped over times (fct/@times) and it is not clear what causes this buggy behavior. However, an easy way to deal with it and at the same time speeding up the calculation significantly (more than 3.5 times on my 4 core PC) is to use ParallelMap:

results = ParallelMap[fct, times, Method -> "FinestGrained"];

The values and positions can then be extracted using

values = results[[All, 1]];
positions = {x, p} /. results[[All, 2]];

and plotted with

ListPlot[Transpose[{times, values}], PlotRange -> All, Frame -> True, FrameLabel -> {"time", "value"}]

and

ListPlot[positions, Frame -> True, FrameLabel -> {"x", "p"}]

or

ListPlot[{Transpose[{times, First@Transpose[positions]}], 
          Transpose[{times, Last@Transpose[positions]}]}, 
          PlotStyle -> {Black, Blue}, PlotRange -> All, Frame -> True, 
          FrameLabel -> {"time"}, PlotLegends -> {"x", "p"}]

---

But why is the kernel crashing using Map and how can ParallelMap still successfully perform the calculations?

As the messages

LinkObject::linkd: Unable to communicate with closed link LinkObject["C:\Program Files\Wolfram Research\Mathematica\10.0\MathKernel" -subkernel -noinit -mathlink -noicon,331,6]. >>
KernelObject::rdead: Subkernel connected through KernelObject[3,local] appears dead. >>

reveal, the Subkernels still crash at random frequency due to a lost connection.
But fortunately, the Masterkernel keeps running and is requeueing the failed calculation properly and therefore enables a successfull completion of the calculations:

Parallel`Developer`QueueRun::req: Requeueing evaluations {11} assigned to KernelObject[3,local,<defunct>].
LaunchKernels::clone: Kernel KernelObject[3,local,<defunct>] resurrected as KernelObject[5,local]. >>