Exercises#

The first exercise was to modify the first-order rate law code to account for a y-uncertainty of 5 %.

import numpy as np
import matplotlib.pyplot as plt

Let’s create a new array with the y-uncertainty (c_err).

c = np.array([6.23, 4.84, 3.76, 3.20, 2.60, 2.16, 1.85, 1.49, 1.27, 1.01])
c_err = c * 0.05
t = np.array([2, 4, 6, 8, 10, 12, 14, 16, 18, 20])
plt.errorbar(t, c, c_err)
plt.xlabel('Time/s')
plt.ylabel('[H$_2$O$_2$]/M')
plt.show()
../_images/c159e19a52db4f30f3f5018236fe0a7c517ca80022fed688b72b503c6251cae0.png

Using the model from the example.

def first_order(t, a0, k):
    """
    The first-order rate equation. 
    
    Args:
        t (float): Time (s).
        a0 (float): Initial concentration (mol/dm3).
        k (float): Rate constant (s-1).
    
    Returns:
        (float): Concentration at time t (mol/dm3).
    """
    return a0 * np.exp(-k * t)

But with a modified chi_squared.

def chi_squared(x, t, data, error):
    """
    Determine the chi-squared value for a first-order rate equation.
    
    Args:
        x (list): The variable parameters.
        t (float): Time (s).
        data (float): Experimental concentration data.
        error (float): Uncertainty in concentration data.
    Returns:
        (float): chi^2 value.
    """
    a0 = x[0]
    k = x[1]
    return np.sum((data - first_order(t, a0, k)) ** 2 / (error ** 2))

Then we do the minimisation as before. But with the c_err as an arg and starting positions informed by the previous analysis.

from scipy.optimize import minimize
result = minimize(chi_squared, [7, 0.1], args=(t, c, c_err))
result.x
array([7.08438448, 0.09737419])
x = np.linspace(2, 20, 1000)

plt.plot(t, c, 'o')
plt.plot(x, first_order(x, result.x[0], result.x[1]))
plt.xlabel('Time/s')
plt.ylabel('[H$_2$O$_2$]/M')
plt.show()
../_images/6b7e4d86822b4fd4041e914b683dba5dc27e55b6c9433c031cba4b045205a09a.png

This can also be achieved with the use of the sigma flag in the scipy.optimize.curve_fit function.

from scipy.optimize import curve_fit
popt, pcov = curve_fit(first_order, t, c, sigma=c_err)
/tmp/ipykernel_2482/2319394663.py:13: RuntimeWarning: overflow encountered in exp
  return a0 * np.exp(-k * t)
uncertainties = np.sqrt(pcov)

print(f"[A]_0 = {popt[0]} +/- {uncertainties[0][0]}; k = {popt[1]} +/- {uncertainties[1][1]}")
[A]_0 = 7.084387204820056 +/- 0.17113707059226646; k = 0.09737422484062833 +/- 0.0019530371481279983