Course 3.1Lesson 4 of 5

Simulating radioactive decay

In radioactive decay, each atom decays at a random moment, with a probability p in each second. A simulation follows the number of atoms that are left: in each second, rng.binomial of the atoms and p gives how many decay. The result follows the expected curve n0 times one minus p to the power t, but it fluctuates, and the fluctuation of a count is about the square root of the count, so a small sample is rough and a large sample is smooth. With p equal to 0.05 the expected half-life is ln 2 over minus ln of 1 minus p, which is 13.5 seconds, and a simulation that counts whole seconds finds 14. The exponential law of decay is not a rule for single atoms, whose decay is unpredictable, but the result of the average over a huge number of atoms.

about 35 min1 simulation

Lesson 4 of 5Show all lessons

Try the code and the simulations beside the lesson.

What you will be able to do

After this lesson you can

simulate the decay of a sample of atoms with random numbers,
compare a simulation with the exponential law, and
explain why small samples are noisy and large samples are smooth.

A law made of chance

In chapters 2 and 3, you drew the decay of a sample as a smooth curve, $N = N_0 e^{-\lambda t}$ , and fitted it. But what happens to one atom? Nobody can say when an atom will decay. Each atom has only a probability to decay in the next second. The smooth law is what you see when many atoms, each with its own chance, are counted together. A simulation shows how a smooth law grows out of chance.

The model: each second, each atom that is left decays with the probability $p$ , independently of the others. In a sample of $n$ atoms, how many decay in one second? It is a binomial random number, and NumPy has it: rng.binomial(n, p) is the number of "successes" in $n$ trials, each with the probability $p$ .

The simulation

Predict first. Each second, 5 percent of the atoms decay. How many are left after 10 seconds, on average? (Each second, the number is multiplied by 0.95.)

The output compares the two: at 10 s, 6023 atoms in the simulation and 5987.4 expected; at 30 s, 2176 and 2146.4; at 60 s, 466 and 460.7. The simulation follows the exponential law, with small differences that come from chance. This is a way to see why a law of decay is exponential: each second, the same share of the atoms decays, so the number is multiplied by the same number again and again.

Noise: small and large samples

A graph of the share of atoms that are left against time from 0 to 80 seconds. A black dashed curve, the theory, falls smoothly from 1 to about 0.02. A blue line for 10000 atoms lies on the theory. An orange line for 200 atoms wanders a little above and below it and is rough, with small steps.

Predict first. The expected number after 30 seconds is 21.5 for 100 atoms, 2146 for 10 000, and 214 639 for a million. Which sample will give the most regular result?

A bigger sample is smoother. The picture shows the same story: the orange curve for 200 atoms is rough and wanders around the theory, and the blue curve for 10 000 atoms lies on it.

There is a rule behind this. The scatter of a count is about the square root of the count, so the relative scatter is $\dfrac{\sqrt{N}}{N} = \dfrac{1}{\sqrt{N}}$ . For 21 survivors, it is $1/\sqrt{21} = 0.22$ , which is exactly the number that the program found. This is the same $1/\sqrt{N}$ law as in the Monte Carlo method, and it explains why a Geiger counter that counts only 100 particles is accurate to about 10 percent.

The dots are the expected counts. The code's simulated counts scatter around them by about the square root of the count.

The half-life

Run it. The half-life is the time in which half the atoms decay. With $p = 0.05$ , the exact number is found from $0.95^t = 0.5$ , which gives $t = \dfrac{\ln 2}{-\ln(1 - p)} = 13.51$ s. The simulation counts how many whole seconds are needed until half the atoms are gone, and averages 200 runs. It gives 14.01, not 13.51. Why? The simulation can only look at whole seconds. After 13 seconds, there are still more than half the atoms (about 5133), and after 14 seconds fewer than half (about 4877). So the first whole second at which the count is below half is almost always the next whole number, 14. The difference of 0.5 is not an error of the method: it comes from the time steps. A smaller time step gives a result nearer to 13.51.

Common mistakes

Taking a random number of decays for each atom in a loop. It is slow. Use rng.binomial for the whole sample.
Forgetting to update the count of the atoms that are left. The next step must use the new count.
Comparing one run with the theory and calling the difference an error. Run many times.
Thinking that a small sample must follow the law exactly. Its scatter is large.
Forgetting the time step when you compare with a half-life.

Practice

Start with two quick questions. Each New numbers gives you a fresh one, so you can practise as long as you like.

Then the labs and the problems:

Warm-up. Open the lab Simulate the decay. You write the simulation as a function, and a function for the half-life.
Build. Open the problem Simulated half-life and the problem Noise in a count.

Remember these

An atom decays at a random moment. A sample decays by a smooth exponential law, because many chances are averaged.
rng.binomial(n, p) is the number of decays among $n$ atoms in one step.
The scatter of a count of $N$ is about $\sqrt{N}$ , so a larger sample is relatively smoother, by $1/\sqrt{N}$ .
A simulation with whole time steps gives a half-life that is rounded up to the next step.

Check yourself

Try each one on your own first. If you are stuck, open a hint. Then check with Show the answer.

10 000 atoms decay with $p = 0.1$ per second. About how many are left after 3 seconds?

A count of 400 particles has what usual scatter, and what relative scatter?

Why is the decay of a sample exponential?

Why does the simulation give 14 and not 13.51?

Enroll free in Simulation and Randomness

Every lesson is open to read. Enrolling also lets you:

Mark lessons done and carry on where you stopped.
Take the tests of this lesson and its course.
Count the course towards the certificate for its chapter.

No payment. You can leave at any time.

Problems (2)

Each atom of a sample decays with the probability 0.02 in each second.Medium · High School (11–12) · by StemPy Lab

A sample has 400 atoms, and each decays with the probability 0.05 per second.Hard · High School (11–12) · by StemPy Lab

Python Labs (1)

Next: lesson 5Project: how sure is a measurement?

← Simulation and Randomness

← Random walks and diffusion Project: how sure is a measurement? →

Discussion (0)

Loading…