In the previous section we have introduced the
Euler method
for the solution of the simple differential equation of first order
that describes the radioactive decay of a number of atoms. In this
section we discuss how to translate this method into an actual
program. In fact, there are many different programming languages that
would do the trick; here and in the following we will concentrate on
programming in Python. For the structuring and design of the programs,
however, we will write some "pseudo-code". This is a rather descriptive
way of formulating the solution of the problem before going into the
linguistic details. Therefore, at the level of pseudo-code, the
language-specific semantics is irrelevant.
This pseudo-code should:
Structure the solution of the problem into sections
that may be represented by functions, routines, or any
other language-specific structures.
Identify all necessary ingredients (algorithms, equations,
parameters etc.).
Detail all steps necessary at each stage of the program,
represented by the aforementioned sections. This then serves
as a guideline for the actual programming work.
In particular this step involves details on all parameters etc.
needed to solve the problem, and it also indicates the flow of
these parameters through the final program.
Enforce the essential rule: Think before you write!
In the case at hand, this is comparably simple.
The essential equation for the solution of our decay problem is
This equation already contains all necessary parameters and variables:
The number of atoms N(t), discretised in time:
N(0), N(Δt), N(2Δt),
N(3Δt), ... .
This will be put into a list or an array (native Python's strength
definitely lies in list-manipulations, but very good arrays are also
provided by numpy) - we will choose arrays.
The discretised time t = 0, Δt,
2Δt, ... .
Again, also this will be put into a list or array (in the latter,
the integer multiplying Δt - 0, 1, 2, ... - would
indicate the position in the array).
The time step Δt, a floating point number.
The time constant τ, a floating point number.
From an algorithmic point of view it is clear that after fixing the
initial values t0=0 and
N(t0=0) =N0
the equation above is repeatedly applied until some final time has been
reached. This automatically fills the two lists of numbers of atoms
and times.
Main program
In the main program, the steering unit, the following
tasks will have to be done:
Input the parameters of the program
Initialise a class Radioactive, encapsulating the physics
problem and a class DEq_Solver for the actual calculation.
Use the latter to solve the differential equation and produce
results.
Store and plot the results.
Initialisation of the physics problem
(encapsulated in the class "Radioactive")
The following parameters etc. need to be initialised here:
the initial time $t_0$, with the obvious default value 0;
the initial number of atoms, $N_0$;
the final time $t_{\rm end}$;
the time constant $\tau$.
The time step δt is not part of the physics
problem.
The class should contain a method to yield the r.h.s. of the defining
differential equation
and for testing purposes also an analytical solution will have to be
provided.
Calculation: Solution of a general differential equation
of first order (encapsulated in the class "DEq_Solver")
The calculation amounts to the following algorithm:
Repeat time steps with size Δt until the
final time is reached.
In the case at hand here, identify the number of atoms $N$
with the $x$, which of course is not a vector then, i.e.
$N\leftrightarrow x$, and the function $f$ on the r.h.s.
becomes $f(N,t)=-N(t)/\tau$.
Store the results of each step and return an array with
times $t_i$ and variables $x_i$.