Finite Difference Simulation: Voltammetry in Python

As a programming language Python is many things... but one thing it is not is fast. Does this mean we can't use it in electrochemical simulation? The book `Understanding Voltammetry: Simulation of Electrode Process' presents example programs written in C++ that enable the simulation of electrochemical processes using finite difference methods. The question is, can we translate these example programs into Python and allow the problems to be solved in a reasonable amount of time?

C++ vs. Python

The article is going to (albeit crudely) focus on benchmarking the execution time for solving the one-dimensional diffusion equation through the use of a finite difference method. Here we take the 1D diffusion equation, discretise it using a central finite difference method and solve the problem fully implicitly, see this book for details on the used method. There are certainly examples of people using Python to simulate cyclic voltammetry using explicit finite difference methods. Explicit procedures are less computationally expensive but are only conditionally stable, this stability issue can be problematic when simulating voltammetric responses of an electroactive species with fast coupled homogeneous kinetics. As such in the following we briefly consider if we can perform fully implicit 1D voltammetric simulations in Python, a further post will consider two 2D problems.

The simulation solves the Diffusion Equation:
$$\frac{\partial c}{\partial t} = \nabla^2 c$$

subject to the appropriate boundary conditions. After discretisation of this expression, then for a one-dimensional problem, we are left with a set of equations that can be represented by a tridiagonal matrix. In the textbook example C++ code this tridiagonal matrix is solved for directly using an implementation of the Thomas Algorithm. This implementation has not been optimised but even on an antiquated laptop (Intel Core i5-5200U, 2.2 GHz) the simulation only takes 2.7 seconds.

What happens if we take a naïve approach and directly translate this code to Python?


This isn't good... the program takes ~75 s to run; almost 30 times slower than the C++ implementation. Happily there's a lot of room for improvement. We essentially have two options, first we could use a Python library to do the heavy lifting for us or second we could write our own pre-compiled extension.

As we have a tridiagonal matrix lets look at what the Scipy library has to offer. Scipy provides access to the Lapack banded matrix solver. Using this solver requires use to modify how we set up the matrix but this isn't a major undertaking. If we use this Fortran based package to solve the matrix we can solve the one-dimensional mass-transport problem in only 3.5 s, a big improvement and competitive with the non optimised C++ code.


Given how simple this implementation is... it is pretty good. Even better the Lapack routine is a more general solver allowing us to easily tackle pentadiagonal or even denser matrices. In many cases just using a Python library can sufficiently decrease a programs execution time enough to justify the use of an easier to understand programming language.

We could finish here but lets see if we can get any further with Cython. Cython gives us an easy way by which we can write a compiled extension for Python but in doing so we loose some of the point of using Python itself. We now have to compile our code and we loose some of clarity of just using Python. Still there are cases in which using Cython seems justifiable; if I have a larger Python program but the runtime is limited by a relative small section of the code and there is no obvious Scipy library that can do the task I want, I will often resort to making a bespoke Python extension. Cython requires us to write a separate .pyx file and to pre-compile it before use. For the present cases if we just transfer a lot of the script directly to a Cython file, we will immediately get some improvements in the runtime. However, the main point is to ensure that we define the types for all of our used variables so that Cython can optimise the compiled code. The script below shows an example Cython function where the some effort has been made to optimise the function providing types for all of the variables. Using this compiled extension the simulation only take 1.7 s, this is faster than the none optimised C++ example given in the textbook!


Given the simplicity of the example, the use of Cython does not really seem justified here. As a final point I just want to highlight that another way in which we can improve the runtime is to implement an expanding spatial grid so as to reduce the size of the matrix that needs solving. As a final example I'm providing a fuller example of a 1D voltammetric implicit simulation program which takes dimensional parameters, converts them to a dimensionless form, uses an expanding spatial grid and employs the Butler-Volmer equation to describe the interfacial electron transfer kinetics. Finally the example script also does some basic plotting and compares the simulated peak height to that expected analytically from the Randles-Ševčík equation. The runtime for this simulation is about 0.5 s.