Using lambdify for plotting expressions

The syntethic isotope Technetium-99m is used in medical diagnostics (scintigraphy): $$ ^{99m}Tc \overset{\lambda_1}{\longrightarrow} \,^{99}Tc \overset{\lambda_2}{\longrightarrow} \,^{99}Ru \\ \lambda_1 = 3.2\cdot 10^{-5}\,s^{-1} \\ \lambda_2 = 1.04 \cdot 10^{-13}\,s^{-1} \\ $$ SymPy can solve the differential equations describing the amounts versus time analytically. Let's denote the concentrations of each isotope $x(t),\ y(t)\ \&\ z(t)$ respectively.

import sympy as sym
sym.init_printing()

symbs = t, l1, l2, x0, y0, z0 = sym.symbols('t lambda_1 lambda_2 x0 y0 z0', real=True, nonnegative=True)
funcs = x, y, z = [sym.Function(s)(t) for s in 'xyz']
inits = [f.subs(t, 0) for f in funcs]
diffs = [f.diff(t) for f in funcs]
exprs = -l1*x, l1*x - l2*y, l2*y
eqs = [sym.Eq(diff, expr) for diff, expr in zip(diffs, exprs)]
eqs

$$\left [ \frac{d}{d t} x{\left (t \right )} = - \lambda_{1} x{\left (t \right )}, \quad \frac{d}{d t} y{\left (t \right )} = \lambda_{1} x{\left (t \right )} - \lambda_{2} y{\left (t \right )}, \quad \frac{d}{d t} z{\left (t \right )} = \lambda_{2} y{\left (t \right )}\right ]$$

solutions = sym.dsolve(eqs)
solutions

$$\left [ x{\left (t \right )} = C_{1} e^{- \lambda_{1} t}, \quad y{\left (t \right )} = \frac{C_{1} \lambda_{1} e^{- \lambda_{1} t}}{- \lambda_{1} + \lambda_{2}} + C_{2} e^{- \lambda_{2} t}, \quad z{\left (t \right )} = - \frac{C_{1} \lambda_{2} e^{- \lambda_{1} t}}{- \lambda_{1} + \lambda_{2}} - C_{2} e^{- \lambda_{2} t} + C_{3}\right ]$$

now we need to determine the integration constants from the intial conditions:

integration_constants = set.union(*[sol.free_symbols for sol in solutions]) - set(symbs)
integration_constants

$$\left\{C_{1}, C_{2}, C_{3}\right\}$$

initial_values = [sol.subs(t, 0) for sol in solutions]
initial_values

$$\left [ x{\left (0 \right )} = C_{1}, \quad y{\left (0 \right )} = \frac{C_{1} \lambda_{1}}{- \lambda_{1} + \lambda_{2}} + C_{2}, \quad z{\left (0 \right )} = - \frac{C_{1} \lambda_{2}}{- \lambda_{1} + \lambda_{2}} - C_{2} + C_{3}\right ]$$

const_exprs = sym.solve(initial_values, integration_constants)
const_exprs

$$\left \{ C_{1} : x{\left (0 \right )}, \quad C_{2} : \frac{1}{\lambda_{1} - \lambda_{2}} \left(\lambda_{1} x{\left (0 \right )} + \lambda_{1} y{\left (0 \right )} - \lambda_{2} y{\left (0 \right )}\right), \quad C_{3} : x{\left (0 \right )} + y{\left (0 \right )} + z{\left (0 \right )}\right \}$$

analytic = [sol.subs(const_exprs) for sol in solutions]
analytic

$$\left [ x{\left (t \right )} = x{\left (0 \right )} e^{- \lambda_{1} t}, \quad y{\left (t \right )} = \frac{\lambda_{1} x{\left (0 \right )} e^{- \lambda_{1} t}}{- \lambda_{1} + \lambda_{2}} + \frac{e^{- \lambda_{2} t}}{\lambda_{1} - \lambda_{2}} \left(\lambda_{1} x{\left (0 \right )} + \lambda_{1} y{\left (0 \right )} - \lambda_{2} y{\left (0 \right )}\right), \quad z{\left (t \right )} = - \frac{\lambda_{2} x{\left (0 \right )} e^{- \lambda_{1} t}}{- \lambda_{1} + \lambda_{2}} + x{\left (0 \right )} + y{\left (0 \right )} + z{\left (0 \right )} - \frac{e^{- \lambda_{2} t}}{\lambda_{1} - \lambda_{2}} \left(\lambda_{1} x{\left (0 \right )} + \lambda_{1} y{\left (0 \right )} - \lambda_{2} y{\left (0 \right )}\right)\right ]$$

Exercise: Create a function from a symbolic expression

We want to plot the time evolution of x, y & z from the above analytic expression (called analytic above):

from math import log10
import numpy as np
year_s = 365.25*24*3600
tout = np.logspace(0, log10(3e6*year_s), 500)  # 1 s to 3 million years

%load_ext scipy2017codegen.exercise

Use either the %exercise or %load magic to get the exercise / solution respecitvely:

Replace ??? so that f(t) evaluates $x(t),\ y(t)\ \&\ z(t)$. Hint: use the right hand side of the equations in analytic (use the attribute rhs of the items in anayltic):

# %exercise exercise_Tc99.py
xyz_num = sym.lambdify([t, l1, l2, *inits], [eq.rhs for eq in analytic])
yout = xyz_num(tout, 3.2e-5, 1.04e-13, 1, 0, 0)

import matplotlib.pyplot as plt
%matplotlib inline
fig, ax = plt.subplots(1, 1, figsize=(14, 4))
ax.loglog(tout.reshape((tout.size, 1)), np.array(yout).T)
ax.legend(['$^{99m}Tc$', '$^{99}Tc$', '$^{99}Ru$'])
ax.set_xlabel('Time / s')
ax.set_ylabel('Concentration / a.u.')
_ = ax.set_ylim([1e-11, 2])