I'm seeking to find solutions to a numerical integration with a large set of parameter combinations (basically, I'm doing a brute parameter sampling). Yet, I believe the memory of the computer is storing unnecessary information, but I can't understand why: the memory stores near 16Gb of memory coming from Mathematica Kernel, while the output is only ~25mb. This sometimes makes the computer run out of memory and Kernel quitting. I am also afraid this may slow my computations.

I will abbreviate a lot the code because I just want to know if I'm doing something fundamentally wrong with my rationale, specially on point 3 below. Basically, I am numerically integrating some equations for many parameter values. The output of the integration is then used to compute values of the dependent variables at specific time points. These values are my output ~25 mb.

Note that I simplified my code a lot, and also left some of it blank because it is unnecessary to explain what I am doing. What I'm doing is the following:

1- Considering a very large table

Table4 =
     {10^b, 10^g, 10^v, 10^k}, {b, -5, 5, 0.5}, {g, -7, 3, 
      0.5}, {v, -11, -1, 0.5}, {k, -8, 2, 0.5}

2- I have a function deqns4 that is simply a system of ODEs to be numerically integrated and is called several times in point 3 (This is a random system of equations, my real one has 13 dependent variables).

deqns4[b_, g_, v_, k_] :=
  {A'[t] == ki H[t] + k X[t] - b A[t],
   X'[t] == g A[t] - v X[t] - kc X[t],
   A[0] == 2,
   X[0] == 1

3- Then, the following function calls deqns4 several times, to **numerically integrate the equations and output the values of the dependent variables at some time points, for the given parameter values in Table4:

lsim4[{b_, g_, v_, k_}, ba_] :=  
Block[{ep, initconds, sssol0, tsim4, tsim3, tsimc50, tsimc100, c1,     c2, c3, c4, c5, c6, c7},
     (*This finds initial conditions*)   
sssol0 = Cases[
           ki H + k X - b A == 0,
           g A - v X - kc X == 0} /. ki -> 2 /. kc -> 0.1
      {A, X}
     {Rule[_, _?NonNegative] ..}
     ] /. {A -> A0, X -> X0};

(*This and the following commands use the initial conditions above to numerically integrate the system of equations "deqns4"*)
tsim4 = Flatten[
          {H[t] == 0.009 Exp[-0.1 t] + ba },
          deqns4[b, g, v, k]
          ] /. kc -> 0.1 /. sssol0],
     {A, X, H},
     {t, 0, 9720}, MaxSteps -> 1*^12

(*tsim3, tsimc50, and tsimc100 are also numerical integration of deqns4,    similarly to tsim4*)
tsim3 = Flatten[NDSolve[(*...left blank to simplify...*)], 1];
tsimc50 = Flatten[NDSolve[(*...left blank to simplify...*)], 1];
tsimc100 = Flatten[NDSolve[(*...left blank to simplify...*)], 1];

(*then I want the output to be a set of values computed from the \ numerical integration above*)   c1 = (A[324] + X[324]) 100 /. tsim4;
     c2 = (A[9720]) 100/(A[0]) /. tsim4;
     c3 = (X[3240]) 100 /. tsim4;
     c4 = (X[9720]) 100 /. tsim4;
     c5 = (A[3240] + X[3240]) 100 /. tsim3;
     c6 = (X[3240]) 100 /. tsimc100;
     c7 = (X[3240]) 100 /. tsimc50;
     {c1, c2, c3, c4, c5, c6, c7}

lsimsols = ParallelMap[lsim4[#, 1.2] &, Table4];

(*output is a list with the form lsimsols={{0.1,0.99,100,0.52,76,100,100},{0.01,0.99,90,53,75,100},etc}*)

My questions are the following:

  1. Is NDSolvestoring information each time it runs, though it is within a Block function?
  2. The function lsim4 is only outputing c1 to c7 for each element of Table4. This is what is being stored by lsimsols, and is the only information I am interested in. Yet, a lot more data is being stored by the kernel. Is there any way for me to avoid storage of so much information?
  3. Even without specifying MaxStepSize, takes me ~4 hours to run all the values of parameters. Is there any more efficient way of doing this parameter sampling other than ParallelMap?

Thanks in advance!


  • as suggested in the comments, I set $HistoryLength=0 which did not change all the memory allocated to Mathematica - I still get about 16 Gbs of stored information.
  • Also, ParametricNDSolve seems to produce faster results than NDSolve.

Edit #2:

  • Someone suggested to use Map instead of ParallelMap, as this would possibly allow less memory usage. This didn't make any improvements.
  • The only way I was able to attain a partial solution to the problem above is by splitting Table4 in multiple parts and run deqns4 for each part separately and then joining everything.
  • 2
    $\begingroup$ Have you tried setting $HistoryLength to something other than Infinity (the default)? $\endgroup$
    – bill s
    Commented May 2, 2013 at 14:47
  • $\begingroup$ I have had that problem twice. The kernal stops when you run out of swap file space, so you might get it done by increasing the swap file. If this is something to be done often, you will need a better algorithm. $\endgroup$ Commented May 2, 2013 at 15:06
  • $\begingroup$ @Sosi, you could use ParametricNDSolve. As pointed out above try $HistoryLength=0 $\endgroup$
    – user21
    Commented May 2, 2013 at 15:23
  • $\begingroup$ @bills - I've never heard of that. I'll look into it! @FredKline - sorry to ask, but how do you do that? @ruebenko - I'll check ParametricNDSolve. I think I've used it before. $\endgroup$
    – Sos
    Commented May 2, 2013 at 15:47
  • $\begingroup$ Ok, I set $HistoryLength=0 but Mathematica still steals all memory. Using ParametricNDSolve made the computations faster though, but I then have to export the output, quit the kernel, import, and only then I can actually do something in mathematica $\endgroup$
    – Sos
    Commented May 3, 2013 at 15:02


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.