# Best step in Runge Kutta 4

I'm simulating radioactive decay of an atom and I have a system of 6 differential equations that i'm solving numerically with runge kutta 4.

I have to simulate that decay for about 40 days which gives ~3e+6 seconds (independent variable). The thing is: of course I'm not gonna use a step of, let's say, 1 second, otherwise I'd have too many points. On the other hand, I'm afraid that if I use a step of ~2 hours (7200 seconds) I will increase error.

So my question is: How do I calculate the best step? Ie, the step that gives me a good precision and a moderate number of points (given the total time that I'm studying).

Thanks!

• you could use $\ln(3,000,000)\approx 15$
– JMP
Dec 21, 2015 at 16:58
• Why the logarithm?
– AJHC
Dec 21, 2015 at 17:04
• so the data fits in with the PNT
– JMP
Dec 21, 2015 at 17:05
• Has it to be RK4? For stiff problems, an implicit method usually works better for larger step sizes. Try if something improves when changing to an exponential parametrization $x_k=\exp(u_k)$, $\dot u_k=\dot x_k/x_k$. This will prevent the $x_k$ crossing into negative territory and could improve the numerical results by reducing the gap in scales of the variables. Dec 22, 2015 at 10:58

That ignores the point that there is nothing special about one second unless you know something about the time scale of your problem. You presumably have six different half lives, which are what set the time scale of the problem. If they are all of the same magnitude, you can take a few dozen steps per half life and be in good shape. The problem comes if you have some very short half lives. If you have a species that lasts $10^{-12}$ seconds, you would be up to $10^{18}$ loops. One approach is to remove those from the model and compute the composition analytically, assuming they disappear "instantly". Say you have decay chains $A \to B \to C, \text { and } D \to B \to C$ to model and $B \to C$ is very fast. You can just model $A \to C, D \to C$ and the amount of $B$ is just the amount of $A,D$ that disappeared in the time step times the ratio of the B decay time to the time step. Now your step size can stay reasonable.