A short introduction to digital simulations in electrochemistry : simulating the Cottrell experiment in NI LabVIEW

A brief introduction to the use of digital simulations in electrochemistry is given by a detailed description of the simulation of Cottrell’s experiment in the LabVIEW programming language. A step-by-step approach is followed and different simulation techniques (explicit and implicit Euler, Runge–Kutta and Crank–Nicolson methods) are applied. The applied techniques are introduced and discussed on the basis of Padé approximants. The paper might be found useful by undergraduate and graduate students familiarizing themselves with the digital simulation of electrochemical problems, as well as by university lecturers involved with the teaching of theoretical electrochemistry.


Introduction
At the 6 th Regional Symposium on Electrochemistry for South-East Europe (6 th RSE-SEE), I gave a keynote lecture [1] on the kinetics of hydrogen evolution in mildly acidic solutions.Under these chemical conditions -at low to moderate pH values (1 < pH < 4) -the reduction of H + ions and that of water molecules both give a significant contribution to the measured cathodic current.In my talk, I presented digital simulations (and a simplified analytical model) that can describe polarization curves measured on rotating disk electrodes immersed into mildly acidic solutions.
Results presented in my talk at the 6 th RSE-SEE were published [2] just a few weeks before the conference.Instead of further discussing these results, I chose to focus more in this paper on the applied method of digital simulation (DS).
DS presents a powerful tool that can be used for the qualitative understanding, as well as for the quantitative description of electrode processes.When carrying out DS, we create a numerical model of the electrochemical system within a computer and allow this model to evolve according to a set of simple algebraic equations.These equations are derived from the physical laws that govern the electrochemical system under study.In effect, we carry out a simulation of the experiment with the aim to extract numeric representations of current, concentration profiles, potential transients, and so on [3].The results of digital simulations can then be compared to measured quantities, so that the physical parameters of the system can be fine-tuned in subsequent simulations, until the measured data are reproduced at a desired accuracy.
DS thus allows the determination of kinetic parameters, such as charge transfer rates and diffusion coefficients, etc. DS can yield important quantitative information even for complex electrochemical systems where the often-complicated set of partial differential equations, describing the transformation and movement of material, can be solved in closed form with difficulty or not at all [3].
In this paper I give a very brief and practice-oriented introduction to digital simulations and their use in electrochemistry.Instead of reviewing the topic of digital simulations in detailseveral good textbooks are dedicated to the subject and are recommended to the reader [3,4] -I choose a rather straightforward approach: that is, giving a step-by-step guide to the simulation of a simple electrochemical problem, that is known as Cottrell's experiment [5].
The paper is intended to be a starting point for undergraduate or graduate students, interested in digital simulations.Furthermore, the paper may also be found useful by university lecturers, as it contains a simpler than usual description and classification of different partial differential equation solving methods, based on Padé's approximants.

Experimental
Cottrell's experiment, first described in 1903 by F.G. Cottrell, is one of the oldest problems in electrochemistry that has a well-known analytical solution.The simulation of Cottrell's experiment is thus expedient, since the success of any simulation can directly be tested [6] by comparison to analytically obtained results.
In order to understand Cottrell's experiment, consider [7] a large planar electrode, of surface A, initially at rest, in contact with a semi-infinite layer of unstirred solution that contains some excess electrolyte and a small amount of a redox-active species R with a bulk concentration of c  .At the instant t = 0, the electrode potential is suddenly changed to a value at which the species R is oxidized to some product O in the one-electron reduction and the concentration of R, close to the electrode surface, is brought to essentially zero.In this experiment, the reaction rate (and thus, the current) is determined by the rate of transport at which the species R is resupplied at the electrode surface.That is, the partial differential equation governing the system is that describing Fick's law, Equation ( 2) is a second-order partial differential equation where c denotes the concentration and D the diffusion coefficient of the species R, while x and t denote the space and time variables, respectively.To obtain an analytical solution of Equation (2), we take two boundary conditions into consideration: i.) that the concentration of the species R is zero at the electrode surface at any t > 0 and ii.) that the concentration at infinite distance from the electrode surface never changes and remains equal to c  Taking the above boundary conditions into consideration leads to the solution: where the function erf(z) denotes Gauss's error function, defined by the integral: The current density j(t), passing the electrode surface at any time t > 0 is thus: Equation (3) makes it straightforward to calculate the concentration of species R at any given time t > 0 at a distance x measured from the electrode surface; while from Equation ( 4), the wellknown Cottrell equation, the current density can be calculated for any t > 0. These analytical formulae will be used below for testing the results of digital simulations.

Simulation approaches
LabVIEW, developed by National Instruments, is an easy-to-use, interactive, graphical programming language that helps users write sophisticated programs and applications in a short amount of time without needing a computer science degree [8].The popularity of LabVIEW continuously increases amongst students and scientific researchers, also in electrochemistry, mostly due to its versatility in measurement automation.However, LabVIEW also provides simple means for digital simulations and I therefore chose to present the basics of simulation algorithms in this programming language.Figure 1, discussed extensively below, shows the user interface of a program written for the simulation of Cottrell's experiment, as well as some details of the code.The main routine is that shown by the simulate.visnippet.
In most digital simulation schemes, we consider the electrolyte solution in terms of discrete and small (Δx wide) volume elements, as illustrated by Figure 2. When we apply this approach to the simulation of Cottrell's experiment in NI LabVIEW, we create a 1-dimensional array of concentrations, containing n entries, which we denote by c.At the start of the simulation process, all entries of the array c will be set equal to c  (the bulk concentration).The initialization of the c array in LabVIEW is shown by the create c array.vi snippet of Figure 1.
Each i entry of the array c (the index i can take values between 0 and n -1 in LabVIEW) corresponds to the concentration of the species R at a certain distance from the electrode surface, as shown by Figure 2.For accurate simulations, the Δx width of the control volumes must be small, thus to include a big enough space in the simulations, a sufficiently high number of the control volumes is required.(See Table 1 for numeric parameter values, used for the simulations presented in this paper.) When carrying out digital simulations, we iteratively modify the entries of the c array, as shown by the simulate.visnippet of Figure 1.Each iteration step represents a Δt time-step in which we replace the first element of c with zero (thereby realizing the boundary condition of the experiment) and then multiply the array (vector) c with a so-called stepper matrix, denoted here by S. As the iteration proceeds, we continuously "update" the c array, which at the end of the iteration will turn to be a discrete representation of the concentration profile, corresponding to the time t.The multiplication of the c array with the stepper matrix S in each iteration step mimics the role of transport in the simulation process, thus the construction of the stepper matrix has a deep impact on both the accuracy and the speed of the simulation.
To understand the role of the stepper matrix S in digital simulations, let us consider Figure 2, and assume that in each control volume shown, the concentration values are different.Let us then try to estimate the changes of concentration in the control volume marked by the index i during a small time-step Δt.Fick's first law of diffusion allows us to calculate the J left and J right fluxes of material from the left-hand (i -1) and the right-hand (i + 1) neighbours, directed to the i th control volume as: and where Δn left and Δn right denote the amount of substance arriving to the i th cell from its left and right neighbours, respectively.In each time-step Δt, the amount of substance that enters the i th cell is and taking into consideration the volume of the cell (AΔx), this results in a ( ) concentration change in the i th control volume.We note that the fraction on the right-hand side of Equation ( 6) approximates the second spatial derivative of the concentration profile, and that Equation ( 6) is in fact a discretized version, with respect to both space and time, of Equation ( 2), Fick's second law.
Equation (6) dictates that in a time-step of Δt duration, the change of concentration in the i th control volume depends on the concentrations in the i th , (i -1) th and (i + 1) th control volumes, as well as on the model parameters D, Δx and Δt.Equation ( 6) may be re-written in the form of a matrix-vector equation as: (7) by defining the so-called Laplacian matrix L, an n × n matrix as: ( ) where ( ) deg i is the number of neighbours of the i th control volume.For the simulation scheme shown in Figure 2, the Laplacian -a negative semi-definite matrix -takes the following form: Note that Equation ( 7) contains a finite difference (with respect to time) on its left-hand side, approximating a time derivative.For infinitesimally small time-steps, Equation ( 7) may be rewritten as: (10) which is a differential equation for the vector c.Assuming that the concentration vector () k c is known at any k th step of the simulation, the + ( 1) k c vector in the next step, Δt time later, could be calculated by solving Equation (10) in the form: where L w denotes the so-called weighted Laplacian, defined as: The LabVIEW code used for the calculation of the weighted Laplacian is shown by the weighted Laplace matrix.visnippet of Figure 1.
As can be seen in Equation (11), the "ideal" choice for the stepper matrix S we introduced before would be that the multiplication of the concentration vector c in each simulation step with this matrix could model the time evolution of the system, leaving the resolution of spatial discretization the only factor that limits the accuracy of calculations.
The problem is, however, that the exponential of the weighted Laplacian is not known and can only be approximated.One feasible way of approximating the exponential of a matrix is to truncate its Taylor series at a given order; another, more accurate way, is to use its Padé approximant [9].
It is this latter approach which I will follow here, as many digital simulation techniques -such as the explicit and implicit Euler, the Runge-Kutta and the Crank-Nicolson methods [4] -may be interpreted as techniques relying on the use of different Padé approximants.
Given a function f(x) and two integers m ≥ 0 and n ≥ 0, the Padé approximant of the function f(x) of order [m/n] is the rational function: which agrees with f(x) to the highest possible order, which amounts to: Equation ( 14) is a system of Equations, by solving which the parameters a j and b j of Equation ( 13) can be determined up to one degree of freedom; it is usually assumed that b 0 = 1.
In what follows, I will use P m,n (x) to note the Padé approximant of the exponential function exp(x).Exact expressions for P m,n (x) are tabulated in Table 2 for different values of m and n.Some Padé approximants to the exponential function are plotted, together with exp(x), in Figure 3 to show the "goodness" of these approximations.
Padé approximants, when applied to matrices instead of scalars, form the basis of many known simulation techniques, as will be discussed below.The P 1,0 (explicit Euler) method The explicit Euler method is the most straightforward one applied for the solution of partial differential equations.It is based on the P 1,0 (x) Padé approximation ( ) where I is the n × n identity matrix (see the identity matrix.visnippet of Figure 1).When plugged into Equation (11), it yields that The construction of the stepper matrix S, used for simulations based on the explicit Euler method, is shown in the create stepper matrix.visnippet of Figure 1.
Note that as 16) can also be directly obtained by rearranging Equation (10).The explicit Euler method can thus be interpreted as approximating the differential d / dt c in Equation (10) with the finite difference Δ / Δ .
Consequently, the explicit Euler method, a linear approximation, may only yield accurate results when a small enough time-step Δt is used.In-fact, when  2 Δ Δ 2 , t x D the explicit Euler method fails to converge and gives erroneous results, as will be discussed later.
The P 2,0 , P 3,0 and P 4,0 (Runge-Kutta) methods The overall error of the explicit Euler method can be decreased by taking higher order terms into consideration.This can be achieved by using the P 2,0 , P 3,0 and P 4,0 Padé approximants, which are in-fact Taylor series expansions of different orders to the exponential function.This forms the basis of the so-called 2 nd , 3 rd and 4 th order Runge-Kutta methods [10].While the explicit Euler method could be interpreted as a technique that takes the slope (but not the curvature) of the concentration vs. time dependence into account, Runge-Kutta methods aim to correct this error.
When applying Runge-Kutta methods, the stepper matrix S is constructed as The construction of the stepper matrix S, used for simulations based on the 4 th order Runge-Kutta method, is shown by the respective create stepper matrix.vi(detail) snippets of Figure 1.Note that -similarly to the Euler method discussed above -Runge-Kutta methods are still explicit; that is, in each iteration step, they calculate the state of a system at a later time from the state of the system at the current time.
The P 0,1 (implicit Euler) method The implicit (also called backward [4]) Euler method differs largely from explicit methods.The method is implicit, meaning that it finds a solution for the state of the system at a later time by solving an equation that involves both the current state of the system and the later (yet unknown) state.The stepper matrix used by the implicit Euler method, as also shown by the respective create stepper matrix.vi(detail) snippet of Figure 1, is ( ) Note that due to the implicit nature of the method, the creation of the stepper matrix involves matrix inversion.This clearly requires some extra computation, yet implicit methods can become very useful for the solution of stiff problems where explicit methods require the application of impractically small Δt values.As implicit methods are numerically stable, they can be applied with larger time-steps, and as usually the matrix inversion in Equation (20) has to be carried out only once (at the beginning of the iteration), implicit methods still require small computation times.

The P 1,1 (Crank-Nicolson) method
The method named after Crank and Nicolson [11] was developed for the numerical solution of partial differential equations, but it lies on the basis of the trapezium method of solving ordinary differential equations.The method is semi-implicit and is unconditionally stable; although when applied for stiff systems with large time-steps it may still be prone to spurious (decaying) oscillations.The stepper matrix used by the Crank-Nicolson method, as also shown by the respective create stepper matrix.vi(detail) snippet of Figure 1, is ( )

Comparison of the different simulation methods
As expected, applying Padé approximations of different order for the construction of the stepper matrix S yields solutions that also behave differently.The efficiency of the methods was tested by simulating concentration profiles corresponding to t = 5 s, using the parameter values of Table 1 and Δt values ranging between 1 ms and 5 s.Each simulated concentration profile was compared to the theoretically expected one, Equation (3), and the square-root of the mean squared deviation was plotted as a function of the used time-step in Figure 4.
Figure 4 clearly shows a major difference with respect to the stability of fully explicit (explicit Euler and Runge-Kutta) and implicit (implicit Euler and Crank-Nicolson) methods; at about  2 Δ Δ 2 , t x D all explicit methods begin to diverge (see also Figure 5).Compared to the explicit Euler, the 2 nd order Runge-Kutta method brings a significant improvement of the error; this improvement is less significant for the 3 rd and 4 th order Runge-Kutta methods.
The stability of the implicit Euler and the Crank-Nicolson methods is better (see also Figure 5), although the error of these methods is also significant at larger time-steps.

Summary
In this paper I attempted to describe numerical methods used for the digital simulation of a rather simple, however instructive electrochemical problem, Cotrell's experiment.The same techniques may also be used -by modifying the near-electrode boundary condition -for the simulation of more complicated electrochemical experiments, such as cyclic voltammetry.The described simulation strategies may also be extended to take into account homogeneous reactions [12], or effects related to ohmic drop [13,14].The further discussion of more complicated systems is, however, beyond the scope of this paper: my only intention here was to give a starting point for undergraduate or graduate students who decided to familiarize themselves with digital simulations.Accordingly, I attempted to describe basic numerical procedures in a simpler than usual manner, following a classification scheme based on Padé's approximants.Readers interested in the topic of digital simulation in more detail are referred to other sources [3,4].

Figure 1 .
Figure 1.The graphical user interface (top) of a simple LabVIEW program that simulates Cottrell's experiment.Snippets (code details) are shown in the blue boxes; see the text for explanation

Figure 3 .
Figure 3.Some Padé approximants P m,n (x) of the function exp(x) for the 4 th order method.(19)

Figure 4 .
Figure 4.The error of six different digital simulation methods, as a function of the applied timestep.Other simulation parameters are tabulated in Table 1

Figure 5 .
Figure 5. Concentration profiles simulated using six different methods for t = 7 s, using a relatively large time-step (Δt = 700 ms).Faded thick curves show theoretical profiles, as calculated from Equation (3).Forthe values of other simulation parameters, see Table1
 c bulk concentration 10 -6 mol cm -3 Δt time-step varied between 1 ms and 5 s Figure 2. Scheme for the digital simulation of Cottrell's experiment

Table 2 .
Some Padé approximants P m,n (x) of the function exp(x)