[sundials-users] Numerical solution decreases even though ODE RHS is non-negative

[Eliott Tixier]

Hello everyone!

Consider the following example:

dy1/dt = max(0, 10*y2-y1)
dy2/dt = (1/3600) * (0.8 * sin(t/3600/24) - y2
y1(t=0) = 0.45
y2(t=0) = 0.5

Solving this with Sundials CVODE with tolrel=tolabs1=tolabs2=1e-6 and tsteps = 0, 3600, ..., 3*24*3600, one gets this solution for y1:

As you can see, y1 decreases at some point even though the RHS of its ODE is non-negative at all times.
With tighter tolerances (tolrel=tolabs1=tolabs2=1e-8 for instance), one gets the expected result, i.e. y1 plateaus at 8.

My interrogation is two-fold:

  1. Is this expected? The error on y1 after 72 time steps (264 internal ones) of the simulation is considerably higher than the tolerance of 1e-6. In my (possibly flawed) understanding, the global error is not expected to be many orders of magnitude greater than (abstol + |y|*reltol)*num_steps.

  2. If the answer to 1 is yes, what set of measures can a user of CVODE apply to understand the source of numerical error? In the attached source file we did this:
```
    retval = CVode(cvode_mem, tout, y, &t, CV_NORMAL);
    CVodeGetCurrentTime(cvode_mem, &tcur);
    CVodeGetLastStep(cvode_mem, &tstep);
    CVodeGetCurrentOrder(cvode_mem, &qcur);
    CVodeGetErrWeights(cvode_mem, eweight);
    CVodeGetEstLocalErrors(cvode_mem, ele);
    printf("At asked_t = %0.4e y1 = %0.6e current_t = %0.4e step = %.1e order = %d err1 = %0.1e \n",
      t, Ith(y, 1), tcur, tstep, qcur, Ith(ele, 1));
```
which prints for each time step:
```
At asked_t = 1.6920e+05 y1 = 7.520222e+00 current_t = 1.7744e+05 step = 8.4e+03 order = 4 err1 = 2.6e-06
```
Beside the fact that the internal time step seems large and the order high, this does not give out many clues. 

Thank you in advance for your help and many thanks for developing and maintaing Sundials :)

To reproduce the above example (source file in the attachment), using sing Sundials 7.1.1:
$ gcc example.c -lsundials_cvode -lsundials_core -lm
$ ./a.out

This message is intended only for the personal and confidential use of the designated recipient(s) named above. If you are not the intended recipient of this message you are hereby notified that any review, dissemination, distribution or copying of this message is strictly prohibited. Email transmission cannot be guaranteed to be secure or error-free. Therefore, we do not represent that this information is complete or accurate and it should not be relied upon as such. All information is subject to change without notice.

---

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

[Eliott Tixier]

Hi Dave,

You are right, I did accidentally omit the initial step in the plot. Here is the corrected version of the first plot (tolrel = tolabs1 = tolabs2 = 1e-6).

For the sake of completeness, I've also added the plot for tolrel = tolabs1 = tolabs2 = 1e-8.

On Thu, 12 Sept 2024 at 19:40, David Lorenzetti <dmlore...@lbl.gov> wrote:

Your plot of y1 doesn't honor the initial condition, either (i.e., plot shows y1 starts at 5, not at 0.45).  This suggests there's something up with your code.

-Dave

[Roberts, Steven Byram]

Hi Eliott,

 

The discontinuity in the right-hand side from the max function may be the culprit. The underlying multistep methods in CVODE assume the solution is several times differentiable which is not the case. In the attachment, I modified your code slightly to enable rootfinding on that discontinuity, i.e., when 10*y2-y1=0. When the root is detected, CVODE is reinitialized. This gives the correct asymptotic value even for loose tolerances.

 

I’ll also not that most integrators do not give a monotonic solution when the RHS is nonnegative. These deviations from positivity will be on the order of the local truncation error, though, which is much smaller than the deviation you saw from the discontinuity.

 

Steven

 

From: sundials-users <sundial...@llnl.gov> on behalf of Eliott Tixier <eliott...@NOVAINSILICO.AI>
Date: Friday, September 13, 2024 at 8:34 AM
To: sundials-users <sundial...@llnl.gov>
Subject: Re: [sundials-users] Numerical solution decreases even though ODE RHS is non-negative

Hi Dave,

You are right, I did accidentally omit the initial step in the plot. Here is the corrected version of the first plot (tolrel = tolabs1 = tolabs2 = 1e-6).

For the sake of completeness, I've also added the plot for tolrel = tolabs1 = tolabs2 = 1e-8.

 

 

On Thu, 12 Sept 2024 at 19:40, David Lorenzetti <dmlore...@lbl.gov> wrote:

Your plot of y1 doesn't honor the initial condition, either (i.e., plot shows y1 starts at 5, not at 0.45).  This suggests there's something up with your code.

 

-Dave

This message is intended only for the personal and confidential use of the designated recipient(s) named above. If you are not the intended recipient of this message you are hereby notified that any review, dissemination, distribution or copying of this message is strictly prohibited. Email transmission cannot be guaranteed to be secure or error-free. Therefore, we do not represent that this information is complete or accurate and it should not be relied upon as such. All information is subject to change without notice.

---

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

[Eliott Tixier]

Thank you Steven for your answer and the updated code!

> I’ll also not that most integrators do not give a monotonic solution when the RHS is nonnegative. These deviations from positivity will be on the order of the local truncation error

Yes we had that in mind indeed. Like you pointed it out, our bewilderment was caused by the magnitude of the deviation.

> I modified your code slightly to enable rootfinding on that discontinuity, i.e., when 10*y2-y1=0. When the root is detected, CVODE is reinitialized

Thank you for that. Since we work with symbolic expressions, we will automatically add such root functions based on the non-differentiable terms of the ODE RHS. 

Eliott