Short version: Do the methods (or other things the function uses to solve the equations like time steps or error estimates) that NDSolveValue uses by default depend on the types and quantity of equations passed to it? If so, is there an easy way to see what is automatically being used in different cases? I really just need to be more aware of the ways I can look under the hood at how Mathematica is solving my equations. Below is more background and specifics if needed.

EDIT The answer by MichaelE2 below gives good things I can edit in NDSolve, but what if I want to look at what NDSolve is doing automatically in each case? How can I see what methods, time steps, errors, etc. are being used at specific points in the time integration?

I am attempting to solve multiple coupled, first order, ordinary differential equations using NDSolveValue. Here is the background. I am modeling some processes that take place in a biological cell. I have a set of 11 coupled differential equations representing various quantities for a single cell. I am then looking at a "grid" of cells (a cell tissue). Currently, each cell is independent from all other cells, but each cell does differ with respect to a single parameter that depends on the location of the cell in the tissue. Since we have independent cells, what I have been doing is just using NDSolveValue to solve the set of 22 equations (11 deqs and 11 initial conditions) multiple times, once for each cell. Currently I am looking at about 250 cells, so I am calling NDSolve 250 times, each time giving it 22 equations.

I now want to include a term in one of the differential equations that actually couples certain cells together. The basic form of the term just looks like


where fluxIP3 is what is coupling the cells together, and gjPerm is just a parameter that sets the strength of this coupling. Because of this new term, I now have to specify 11 separate dependent variables for each cell. Therefore, the new way I am using NDSolveValue is giving it 250*22 equations (number of cells*number of equations) all at once, but only calling it one time. To check that my new code is the same as the old code, I have set the gjPerm parameter to 0. So the only difference between the old and new methods is that I am solving the equations for all of the cells at once rather than solving the equations multiple times for each cell.

The issue I am having is that some (not all) of the solutions for a cell are different between the old and new methods. Below are three examples. One where the solutions match, and two where they do not. The plots are a signal value with respect to time (the independent variable in the differential equations).

Good example Bad Example 1 Bad Example 2

Plotting the maximum difference between the new and old solutions for each cell shows how some solutions match and some do not. The parameter that differentiates one cell from another is determined randomly, so there is nothing to be gained from looking at which cells are messing up just in terms of the cell number on the horizontal axis.

Comparison of solution differences

I have tried making the new code from the old code multiple ways, but I have gotten the same results. Therefore, I do not believe that I have made a mistake in the code. The next thought I am having is that Mathematica is solving the equations differently since the old way only involves solving 22 equations whereas the new version solves 250*22 equations. As evident from the examples, it seems like differences arise in cells where this signal is too high for too long. How can I check how Mathematica is treating these equations, and does anyone have any advice as to solving multiple differential equations? I just need a good way to "look under the hood" I guess. But maybe I am approaching this all wrong, so really anything is appreciated.


I am also wondering which numeric solution is closer to the actual solution. I am guessing that the solution where I pass fewer equations to NDSolve is better, but if someone thinks otherwise let me know.

  • 2
    $\begingroup$ It is possible that there is a difference in the step of integration. For a system from a large number of equations, the step is selected more coarse to save the resource. It is necessary to equalize the methods and steps, then it will be possible to compare. We do not see the code, so we can not tell, that it is necessary to change. $\endgroup$ Commented Jul 5, 2018 at 16:28
  • $\begingroup$ @AlexTrounev This makes sense. Is there a way I can explicitly see what step size is being used automatically in each case, or do I just need to explicitly state it to be the same in each case? As for seeing the code, there is way too much to put on here. $\endgroup$ Commented Jul 5, 2018 at 16:34
  • $\begingroup$ Use the option StepMonitor in each case. $\endgroup$ Commented Jul 5, 2018 at 18:14
  • $\begingroup$ @AlexTrounev Ok I will look into it thanks! Do you think there could there be other issues at play besides different step sizes? $\endgroup$ Commented Jul 5, 2018 at 18:24
  • $\begingroup$ @AlexTrounev Actually, I think in the difference in the steps tends to be the reverse, that the system with the large number of equations has more steps, although I believe it could go either way depending on the distribution of errors. However, I do not think it tries to save resources; the docs for NDSolve indicate otherwise. $\endgroup$
    – Michael E2
    Commented Jul 5, 2018 at 23:18

1 Answer 1


General reasons for a difference

There are (at least) two differences why solving independent systems with one call might result in different solution than solving them separately:

  1. Independent vs. common time-stepping. When solved separately, the time steps are adapted to each system according the error estimates of each system separately.

  2. Different error estimations. NDSolve uses ScaledVectorNorm to estimate the error of a time step of the whole system. The norms of subsets of the errors are different than the norm of the errors altogether.

These are easy to verify with MWEs.


1. Different time-stepping. Here's a nonlinear system designed to oscillate at rates depending on the magnitude of a variable's next neighbor. The system consists of 10 groups of 11 variables. In each group, a variable depends on the next one in cyclical order (cf. Mod[i+ 1, 11, 1] in the code).

One can see in the graphic below that the solutions for x[k] in each group have the same number of time steps. When the groups are integrated separately, the number of time steps varies across groups.

eqs = Table[
   x[11 j + i]''[t] + i/2*x[11 j + Mod[i + 1, 11, 1]][t]^2 x[11 j + i][t] == 0,
   {j, 0, 9}, {i, 11}];
ics = Table[{x[11 j + i][0] == j/10 + i/11, x[11 j + i]'[0] == 0},
   {j, 0, 9}, {i, 11}];

(* Default precision/accuracy parameters *)
pg = Automatic;
wp = MachinePrecision;
ms = Automatic;

solEA = Flatten@MapThread[
    NDSolve[{##}, #2[[All, 1, 1, 0]], {t, 0, 10}, 
       PrecisionGoal -> pg, WorkingPrecision -> wp, MaxSteps -> ms] &,
    {eqs, ics}];

solALL = Flatten@
   NDSolve[{eqs, ics}, Flatten@ics[[All, All, 1, 1, 0]], {t, 0, 10}, 
    PrecisionGoal -> pg, WorkingPrecision -> wp, MaxSteps -> ms];

  Table[x[k]["Grid"], {k, 110}] /. {solEA, solALL}, {2}], 
 PlotLegends -> {"Separately", "All together"}]

Mathematica graphics

2. Different error estimations. The default measure of error is the ScaledVectorNorm based on a $p$-norm with $p=2$. Thus it averages out the errors of the components. It is "scaled" so that an error norm of less than 1 indicates that the error is within the goal set by PrecisionGoal and AccuracyGoal. In this example, I have the equivalent of PrecisionGoal -> 8 and AccuracyGoal -> 8 (which correspond to the 10^-8).

svn = NDSolve`ScaledVectorNorm[2, {10.^-8, 10.^-8}];  (* {pg, ag} *)
ref = RandomReal[{10, 50}, 110];
err = 10^RandomReal[{-9, -8}, Length@ref] ref;
svn[err, ref]
ref = Partition[ref, 11];
err = Partition[err, 11];
MapThread[svn, {err, ref}]
Mean[%^2] // Sqrt
  0.415799                                             <-- Overall error
  {0.442989, 0.367042, 0.298286, 0.415273, 0.469741,   <-- Group errors
   0.304592, 0.32835, 0.543487, 0.487843, 0.426483}
  0.415799                                             <-- Avg. of gps. == overall

What to do about it?

It's taken some time to think about how best to present some ideas for investigating the OP's problem. I learned some things which made the task of explaining more difficult. While my example system shows some superficial resemblance to the OP's images (varying oscillatory behavior), there's little other evidence to suggest that understanding my system's particular problems have any bearing on the OP's system. So I will keep my suggestions general.

Increase WorkingPrecision. A good, accurate method of solution should be stable and converge under increasing settings of PrecisionGoal and AccuracyGoal. In Mathematica, this can be conveniently accomplished by raising WorkingPrecision (the goals are set to one half the value). I tried these on my system (by "both," I mean the separate solEA approach and the all-at-once solALL approach):

WorkingPrecision -> 32  (* both different from MachinePrecision *)
WorkingPrecision -> 40  (* both very close to W.P. -> 32 *)
WorkingPrecision -> 60  (* error test failures NDSolve::nderr *)

The last result is disconcerting, but I ultimately accepted the W.P. -> 40 solution as a good approximation of the true, exact solution. One of the things I learned comes from the docs for NDSolve::nderr:

If you see this message in an example with WorkingPrecision greater than MachinePrecision try using the option Method->"StiffnessSwitching" in NDSolve.

Try a different integration Method. The default ODE method is LSODA. I suspected two things, small errors in early oscillation might be magnified over a long integration period and stiffness. The first should be addressed by WorkingPrecision. My system's stiffness seems not so bad to LSODA, which switches between a nonstiff and stiff solvers (Adams and BDF resp.), that it detects stiffness it cannot handle. I don't know that much about LSODA. The docs state concerning the "NormBound" approach that LSODA uses for "NonstiffTest", "Norm bounds for dense matrices overestimate and the bounds become worse as the dimension increases." Maybe that applies here, since we have a high-dimensional system. I tried the following alternatives:

Method -> "ExplicitRungeKutta"  (* fails, NDSolve::nderr suggests stiffness *)
Method -> "StiffnessSwitching"  (* works well, both agree *)

"ExplicitRungeKutta" does poorly in the presence of stiffness, and it failed quickly on my system. (Only later I discovered NDSolve::nderr had documentation.) "StiffnessSwitching" gave good results for both the solEA and solALL approaches: They were close to each other and close to the "exact" solution obtained with WorkingPrecision -> 40.

Use a different NormFunction. Using the $\infty$-norm suggested itself from the explanation of the difference in error estimation. With the $\infty$-norm, the error estimate would be the greatest scaled error of the components instead of the RMS average. The step size would be governed by the worst component at any given time. However, stiffness turned out to be the dominant factor in my system, and the $\infty$-norm helped only a little. Interestingly, it reversed which approach, separate integration or all together, did better on the worst solution with LSODA, x[83]. Use the $\infty$-norm with the option:

NormFunction -> Infinity

Note NormFunction has nothing to do with the "NormBound" approach LSODA uses for stiffness switching mentioned above.

My example. It's not important, except perhaps to hint there are things to explore, but of all the things I tried using MachinePrecision, the best solution was produced with the following options:

Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}},
 NormFunction -> Infinity
  • $\begingroup$ What are MWEs? I can't seem to find this in the link. $\endgroup$ Commented Jul 5, 2018 at 22:21
  • $\begingroup$ @AaronStevens Minimal Working Example (second bullet point of the link in the answer). Sorry about the jargon. Basically I meant to imply that I cannot verify what the specific issues are with your code, but that these two differences can be observed with simple examples. $\endgroup$
    – Michael E2
    Commented Jul 5, 2018 at 22:55
  • $\begingroup$ Oh ok. Yes I understand this. The problem is that I really can't reduce this down to a simpler example. $\endgroup$ Commented Jul 6, 2018 at 0:25
  • $\begingroup$ How do I look into how NDSolve is using these things in each case? $\endgroup$ Commented Jul 6, 2018 at 18:14
  • $\begingroup$ Thank you for the addition. Things are making a little bit more sense now. I have been messing around with these things but still cannot get the solutions to match. I look forward to what you have to say towards investigating this. $\endgroup$ Commented Jul 9, 2018 at 17:37

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.