FD calculation of radioactive decay

We will solve for radioactive using explicit and implicit FD schemes.

Load libraries and set resolution of output figures

[1]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rcParams['figure.dpi']= 300

Define physics constants, here we use arbitrary values without any physical meaning

[2]:
c0      = 1      # initial concentration
k       = 0.1    # decay constant

Setup numerics

[69]:
dt      = 50                            # time step
steps   = 11                            # number of time steps
Cana    = np.zeros(steps)               # analytical solution
Cana[0] = c0                            # initial condition

Cexp    = Cana.copy()                   # make a copy
Cimp    = Cana.copy()
Ccen    = Cana.copy()

Time    = np.linspace(0, (steps-1)*dt, steps)  # time stepping, and Time vector

Perform calculation

[70]:
# hide: the code in this cell is hidden by the author
[71]:
# write update rules for the explicit and implicit cases!
# Bonus: write an update rule that uses the concentration at the middle of the time step
#for n in range(0,len(Cana)-1):

    # Cana[n+1] = np.exp(-k*Time[n+1])            #analytical
    # Cexp[n+1] = ???
    # Cimp[n+1] = ???

Plot different solutions

[72]:
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=(10,5))
fig.clf()
ax  = plt.axes(xlim=(0, dt*(steps-1)), ylim=(-c0, c0))
line, = ax.plot([], [], lw=1)
ax.set_xlabel('Time')
ax.set_ylabel('Concentration')
plt.plot(Time, Cana, label='Analytical solution')
plt.plot(Time, Cexp, label='Explicit solution')
plt.plot(Time, Cimp, label='Implicit solution')
plt.plot(Time, Ccen, label='Centered solution')

plt.legend()

plt.show
[72]:
<function matplotlib.pyplot.show(close=None, block=None)>
../../_images/lecture1_jupyter_radio_decay_11_1.png
[ ]:

[ ]: