Improve repeatability in NDSolve?

Questions
,
1

I am working with a system of differential equations in an algorithm which requires repeated solution using NDSolve. There are currently over 300 equations which are both differential and algebraic. The system is highly nonlinear and can be somewhat spikey in response. Most of the time the system solves in about 4 seconds. However, now and again it takes minutes rather than seconds. When this occurs stand-alone, and I see it taking so long, I can often abort the execution and start it again with shift-enter, and it often then completes in the usual 4 seconds, with exactly the same parameters.

I thought it may be taking different paths at different times based on numerical error, especially as it tries to adapt step sizes. So I gave it a MaxStepSize which is significantly smaller than the smallest expected features in the solution. That seems to have helped, but it has not completely alleviated the problem. I would appreciate any advice as to how I can steer the process toward a repeatable 4 seconds.

Kind regards,

David

2

Here is an example. The function ChainOutput[] uses NDSolve to produce its result. I executed it 4 times in a row.

  1. It is wrapped in a set of statements that show that it took 98 seconds cpu time within 355 second elapsed time to produce a solution. During execution Task Manager showed all 4 cores at 100%

  2. This time it solved in 4.6 seconds.

  3. This time 64 cpu seconds in 76 seconds elapsed time.

  4. This time 41 seconds.

I’ve looked to see – I’m not modifying initial conditions between runs. Unfortunately, I can’t share the notebook, but I hope someone can give me a clue as to how I can remove the random number generator from this, with hopefully a solution as quick as the fastest.

In[153]:= clock0 = AbsoluteTime[];
Print[Timing[chainOutput[100, 40, 70, Ga -> Ga]]];
timeExpended = AbsoluteTime[] - clock0

During evaluation of In[153]:= {97.828227,106.69089}

Out[155]= 354.7782637

In[156]:= Timing[chainOutput[100, 40, 70, Ga -> Ga]]

Out[156]= {4.5708293, 106.69089}

In[157]:= clock0 = AbsoluteTime[];
Print[Timing[chainOutput[100, 40, 70, Ga -> Ga]]];
timeExpended = AbsoluteTime[] - clock0

During evaluation of In[157]:= {64.475213,106.69089}

Out[159]= 75.5682788

In[160]:= Timing[chainOutput[100, 40, 70, Ga -> Ga]]

Out[160]= {41.262265, 106.69089}