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!
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.