Optimisation of finite difference methods

In this chapter, the main obstacles to efficient and general finite difference simulations are discussed. These include the following:

Each of these are tackled in turn, with an explanation of their origin, significance and methods which may be used to overcome them.

3.1 Errors and Convergence

If iterative methods are used to solve the linear system of equations, there is a 'residual' error in their solution, which decreases with the number of iterations performed (thus the iterative solver converges on the true solution of the linear system). This error is related to the residuals, which are often available and are relatively inexpensive to compute, thus the norm of the residuals is often used as a measure of the error, and iterations continue until this falls below a specified tolerance. Rounding errors will also occur due to a computer's finite-precision arithmetic. The global error is therefore a composite of the truncation error, residual error in the linear system and rounding error. In practice the truncation error is the largest of these. Provided the difference approximation is consistent with the PDE it represents, as the mesh of points is made finer, the truncation error diminishes and the simulated concentration distribution, and the current evaluated from it, tends to the 'true' value. Moreover, the true solution will be reached regardless of the refinement path (the fraction by which the mesh in each dimension is refined by, in each step towards all mesh intervals being zero). This phenomenon is referred to as 'convergence' throughout this thesis. The number of nodes used must be high enough to ensure that the convergence error falls below an appropriate tolerance, often chosen to be 1% for experimental work and 0.1% for theoretical work1. This was achieved in this work by consecutively halving the mesh interval on one co-ordinate whilst keeping the others constant, and then repeating the procedure for the other (spatial and temporal) co-ordinates.

3.2 Extent of simulation space

The diffusion layer is the region where the concentration distribution is perturbed by the electrode - it represents how far, on average, the product diffuses away from the electrode (or the distance for the starting material to reach its bulk concentration value). Often the diffusion layer is confined to a very small region of the cell. This is especially important when simulating a microelectrode or hydrodynamic electrode such as the channel electrode (ChE) where the diffusion layer is compressed by convection. As an example, consider a typical ChE consisting of a 4mm square electrode in a channel 6mm wide by 0.4mm high (height=2h) containing ferrocene in acetonitrile (D=2.3x10-5cm2s-1) 2. At steady-state, the maximum extent of the (non-uniform) diffusion layer above the electrode is less than 30% of the channel height at a flow rate of 0.1cm3s-1. Nodes outside the diffusion layer provide no useful information, so a ChE simulation may be restricted to a height 2fh above the electrode surface, where f is 0.3 for the example system. Since f depends on electrode size, flow rate etc. it must be optimised for each simulation.

Figure 3.1: Steady-state concentration profile at a channel flow cell (simulated using the BI method).

Another spatial consideration in microelectrode simulations is the extent of axial (or radial) diffusion. In the channel, material diffuses from both upstream and downstream to the electrode surface and some space must be allowed in the simulation to account for this (the relative importance of upstream and downstream contributions are examined in Chapter 6). The approach adopted in this thesis is to allow lu electrode lengths upstream and ld downstream in order to accommodate this phenomenon. As the flow rate and electrode size decrease, the degree of axial diffusion increases and lu and ld must be increased accordingly.

3.3 Boundary Singularities

Boundary singularities are caused by discontinuities in boundary conditions. There are two classic examples in electrochemistry. A spatial discontinuity occurs at the edge of a microdisc electrode, where it meets the insulating surround. The diffusion layer above the electrode is hemispherical, so at this edge the diffusion layer thickness tends to zero and the flux is infinite. A temporal singularity occurs in potential step chronoamperometry. At the moment the electrode is 'switched on', the flux is infinite due to the zero diffusion layer thickness. Boundary singularities dramatically impede the rate of convergence of numerical simulations. There are three main ways of tackling a boundary singularity:

3.4 Efficient distribution of nodes

With a simple rectangular mesh, the number of nodes required for a converged solution can be relatively high (meshes of 1000x1000 nodes have been used for BI channel electrode simulations9). Ideally the distribution of the nodes should mimic the concentration distribution of species in the cell. Steep concentration gradients occur in the region above the electrode, so closely-packed nodes near the electrode surface will improve the simulation accuracy.

There are several ways of improving the distribution of nodes. These methods may be subdivided into two categories. The first are known as expanding grids - a mathematical function controlling the node distribution is 'tuned' via one or more adjustable parameters to mimic the concentration profile. The second category is known as a conformal mapping - the co-ordinates are transformed into a space where the concentration profile is smoother.

Whilst such mappings greatly improve the efficiency of simulations, the system is usually more ill-conditioned in the conformal space due to the large range of coefficient magnitudes. Feldberg et al.10 found that whilst the Laasonen implicit method is unconditionally stable on a uniform grid, this is not the case if the co-ordinate is transformed via an exponential mapping.

Figure 3.2: Finite differences formulated over unequal intervals.

3.4.1 Expanding grids

The most intuitive way of distributing the nodes in a non-uniform fashion is to formulate the finite differences over unequal intervals, illustrated in Figure 3.2. Thus:




A simple switch in grid density (from a coarse to fine rectangular grid) was used by Bidwell et al.11 to simulate the concentration distributions arising from double channel electrodes. This technique introduces an error in a manner analogous to upwind differencing since both the first and second derivatives are displaced from the point (j,k).

Figure 3.3: Interpolative expanding grid used for wall-jet BI simulations.

A space-marching method, such as the BI method, used to solve pseudo two-dimensional problems, may modify the extent of its rectangular mesh (in z) as it marches along the electrode (in r), as shown in Figure 3.3. The concentration values must be interpolated from the previous z positions to provide values at the new z positions (all at the old r position). From this the y-vector at the new r position may be calculated. This method was used by Compton et al.12 to simulate the wall-jet electrode (WJE) where the highly non-uniform diffusion layer makes a uniform z grid very inefficient.

A more elegant way of achieving an unequal spacing is through the use of a co-ordinate transformation. The simulation is conducted on a uniform mesh in the transformed co-ordinate, which corresponds to an expanding mesh in 'real' space, thus avoiding displacement errors and instability. Similarly the WJE problem may be solved in a transformed z co-ordinate thus eliminating interpolation errors. This approach was adopted by Ball et al.13 to simulate anodic stripping voltammetry, and Laevers et al.14 to simulate chronoamperometry, at the WJE.

A general transformed co-ordinate:

x = x(x)(3.3)

must be differentiated readily if it is to be used in a co-ordinate transformation:

and (3.4)

These are then substituted into the mass transport equation, but the derivatives of x are initially expressions in terms of x, not x. As Britz15 points out, if the function x(x) is readily inverted to give x(x), then the mass transport may be expressed purely in terms of the transformed co-ordinate, x.

Otherwise values of x must be calculated (and stored) at each value of x across the uniform grid in x. Since


has already been found, this may be achieved by solving the differential equation:


There are several standard methods for this integration such as Runge-Kutta and Adams' methods16 which are available as subroutines in most mathematical libraries (D02BJF and D02CJF in the NAG FORTRAN library). These require a boundary value (i.e. a value of x where x is known) from which to start the integration. This is not usually a problem since the function chosen has a known finite value at the origin. This method was used by Compton et al.17 for a Hale space transformation to simulate an RDE, using a 3rd order Runga-Kutta integration. In this work an Adams method is used via the NAG subroutine D02CJF (since these are generally more efficient than the Runge-Kutta methods18).

Figure 3.4: 'Exponential' expanding grid functions.

A number of 'exponential' functions, shown in Figure 3.4, have been used to expand the co-ordinate normal to the electrode (given here in normalised form) including:

and (3.7)

used by Joslin and Pletcher19,


originally proposed by Feldberg20 and used in DigiSim(tm). Other possible functions include:


based on Pastore's work (see below) and:


used in this thesis.

In addition, a few functions have been used to transform the axial co-ordinate which expand outwards (usually to ±infinity) from a central point. Pastore et al.21 used a hyperbolic function:


to expand the axial co-ordinate in the channel flow cell. The efficiency and applicability of this function is assessed in Chapter 6. Taylor et al.4 used a Fermi-Dirac (sigmoidal) function:


to concentrate nodes near the edge of a microdisc electrode. In Chapter 8, this is used for a similar purpose in a wall-jet simulation. In Chapter 6, a pair of these functions are used for a channel flow-cell co-ordinate transformation.

3.4.2 Conformal mappings

Conformal mappings are virtually identical to expanding grids based on a transformed co-ordinate. The only difference is that in the former, the co-ordinate transformation does not involve an adjustable parameter.

In the 1970s Hale22 used an integral transformation for efficient finite difference simulations of the Rotating Disc Electrode. For a simple electron transfer the concentration profiles are linear in this transformed space (in the 'classical' limit of high Schmidt numbers - see Chapter 7). Later Compton et al.17 conducted explicit finite difference simulations in transformed space. In Chapter 7 this transformation is used, with a more accurate convection approximation than Hale or Compton used, in the general simulator (described in Chapter 5) to simulate kinetics at an RDE.

Michael et al.23 devised a co-ordinate transformation for the microdisc electrode which eliminates the edge singularity. This was later refined into a closed-space form by both Amatore and Fosset24 (based on work by Newman25) who used the Hopscotch algorithm to simulate steady-state and transient processes and Verbrugge and Baker26 who used an explicit method to simulate chronoamperometry. For a simple electron transfer, the concentration profile is invariant with the angular co-ordinate under both transformations. Under Amatore's transformation the concentration profile is also linear in the radial co-ordinate. In Chapter 9 the efficiency of two transformations is compared for simulating microdisc steady-state problems, and the better conditioned of the two (Amatore's) is used in the general simulator for a range of microdisc simulations.

Deakin et al.27 also developed a similar mapping for a microband electrode, based on a Schwartz-Christoffel transformation. This was later extended to a double band arrangement and a triple band array28.

3.4.3 Transformed time

It is also possible to transform the time co-ordinate15. This is particularly useful in chronoamperometry where Cottrell-like behaviour occurs. A simulation with an 'exponentially' expanding time co-ordinate may considerably improve the efficiency of such a simulation. Lavagnini et al.29 used an exponentially expanding time grid with the Hopscotch algorithm to simulate voltammetric curves at a microdisc electrode. Palys et al.30 also used an exponentially expanding time grid for microelectrode simulations. Feldberg and Goldstein31 examined the behaviour of the Laasonen method with an exponentially expanding time grid together with the Richtmyer modification.

A Cottrellian time transformation:

hence (3.13)

has been used to achieve significant savings in CPU time: Prieto32 simulated chronoamperometry at a channel microband electrode and reported at least a 20-fold efficiency improvement; Ball and Compton33 simulated square-wave voltammetry at hydrodynamic mercury electrodes, achieving a 50-fold efficiency improvement over a non-transformed co-ordinate.

3.5 Reaction Layer

When an electron-transfer is followed by a homogeneous chemical reaction, the material escaping from the electrode surface reacts as it diffuses away. Its distribution is therefore mainly confined close to the electrode surface in what is termed a 'reaction layer'. For a fast homogeneous process, the reaction layer may span only a tiny fraction of the diffusion layer. This can lead to poor convergence in simulating the reacting species, since the nodes must span the diffusion layer and satisfactorily describe its curvature in order that the other species may be simulated accurately. Expanding grids and conformal maps improve this situation over a simple Cartesian mesh since they generally concentrate points in the vicinity of the electrode surface. However even in the 'conformal space' the fundamental problem still remains, and can only be addressed by increasing the number of nodes. With an expanding grid, the expansion parameter may be 'tuned' to result in a grid which compromises between the distributions of species. However this is laborious, requiring trial and error at each rate constant, and does not substantially improve the rate of convergence over an expansion parameter based on the diffusion layer. The compromising value of the expansion parameter will be closer to the optimal value for the diffusion layer than the reaction layer, since a grid which expands too rapidly gives far worse results than a grid which expands too slowly, which is reflected in the asymmetric plots of current vs. expansion parameter one obtains when optimising.

3.6 Coupled heterogeneous processes

Reversible chemical reactions lead to coupling of species. Mathematically this means that the material balance equations of these species must be solved simultaneously. Heinze and Störzbach34 claimed to simulate the coupled system using ADI, though did not provide any details. Taylor et al.4 used a modified ADI method to couple species at the boundary (this involved including the boundary nodes in the simulation). A matrix method may be used to solve the entire system in an analogous manner to that discussed below for homogeneous kinetics. Generally, however, if 'standard' finite difference solvers are to be used which only allow a fixed 'stencil' (sparsity pattern), additional methods must be used to allow for coupling through an electrode surface boundary condition.

Compton et al.9 tackled this problem using an unaccelerated successive substitution method to de-couple the species. An initial guess was provided for the surface concentration of one species. This was used to calculate the concentration profile of the other species and vice versa until the system gradually converged to within a specified tolerance. Since simulation of only a single couple was required, the (rapid) BI method could be used to perform simulations. Even so, the CPU time required for so many iterations is in the order of hours, compared to a couple of seconds for a typical single-species BI simulation.

3.6.1 Back-to-Back Grid (L.Mitoseriu, 1995)35

This is an elegant way of overcoming boundary coupling. This works by placing the concentration profiles of the two species in a couple 'back-to-back', allowing the resulting matrix equation to be solved in a single pass.

Figure 3.5: Back-to-back finite difference grid for a channel electrode.

Figure 3.5 shows the back-to-back finite difference grid for a channel electrode, with the x co-ordinate node index, k, running from 1 to NGX in the x direction as usual. The electrode region is denoted by lu.kE+1 to (lu+1).kE. On the back-to-back grid the y co-ordinate index, j runs from 1 to NJ for species A, and NJ+1 to 2NJ for species B. If we denote a general node in Figure 3.5 as gj,k then the concentration profiles of A and B map onto the grid as:




In particular, on the back-to-back grid, J=1 and 2NJ correspond to the node below the wall of the channel opposite the electrode for species A and B, respectively. Hence the point g1,k corresponds to aNJ,k, and g2NJ,k corresponds to bNJ,k. Likewise, on the back-to-back grid j=NJ and j=NJ+1 correspond to the node above the surface of the electrode for species A and B respectively. Thus gNj,k corresponds to Au1,k, and gNJ+1,k corresponds to Bu1,k.

The mass transport equation of species A at one node above the electrode surface is:


Since the simulation grid does not include the electrode surface, the finite difference Butler-Volmer equation (see section 2.5.1):


may be substituted for the electrode surface concentration:


This gives rise to a finite difference equation in which the concentrations of species A and B are interdependent at the node above the electrode surface. If equation (3.18) is transformed onto the back-to-back grid:


the concentrations at one node above the electrode surface are adjacent in the new finite difference grid and therefore may be related to each other within the finite difference stencil. Equation (3.19) simplifies to:


Similarly for species B:


The quasi-reversible electrode boundary conditions may be implemented by combination of the bands in the coefficient matrix. If the three bands in the tridiagonal matrix are denoted a, b, c and the RHS of the linear system as d (see section, for species A the electrode boundary condition may be implemented by:


Note that d must be used to modify c before d itself is modified. Similarly for species B:


3.7 Coupled homogeneous processes

In order to simulate a homogeneously coupled system implicitly, an analogue of the back-to-back grid may be constructed using a 3 dimensional finite difference method. The coupled species may be formulated in adjacent planes, and homogeneous kinetics represented as transport of material between the planes. This is attractive in that a 'standard' finite difference solver may be used such as 3D SIP or ADI (introduced in sections and 2.1.6, respectively). There are 2 major drawbacks:

  • As with the back-to-back grid, the mechanistic complexity is limited by the number of species which may be placed adjacent to each other.
  • Also in common with the back-to-back grid, the resulting system may be ill-conditioned.

    A simulation of a DISP1 process at a channel microband electrode was attempted using 3D SIP. It was not possible to obtain a stable solution, which implies the system was ill conditioned. This may be attributable to the off-diagonal dominance arising from the kinetic terms.

    Alternatively the finite difference equations of all the species may be expressed as a single matrix equation. For each kinetic relationship, an extra band (or two bands for a reversible reaction) is added to the basic matrix defined by the mass transport. A matrix method (capable of accommodating an arbitrary sparsity pattern) may be used to solve such a system. This forms the basis for the general simulation strategy described in Chapter 5.

    In a time-marching or space-marching method, the kinetics may be formulated explicitly (in terms of concentrations at the 'old' time or 'previous' space region). This seems very attractive, since it de-couples the species, allowing them to be solved sequentially using a standard finite difference method. The major drawback is that this approximation breaks down at high rate constants. As the rate constant increases the size of the step must be reduced in the co-ordinate across which the explicit approximation is made. This leads to increasing CPU time, which can become ridiculous at very large rate constants.

    Figure 3.6: BI simulation of an ECE process at a wall-jet electrode.

    Figure 3.6 shows the CPU time vs. the dimensionless rate constant for a BI simulation of an ECE process at a wall-jet electrode, for either an implicit or explicit formulation of kinetics (assuming a minimum of 500 nodes are required for a converged solution in the x co-ordinate). The increasing CPU time for the explicit simulation at high rate constant arises from the increasingly large NGX values required to keep the simulation stable.

    3.8 Neumann boundaries

    Often we need to specify a boundary condition as a derivative, usually a flux. In the previous chapter it was shown how this could be achieved in the finite difference scheme by representing the derivative as a simple difference of two values. It is possible to use a Taylor (or polynomial) expansion15 about the boundary to give a higher order backward difference formula employing more than two points to calculate this derivative and thus obtain a more accurate estimate of its value. This is analogous to the Richtmyer modification used in improving the resolution of time differences discussed in Chapter 2, and the coefficients are the same.

    Higher order difference formulae may be used to specify both the flux to the electrode surface (for example in formulating heterogeneous kinetics using the Butler-Volmer equations) and the 'no flux' boundaries at insulating walls etc.1 Three points have often been used in the literature since they may be incorporated into the standard finite-difference 'stencil'36,37.

    Higher order series have been used to evaluate the flux. Since this takes place after the finite difference simulation, it is merely a matter of using more concentration values to obtain a better estimate of the current15,38. Often the value predicted by a first-order expression is an underestimate of the flux, whilst that predicted by a higher-order expression is an overestimate (or vice versa).

    Figure 3.7: Convergence of a BI channel simulation using a 2-point or 6-point Taylor series to evaluate the flux.

    This provides an error-bound on the current, as shown in Figure 3.7 - a convergence plot for a BI simulation of steady-state electrolysis in a channel flow cell. In this case, the ratio of the error in the 2-point and 6-point currents is approximately constant. If this ratio can be estimated, the weighted average of the 2-point and 6-point functions provides an extremely rapidly convergent formula, as shown in Figure 3.8. In the case of this plot, the ratio was calculated from the converged solution, giving the error ratio:


    Clearly in practice, calculating the ratio from the converged solution defeats the point, but the method does suggest an SOR-type approach where a 'roughly optimal' linear combination could be used in many cases.

    Figure 3.8: Effect of Taylor series or an expanding grid on the rate of convergence of a BI channel flow simulation. The line marked 'Taylor ratio' uses the optimal weighted average of the 2-point and 6-point Taylor series solutions on a Cartesian mesh.

    In order to use a series of more than 3 points to formulate a derivative boundary condition, the solver must be capable of handling sparsity patterns other than the standard finite-difference stencil. In order to achieve this whilst keeping the simulation fully implicit, either a full elimination must be performed, or a sparse iterative method can be used. Bieniasz38 and Gavaghan1 both argue that there is little to be gained from using anything higher than second order boundary conditions since the 5-point central finite difference approximation is itself only second-order accurate.

    3.9 Current Integration

    The flux at a uniformly accessible electrode is constant across the electrode surface, and therefore the current may be calculated from the concentration gradient at the electrode surface. For a non-uniformly accessible electrode this is not the case, as the flux, by definition, varies as a function of at least one co-ordinate across the electrode, across which the current must be integrated. Thus for a band electrode the total flux is given by:


    and for a disc electrode:


    Since the concentration is represented as discrete points along this co-ordinate a form of quadrature39 is necessary to approximate the integral from these discrete values.

    The simplest technique is rectangular integration where:


    Note that since we are summing rectangular 'bars' at each point, the electrode strictly starts and finishes between nodes, as shown in Figure 3.9(a) (i.e Dx = xe/KE). The alternative is to define the electrode starting exactly on a node (Dx = xe/(KE-1)) in which case more accurate trapezoidal integration, shown in Figure 3.9(b), may be used:


    Higher order methods such as Simpson's (3rd order) rule39:


    should give better even better results than trapezoidal integration. Heinze37 used a method for integrating the flux at a microdisc electrode, specially adapted to take account of the edge singularity. Gavaghan1 compared this with trapezoidal and Simpson integration for a microdisc electrode and found it offered little improvement. Cubic spline fitting has also been used by Lavagnini40 to achieve more accurate integration.

    Figure 3.9: Integration of the concentration at a channel microband electrode. The red and blue lines represent the 'true' concentration (simulated using 271 nodes along the electrode), the points are from a coarse-mesh simulation (using 8 nodes along the electrode). The integration areas are shown as green and brown lines for rectangular and trapezoidal integration, respectively.

    Figure 3.10: The current converges more rapidly with the number of nodes in the x co-ordinate due to error cancelling from rectangular integration.

    Figure 3.9 suggests that any integration based on the electrode on a node will consistently under-estimate the current (for an electrolysis). From Figure 3.10 it is evident that by defining the electrode to start between nodes and using rectangular integration, significant error cancelling occurs. However, when integrating across a transformed co-ordinate, the integrals in each space are related by:


    Consider a trapezoidal integration:


    as one would expect this is equivalent to working out the x values at each x value:

    where (Dx)k = x(xk) - x(xk-1).(3.32)

    If the rectangular method were used, since the electrode would lie half-way between points in transformed space, the electrode would change size (in real space) as the number of nodes was altered! Therefore the error cancellation method is best suited to simulations on non-transformed x-meshes.


    1 D.J. Gavaghan, J. Electroanal. Chem., 420, (1997), 147.
    2 P. Sharp, Electrochimica Acta, 28, (1983), 301.
    3 R. Ferrigno, P.F. Brevet, H.H. Girault, Electrochimica Acta., 42, (1997), 1895.
    4 G. Taylor, H.H. Girault, J. McAleer, J. Electroanal. Chem., 293, (1990), 19.
    5 S.W. Feldberg and C.I Goldstein, J. Electroanal. Chem., 397, (1995), 1.
    6 J. Crank, R.M. Furzeland, J. Inst. Math. Applic., 20, (1977), 239.
    7 D.J. Gavaghan, J.S. Rollett, J. Electroanal. Chem., 295, (1989), 1.
    8 J. Galceran, D.J. Gavaghan, J.S. Rollett, J. Electroanal. Chem., 394, (1995), 17.
    9 R.G. Compton, M.B.G. Pilkington, G.M. Stearn, J.Chem.Soc. Faraday Trans. 1, 84, (1988), 2155.
    10 S.W. Feldberg, C.I. Goldstein, M. Rudolph, J. Electroanal. Chem., 413, (1996), 25.
    11 M.J. Bidwell, J.A. Alden, R.G. Compton, Electroanalysis, 9, (1997), 383.
    12 R.G. Compton, C.R. Greaves, A.M. Waller, J. Appl. Electrochem., 20, (1990), 575.
    13 J. Ball, R.G. Compton, C.M.A. Brett, J. Phys. Chem. B., 102, (1998), 162.
    14 P. Laevers, A. Hubin, H. Terryn, J. Vereecken, J. Appl. Electrochem., 25, (1995), 1023.
    15 D. Britz, Digital Simulation in Electrochemistry, 2nd Edition, Springer-Verlag, (1988).
    16 H.R. Schwartz, Numerical Analysis: A Comprehensive Introduction, Wiley, (1989), pp 384-413.
    17 R.G. Compton, M.E. Laing, D. Mason, R.J. Northing and P.R. Unwin, Proc. R. Soc. Lond. A, 418, (1988), 113.
    18 K. Atkinson, Elementary Numerical Analysis, 2nd Edition, Wiley, (1993), p339.
    19 T. Joslin, D. Pletcher, J. Electroanal. Chem., 49, (1974), 171.
    20 S.W. Feldberg, J. Electroanal. Chem., 127, (1981), 1.
    21 L. Pastore, F. Magno, C.A. Amatore, J. Electroanal. Chem., 301, (1991), 1.
    22 J.M. Hale, J. Electroanal. Chem., 6, (1963), 187.
    23 A.C. Micheal, R.M. Wightman, C.A. Amatore, J. Electroanal. Chem., 26, (1989), 33.
    24 C.A. Amatore, B. Fosset, J. Electroanal. Chem., 328, (1992), 21.
    25 J. Newman, J. Electrochem. Soc., 113, (1966), 501.
    26 M.W. Verbrugge, D.R. Baker, J. Phys. Chem., 96, (1992), 4572.
    27 M.R. Deakin, R.M. Wightman, C.A. Amatore, J. Electroanal. Chem., 215, (1986), 49.
    28 B. Fosset, C.A. Amatore, J. Bartlet, R.M. Wightman, Anal. Chem., 63, (1991), 303 and 1403.
    29 I. Lavagnini, P. Pastore, F. Magno, C.A. Amatore, J. Electroanal. Chem., 316, (1991), 37.
    30 M.J. Palys, Z. Stojek, M. Bos, W.E. Vanderlinden, J. Electroanal. Chem., 383, (1995), 105.
    31 S.W. Feldberg, C.I. Goldstein, J. Electroanal. Chem., 397, (1995), 1.
    32 F. Prieto, unpublished work. email dapena@cica.es.
    33 J. Ball, R.G. Compton, J. Phys. Chem. B., 102, (1998), 3967.
    34 J. Heinze, M. Störzbach, Ber. Bunsenges. Phys. Chem., 90, (1986), 1043.
    35 R.A.W. Dryfe, DPhil Thesis, Oxford University, (1995), p82.
    36 J.E.B. Randles, J. Chem. Soc. Faraday. Trans., 44, (1948), 322.
    37 J. Heinze, J. Electroanal. Chem., 124, (1981), 73.
    38 L.K. Bieniasz, O. Osterby, D. Britz, Computers. Chem., 19, (1995), 351.
    39 W.H. Press, W.T. Vetterling, S.A. Teukolsky, B.P. Flannery, Numerical Recipes in C, Cambridge, (1992).
    40 I. Lavagnini, P. Pastore, F. Magno, C.A. Amatore, J. Electroanal. Chem., 316, (1991), 37.