1
$\begingroup$

This is not a question which I estimate is aided in answering by providing complete code, which is excessively long.

I ran an NDSolve[] at MachinePrecision, with the specification MaxSteps -> ∞. I then plotted the step size used:

ListPlot[stepsd, PlotRange -> All]

Step size

I thought the step size was disconcertingly high, especially since I was trying to eliminate numerical errors in the output yet the NDSolve[] ran very fast:

AbsoluteTiming[s = NDSolve[Flatten[{xf, yf, zf,vel, pos}], Flatten[{xt, yt, zt}], {t, 0, time}, MaxSteps -> ∞];]

{0.748079, Null}

Consequently I decided to make a graph of the end numerical error (I know the real values to 17 digit precision at the end point) in relation to MaxStepSize ->α: Note that this is error of the form $$\sigma=\sqrt{\sigma^2_x+\sigma^2_y+\sigma^2_z}$$

ListLogLinearPlot[Transpose[{α, σ}]]

error per max stepsize

As you can see the error plummets to about one third, and stays there when the MaxStepSize is reduced to 0.5 or lower, at which point I suspect the model error starts to dominate (the model is of course not perfect). The timing for 0.5 is:

{1.20144, Null}

As the starting values used for this part are between $O[1]$ and $O[-4]$, and are given with 17 digit precision (though Mathematica registers the dataset at MachinePrecision, I'm surprised this numerical error is allowed as the documentation centre claims that:

NDSolve adapts its step size so that the estimated error in the solution is just within the tolerances specified by PrecisionGoal and AccuracyGoal.

Am I missing something, or did the NDSolve[] fail to accurately pick an appropriate MaxStepSize?

Improve this question
$\endgroup$
9
  • $\begingroup$ Can you at least mention what kind of ODEs are you solving? $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Jul 24, 2016 at 14:43
  • $\begingroup$ It's not clear the 2nd plot indicates the condition "that the estimated error in the solution is just within the tolerances specified by PrecisionGoal and AccuracyGoal" is violated. (1) NDSolve uses NDSolve`ScaledVectorNorm[2, {10.^-pg, 10.^-ag}], (pg, ag = PrecisionGoal, AccuracyGoal settings resp.) to measure the acceptability of the estimated error. (2) The estimated error is an estimate, and not guaranteed to bound the actual error. -- Some clarification about how σ is computed would be helpful. $\endgroup$
    – Michael E2
    Commented Jul 24, 2016 at 15:10
  • $\begingroup$ @J.M. The solar system model mentioned in mathematica.stackexchange.com/questions/121035/… Though I've converted to AU and day units. $\endgroup$
    – Feyre
    Commented Jul 24, 2016 at 15:39
  • $\begingroup$ @MichaelE2 I calculated the distance {at a time period of a body (specifically Earth) } between the expected position, and the position that the model computes at that point. It was my (maybe naive) belief that in theory the error should not exceed Precision which is half WorkingPrecision. Am I to understand that with the available information it cannot be guaranteed that the model delivers 7-digit precision as I expected? $\endgroup$
    – Feyre
    Commented Jul 24, 2016 at 15:48
  • $\begingroup$ @Feyre (1) What controls the step size depends on which is bigger, 10.^-pg * * solution or 10.^-ag: In the first case, relative error will dominate, and in the second, absolute error will dominate. When they're about the same, then the (absolute) error tolerance will be about two times 10^-ag. (This is just a consequence of the formula in the tutorial I linked above.) So the answer to your question depends on whether you're talking about absolute or relative error (if the solution has a numerical magnitude of about 1., then they're roughly the same). $\endgroup$
    – Michael E2
    Commented Jul 24, 2016 at 16:39

0

Reset to default

Your Answer

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