Memory leak with NDSolve

Question

Bug solved in 11.2

---

Bug introduced in 10.3 or earlier and persists through 11.0

---

Edit The technical team confirms the memory leak in NDSolve and has forwarded an incident report to the developers.

---

This is not the first time the question in mentionned (cf in particular this recent question NDsolve memory leak), but there seems to be a serious issue of memory leak - or am I missing something?

First example (no leak) Let's loop 1000 times over the first example in the documentation of NDSolve:

$HistoryLength = 0;
s := NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, 
   y, {x, 0, 30}];
temp = MemoryInUse[];
Table[s, {i, 1000}];
MemoryInUse[] - temp

-1256

There is no significant increase in the memory usage.

Second example (memory leak?!) Now let's loop 1000 times with NDSolve using another example from the same documentation. I see no programming difference in terms of memory, however this time the used memory increases significantly (and for other systems of equations, it can increase much faster).

s2 := NDSolve[{y''[x] + Sin[y[x]] y[x] == 0, y[0] == 1, y'[0] == 0}, y, {x, 0, 30}]
temp = MemoryInUse[];
Table[s2, {i, 1000}];
MemoryInUse[] - temp

293 696

Additionally, Remove[s2] or Clear[s2] does not free any memory (I think that's normal, but what is the memory used for?).

---

How come does the memory usage increases, and how to overcome this serious problem?

Tested on MMA 10.3.0.0 with Ubuntu 15.10 64 bits

---

Update (09 Aug 2016): I just installed version 11.0 and the leak is still present.

Answer 1

Update 2017 11 10 The technical support report this bug has been solved in the 11.2 release.

---

Partial answer

The technical support came back to me with the following workaround: write the ODE(s) in the first order form. Some symbols are apparently leaking when MMA performs the reduction.

For example,

NDSolve[{y'[x] == z[x], z'[x] + Sin[y[x]] y[x] == 0, y[0] == 1, z[0] == 0}, y, {x, 0, 30}]

does not lead to memory leak.

However, the following DAE still leaks, while I see no necessary first-order reduction:

DAE = ({{x1'[t] == x3[t]}, {x2[t] (1 - x2[t]) == 
      0}, {x1[t] x2[t] + x3[t] (1 - x2[t]) == t}});
Clear[old, x, y]
ClearSystemCache[];
MemoryInUse[]
Do[old = MemoryInUse[];
 NDSolve[{DAE, x1'[0] == 1}, {x1, x2, x3}, {t, 0, 1}]; 
 Print[MemoryInUse[] - old];, {10}]
MemoryInUse[]

Present status is that I am waiting for their answer. I'll update the post accordingly.

Edit I was told by the support team that they do not know when this bug will be fixed (but they confirmed that it's a bug).

Answer 2

I have experienced serious memory leaks with NDSolve, NIntegrate, and FindRoot in every version of Mathematica I have owned (still on V9), perhaps over 20 years. They typically show up for me when the routine is called at a deep level in a complex program.

The only cure in many cases is to remove the offending routine and write your own. For NIntegrate, say a Simpson's rule, for FindRoot, say a Newton's method. If using NDSolve on an ODE, again a simple hand-coded solution. (I'm not saying that's the cure for this particular case, which is effectively cured by the HistoryLength trick. I'm saying that you will eventually run across one that seems 'unfixable' after trying $HistoryLength=0, Clearing[] every variable in sight, etc, etc.)

I know this sucks, but it's the only solution I have found in some cases. In my case, I am talking about a half-dozen cases over, say, several thousand programs. So, it's not the end of the world.