Simulating radioactive decay.

[Ludwik Kowalski <KowalskiL@MAIL.MONTCLAIR.EDU>]

This morning, while responding to a message, I wrote:

> An activity with coins (or m%m candies) can be supplemented
> with a similar activity with pencils. Just declare that "a label
> up means a decay". This will give students an opportunity to
> compare two decay curves, one whose probability of decay
> per unit time (one throw) is 1/2 and another with p=1/6.

In trying to be quantitative I just realized that p should not be
identified with the probability of decay per unit time as it is
defined in the theory of radioactive decay. In other words p is
not the same thing as lambda. Both are "probabilities of decay
per unit time" but the definition of lambda calls for dt which is
negligibly small in comparison with the half-life, T. The
definition of p, on the other hand, calls for dt which is equal
to T (I am referring to simulations with coins, or m&m candies).

The well known formula, T=ln(2)/lambda, would produce
wrong results if p was used instead of lambda, except when
p<<1. In that sense simulations with pencils (p=1/6) are better
than simulations with coins (p=1/2). The decay curves are
exactly exponential for all constant values of p but the formula
for calculating T is valid only when very short dt makes p<<1.

1) How many throws of PENCILS will be needed to reduce the
    initial population to 50% ? The answer ln(2)/(1/6)=4.16
    would certainly disagree with experimental data.

2) How many throws of COINS will be needed to reduce the
    initial population to 50% ? The answer ln(2)/(1/2)=1.38
    would certainly disagree with experimental data.

But the discrepancy in the first case will be much smaller
than in the second case. A little program based on random
numbers should illustrate this. A good student project,
I suppose.

Input the value of p, No and Nmin. (such as 1/6, 10000, 1000)
        a) Process each "piece" according to p
        b) Print the number of remaining pieces N
        c) Make new No=N (or stop when N is less than Nmin)
        d) Go back to (a)

The curve N=f(t) can be plotted to determine the "experimental"
value of T. It will NOT be equal to 4.16 units. (1 unit =1 throw).
By the way, one does not need a program to realize that for p=1/2
the answer 1.38 units is larger than the experimental half-life T.
The discrepancy for p=1/6 will be considerably smaller than 38 %,
perhaps only 5 or so. Can somebody predict the exact value of the
discrepancy without using a program (when p=1/6, for example)?
I can not do this.

Ludwik Kowalski