2
$\begingroup$

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!

$\endgroup$
4
  • $\begingroup$ you could use $\ln(3,000,000)\approx 15$ $\endgroup$
    – JMP
    Dec 21, 2015 at 16:58
  • $\begingroup$ Why the logarithm? $\endgroup$
    – AJHC
    Dec 21, 2015 at 17:04
  • $\begingroup$ so the data fits in with the PNT $\endgroup$
    – JMP
    Dec 21, 2015 at 17:05
  • $\begingroup$ 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. $\endgroup$ Dec 22, 2015 at 10:58

1 Answer 1

1
$\begingroup$

How are you doing the computation? If by computer, 3 million loops is nothing. If by hand or a spreadsheet, it is a lot. If you compute 3 million points, you don't have to print/plot them all. You can just print out every 10,000th or so.

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.

The experimental approach is to start with very long time steps, like a day. Repeat the calculation with the time step cut in half. Continue until the results stop changing.

There are two penalties for too short a time step. The first is computation time, which you are thinking about. The other is numeric accuracy. If the time step is too short, the change in number of atoms will be too small to register in the machine accuracy. In this day of 64 bit floats, that is not likely to be a problem, but it used to be.

$\endgroup$
2
  • $\begingroup$ That's exactly my problem. I have 4 half lives of hours and one half life of milliseconds, so if I choose a step of 3600s, I will get negative quantities in the isotope that has half life of ms. I'd like to try out your suggestion of "removing" the fast isotope (to get positive quantities, I hope!), however I did not get this part: "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" Could you explain better? $\endgroup$
    – AJHC
    Dec 21, 2015 at 18:14
  • $\begingroup$ Say you take a time step of 10 seconds, 10^6 atoms of A and 10^5 atoms of D decay in that time step, and the decay time of B is 10^-5 seconds. Then you have 10^5+10^4 atoms of B created every second. That would say you have about 1.1 atoms of B present at each moment. Because the decay time of B is so fast compared to a time step, it decays as fast as it is produced. $\endgroup$ Dec 21, 2015 at 18:19

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .