Numerically accuracy problem in solving ODEs when subtracting two close large numbers

1 view (last 30 days)
Hi,
I'm trying to numerically solve a set of ODEs and got a problem with numerical stability of my solutions. Simplified version of my actual problem (sufficient to illustrate the problem) consists of 2 ODEs for functions n(t) and nc(t). The equations to be solved are:
  • dn/dt = -alpha*(n+nc)*n + nc/tau + gamma*n
  • d(nc)/dt = alpha*(n+nc)*n - nc/tau + gamma*nc
where alpha, tau, gamma and nu are constants set numerically. Both n(t) and nc(t) grow exponentially due to the last term while the first two terms represent a model of conversion between two species. Specifically, I have alpha=1e-9, tau=150, gamma=1e-5, nu=9000. Initial conditions are n(0)=100, nc(0)=0. I'm using a SimBiology toolbox that uses MATLAB's standard ODE solvers.
As far as I can tell, the problem comes from the fact that the first two terms always add to approximately zero (and I want them in the model because later I'd like to choose parameters such that they wouldn't add to zero). This is because the time scale at which the steady state between the first two terms is reached is shorter than 1/gamma.
Just as an example: at some point in time, alpha*(n+nc)*n would be about 10^8 and -nc/tau would be about -10^8. The problem then comes from the fact that we're subtracting two very large numbers that are very close to each other. And as the error in the estimate of their sum grows it eventually becomes so large as to produce non-sensical results: both n(t) and nc(t) eventually stop growing (and they should continue exponentially).
Since it's not a very demanding simulation, ideally I'm looking for a way to improve the precision in calculating n(t) and nc(t). I've tried changing the MaxStep size, Relative and Absolute Tolerance, but without success.
I appreciate any suggestions.
Thank you.
  1 Comment
Jan
Jan on 10 Feb 2013
When you post the Matlab code of the function to be integrated, it would be easier for the readers to perform some tests.

Sign in to comment.

Answers (1)

Arthur Goldsipe
Arthur Goldsipe on 9 Feb 2013
First off, you mention a parameter nu, but it's not referenced anywhere in your equations. I'm guessing it is a parameter that is not necessary for the simplified problem.
Second, I don't think it's true that n and nc both necessarily continue to grow. When I solved the equations numerically, nc continued to grow approximately exponentially while n had a value around 1/(alpha*tau) = 6.667e6. This is the "magic" value of n that causes some cancellation of terms. I wonder if you could use this fact to recast the equations in terms of n2 = n - (1/alpha*tau) and analyze the behavior of the equations around n2 = 0.

Communities

More Answers in the  SimBiology Community

Categories

Find more on Perform Sensitivity Analysis in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!