odes icon indicating copy to clipboard operation
odes copied to clipboard

Example that uses atol to adjust absolute tolerances for CVODE

Open artgoldberg opened this issue 5 years ago • 3 comments

I'm using scikits because it's listed on the CVODE webpage. I'm building a biochemical simulator. Since populations of molecules must be non-negative, I'm following the advice in Hindmarsh, Serban and Reynolds, User Documentation for cvode v5.0.0, 2019, which recommends that tighter absolute tolerances be used to reduce the magnitude of negative values. However, I find that changing atol does not change the solver behavior. I'm passing atol like this:

from scikits.odes import ode
options = {'atol': 1e-11}
CVODE_SOLVER = 'cvode'
solver = ode(CVODE_SOLVER, right_hand_side, old_api=False, **options)

Is this correct? Do you have an example that uses atol to adjust absolute tolerances for CVODE?

Thanks, Arthur

artgoldberg avatar Nov 26 '19 22:11 artgoldberg

This is correct and should work. Note that the default is 1e-12, and rtol 1e-6, so setting 1e-11 might not have effect. To see this action, open the jupyter notebook locally: https://github.com/bmcage/odes/blob/master/ipython_examples/Simple%20Oscillator.ipynb

After the first section, so after [9], you can for example do following command to see obtain a more precise solution:

options1= {'rtol': 1e-6, 'atol': 1e-12, 'max_steps': 50000}      # default rtol and atol
options2= {'rtol': 1e-15, 'atol': 1e-25, 'max_steps': 50000}
solver1 = ode('cvode', rhseqn, old_api=False, **options1)
solver2 = ode('cvode', rhseqn, old_api=False, **options2)
solution1 = solver1.solve([0., 1., 60], initx)
solution2 = solver2.solve([0., 1., 60], initx)

print('\n   t       Solution1       Solution2          Exact')
print('-----------------------------------------------------')
for t, u1, u2 in zip(solution1.values.t, solution1.values.y, solution2.values.y):
    print('{0:>4.0f} {1:15.8g} {2:15.8g} {3:15.8g}'.format(t, u1[0], u2[0], 
           initx[0]*np.cos(np.sqrt(k/m)*t)+initx[1]*np.sin(np.sqrt(k/m)*t)/np.sqrt(k/m)))

Output would be:

   t       Solution1       Solution2          Exact
-----------------------------------------------------
   0               1               1               1
   1     -0.37069371     -0.37068197     -0.37068197
  60       0.8430298      0.84321153      0.84321153

bmcage avatar Nov 28 '19 13:11 bmcage

Thanks, that's helpful. I evaluated the sensitivity of 17 test cases to the tolerances. 2020-01-08_test_tols_all.txt.pdf

artgoldberg avatar Jan 11 '20 19:01 artgoldberg

I'm confused on how to interpret this. Rel tol of 0.0001 has highest compute time ? For failure rate, do you consider max steps reached a failure ? For advanced computation, one should use cvode stepwise, so compute forward one step at a time without using max_steps.

bmcage avatar Jan 12 '20 10:01 bmcage