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

4 次查看（过去 30 天）
I
I 2013-2-9
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 个）

Arthur Goldsipe
Arthur Goldsipe 2013-2-9
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.

社区

更多回答在  SimBiology Community

类别

在 Help Center 和 File Exchange 中查找有关 Scan Parameter Ranges 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by