Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Example that uses atol to adjust absolute tolerances for CVODE #110

Open
artgoldberg opened this issue Nov 26, 2019 · 3 comments
Open

Example that uses atol to adjust absolute tolerances for CVODE #110

artgoldberg opened this issue Nov 26, 2019 · 3 comments

Comments

@artgoldberg
Copy link

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

@bmcage
Copy link
Owner

bmcage commented Nov 28, 2019

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

@artgoldberg
Copy link
Author

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

@bmcage
Copy link
Owner

bmcage commented Jan 12, 2020

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants