Course 2.1Lesson 5 of 5

Project: the half-life with SciPy

SciPy has ready tools for the methods of this chapter. scipy.optimize.brentq solves an equation f of x equal to 0 on a bracket with a change of sign, as fast as the best methods and as safe as bisection. scipy.optimize.curve_fit finds the numbers of any model function that fit measured data best, and returns the best values and a covariance matrix whose diagonal gives their uncertainties. With a sigma of the square root of the counts, a fit of the exponential decay n0 times e to the minus lambda t to the counts of a radioactive sample gives the decay constant, and the half-life is ln 2 over lambda, about 20 seconds in the example. The same half-life comes from a straight line through the logarithms of the counts, and brentq finds the time when the counts reach 100.

about 35 min1 simulation

Lesson 5 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

solve an equation on a bracket with brentq,
fit any model to data with curve_fit and read the uncertainties, and
combine the tools of the whole chapter in a project.

SciPy

SciPy is the library of scientific tools that is built on NumPy. It has methods for integration, optimisation, statistics and much more. You use it in the same way as the other libraries, and on this site it is ready. The first time that your code imports it, the site needs a moment to load it. We use the part scipy.optimize.

Solving an equation: brentq

Predict first. The first equation is the one that you solved by bisection. What will the answer be? Then run it.

brentq(f, a, b) finds a root of f on the bracket $[a, b]$ . Like bisection, it needs a change of sign. Unlike bisection, it uses faster methods inside, and takes few steps, but it keeps the root inside the bracket, so it cannot fly away like Newton's method.
The first call finds the root 2.094551 of $x^3 - 2x - 5$ . The second finds the time when the ball lands, 4.179 s. The third solves $\cos x = x$ , which has no formula: $x = 0.739085$ , and the check shows that $\cos x$ is also 0.739085. This number is a famous constant: press the cosine button on a calculator, in radian mode, again and again, and the display goes to it.
A lambda is a short function without a name, written in one line: lambda x: x ** 3 - 2 * x - 5.

Rule of thumb: to solve an equation in one unknown, find a bracket, and call brentq.

Fitting a model to data: curve_fit

In the last lesson, polyfit fitted a polynomial. But many laws are not polynomials. Radioactive decay follows $N = N_0 e^{-\lambda t}$ , where $N$ is the count at time $t$ , $N_0$ is the count at the start, and $\lambda$ is the decay constant. curve_fit fits any formula to data. You write the formula as a function, and give it the data.

Predict first. The counts fall from 1000 to about half in 20 seconds. What do you expect for the half-life? Then run it.

Read the program.

- model(t, n0, lam) is the formula. The first argument is the variable ( $t$ ), and the others are the numbers to be found ( $N_0$ and $\lambda$ ).

- curve_fit(model, t, counts, p0=[1000, 0.03], ...) finds the best $N_0$ and $\lambda$ . p0 is a starting guess, because the fit is made by successive improvements, in the manner of Newton's method, and it needs a place to start. A guess of the right size is enough.

- sigma=np.sqrt(counts) tells the fit the uncertainty of each count, as in chapter 2: a count of 1000 is uncertain by about 32, and a count of 125 by about 11. The points with the smaller uncertainty count more. absolute_sigma=True says that these are real uncertainties.

- The function gives two things. params holds the best values. cov is the covariance matrix, and the square roots of its diagonal, np.sqrt(np.diag(cov)), are the uncertainties of the numbers.

- Results: $N_0 = 999.1 \pm 25.3$ and $\lambda = 0.0344 \pm 0.0010$ per second. A fitted number is always given with its uncertainty. The half-life, the time in which the count halves, is $\ln 2 / \lambda = 20.2$ s.

The project: the half-life of a sample

Now we put the pieces together to answer three questions about a radioactive sample: what is its half-life, does a second method agree, and when will the count fall to 100 per second?

Predict first. The curve falls to half in about 20 seconds, so after 60 seconds it should be about one eighth of the start. When will it be at 100? Then run it.

A graph titled Decay of a sample. Seven blue dots with small error bars fall from 1000 counts per second at time 0 to 125 counts at 60 seconds. An orange curve, the fit, passes through the error bars and continues to 80 seconds. Dotted grey lines mark 100 counts and the time of about 67 seconds, where the curve reaches 100 counts.

Read the program.

- The first part is the fit, as before. The half-life is 20.2 s.

- A second way: since $\ln N = \ln N_0 - \lambda t$ , the logarithms of the counts lie on a straight line whose slope is $-\lambda$ . np.polyfit(t, np.log(counts), 1) finds it. The half-life from this line is $-\ln 2 / \text{slope} = 20.1$ s. Two different methods agree to within 0.1 s, and that is a good sign. (The log method does not use the uncertainties of the points, so a small difference is normal.)

- The last question is an equation: when is the model equal to 100? brentq(lambda s: model(s, n0, lam) - 100, 0, 200) finds the root of "model minus 100" on the bracket 0 to 200 seconds. At 0 the difference is positive (about 900), and at 200 it is negative (about $-99$ ), so there is a bracket. The answer is 66.9 s.

- The graph shows the data with error bars, the fitted curve, and the two dotted lines that meet at the answer.

Notice how the chapter fits together: the fit gives the model, the solver answers a question about the model, and the graph shows the result. This is how a real analysis goes.

brentq finds the time where the count is half by searching a bracket, like bisection, only faster. It agrees with ln 2 / λ.

A warning about fits

A fit always gives numbers, even for a bad model. Always ask three questions. Does the curve go through the error bars on the graph? Are the uncertainties of the numbers small compared to the numbers? Does the result make sense for the physics (is the half-life positive and in the range that you expect)? A fit that gives a half-life of 20.2 s with an uncertainty of 0.5 s is a result. One with an uncertainty of 20 s says that the data cannot tell.

Common mistakes

Giving brentq an interval without a sign change. It raises a ValueError.
Writing the model with the variable second. In curve_fit, the variable is the first argument.
A starting guess far from the truth. The fit can fail, or give a wrong answer.
Reporting a fitted number with no uncertainty.
Taking the logarithm of a count of zero. It does not exist.
Trusting the fit without looking at the graph.

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 Project: the half-life. You write the fit, the half-life and the time to reach a given count.
Build. Open the problem Half-life of a sample and the problem A different law.

Remember these

brentq(f, a, b) solves $f(x) = 0$ on a bracket. It is safe and quick.
curve_fit(model, x, y, p0=...) fits any model. The variable is the first argument of the model.
The uncertainties of the fitted numbers are np.sqrt(np.diag(cov)).
The half-life is $\ln 2 / \lambda$ .
Fit, solve and plot: this is the pattern of a real analysis. Check each result in another way.