[sundials-users] Trying to understand error control

[John Weeks]

I am trying to understand how CVode()'s error control relates to actual errors. To that end, I am running CVode() on the simplest system I could think of:
dydx = -0.05 * y
Since I can solve it analytically, I can compare the output to the analytic solution.
I call CVodeSStolerances(cvode_mem, 1e-10, 1e-10).
I am using the Adams-Moulton method.
I set the initial value to 10. I then call CVode() repeatedly to deliver solutions at x = 1, 2, ..., 100.
The analytic solution is y = 10*exp(-0.05 * x). I computed the differences between the CVode() output and the analytic solution, and divided by the analytic solution for a relative error. Here is the result, which shows numerical errors quite a bit larger than 1e-10:

Is there any statement that can be made to connect these real errors with the settings of atol and rtol?

We use CVode() as the engine for a feature in our scientific graphing and numerical analysis application Igor Pro.

Thanks much,
John Weeks
WaveMetrics, a division of Sutter Instrument, inc.

---

To unsubscribe from the SUNDIALS-USERS list: write to: mailto:SUNDIALS-USERS-...@LISTSERV.LLNL.GOV

[Daniel Reynolds]

Hi John,

Unlike the case with hand-constructed tests like yours, the time integrators in SUNDIALS have no knowledge of the solution to an IVP, nor of the true error incurred in a time step.  Instead, they perform a process called "local error control," wherein they must approximate the error incurred at each time step.  If that single-step error exceeds the user-specified tolerances then the step is rejected.  However, if the single-step error is within tolerances the step is accepted, and the adaptivity strategies choose a step size for the next attempt (and in the case of CVODE and IDA, the order of the linear multistep method to use).  Details on how this works in CVODE is given in the documentation (https://sundials.readthedocs.io/en/latest/cvode/Mathematics_link.html#local-error-test).

Thus, while you indeed specify both a scalar-valued relative tolerance, and a scalar- or vector-valued absolute tolerance, those are not requests on the accuracy of the global-in-time solution (as you plot in your graph) -- those are requests for the estimated single-step errors.  Depending on the IVP in question, local errors can compound, resulting in global errors that exceed tolerances, or they can even cancel/decay, leading to global errors at (or even below) tolerances.

To me, your plot looks entirely reasonable given the tolerances, in that the green curve is always well below 1e-9, indicating that the global relative error stays within a factor of 10 from the requested relative tolerance.  I'll note that the way that these tolerances are combined in SUNDIALS (see https://sundials.readthedocs.io/en/latest/cvode/Mathematics_link.html#ivp-solution), a step can essentially be successful when it satisfies either the requested relative or absolute tolerance -- from your results it looks like the relative tolerance is "active" throughout most of the time interval.

Daniel R. Reynolds (he/him)

Professor, Dept of Mathematics and Statistics

University of Maryland Baltimore County, MP 429

410.455.2429

https://drreynolds.github.io

[John Weeks]

Thank you, Daniel, for your thoughtful answer to my question.

I do understand the difference between local errors and global errors. I was hoping that they might be connected by some predictable relationship, like multiplying by step size or number of steps. And it occurs to me now that each point on my graph is the result of integration from the previous point. That is, the previous point is the initial condition for the next.

My goal is to write something like unit tests that will identify when we have broken something about the solver code, or if an update of SUNDIALS breaks something. So I guess the answer is that I can't really count on a given actual global error, but I can count on being able to test the global error in a more generalized way. And perhaps the test has to be written based on empirical evidence, as opposed to what I was trying which was to compute a limit.

John Weeks