A General Method for Electrochemical Simulations. 1. Formulation of

J. Phys. Chem. B 1997, 101, 8941-8954

8941

A General Method for Electrochemical Simulations. 1. Formulation of the Strategy for Two-Dimensional Simulations John A. Alden and Richard G. Compton* Physical and Theoretical Chemistry Laboratory, Oxford UniVersity, South Parks Road, Oxford OX1 3QZ, United Kingdom ReceiVed: June 13, 1997; In Final Form: August 8, 1997X

A strategy for a general electrochemical simulator is presented based on formulating the finite difference form of coupled mass transport and kinetic equations as a single sparse matrix. A treatment of a general mechanism consisting of first- and second-order homogeneous steps and quasi-reversible heterogeneous couples is developed for both steady-state and transient conditions. This is illustrated for two-dimensional simulations with potential applications including microdisk and channel microband electrodes, though the method may easily be applied to higher or lower dimensions and could include terms describing mass transport due to migration. Nonlinearities due to second-order kinetics are treated using Newton’s method with an analytically derived Jacobian. ILU and MILU preconditioned Krylov subspace methods (CGS, BICGSTAB(l ), and RGMRES) are used to solve the resulting nonsymmetric sparse linear system. The efficiency of these solvers is compared to “brute force” LU factorization, SIP, and multigrid methods for a typical steady-state channel microband problem.

1. Introduction To predict the current flowing due to particular electrochemical mechanism at a given electrode geometry as a function of time or potential, the combined mass transport and kinetic equations must be solved for most of the species involved in the reaction. For many macroelectrode geometries, the mass transport can be described by one space dimension and may be solved analytically for simple mechanisms, or by relatively simple time-marching finite difference numerical methods1-3 for general mechanisms.4 Such simulations are the basis for commercial packages such as DigiSim5 and ELSIM.6 The mass transport to almost all microelectrodes is inherently two-dimensional and often tends to a steady-state, so more complex numerical modeling techniques are required. Electrochemical finite difference simulations have been conducted using a range of methods: explicit methods7 were rapidly superseded by the more efficient quasi-explicit methods such as the hopscotch8 method, although Feldberg9 found this method to be inaccurate at early times. Other semi-implicit relaxation methods which are more efficient, accurate, and stable than hopscotch include ADI10 and successive over-relaxation (SOR), which Gavaghan11 has used to simulate the steady-state current at a microdisk electrode. Implicit methods such as SIP12,13 and multigrid14,15 methods are more even more efficient and stable2 and allow direct solution of steady-state problems. Linear matrix methods may also be used to solve the simultaneous equations arising from the application of the finite difference method. Friesner, Bard, and co-workers16,17 showed how a one-dimensional diffusion equation may be written as a matrix equation and, because the matrix in this case was symmetric, solved using a Lanczos (early Krylov subspace) solver. Later they applied this method to solve the chronoamperometric response of a microdisk in a scanning electrochemical microscope.18 Since the mass transport occurs in two space dimensions, they separated the equation into ordinary differential equations in each coordinate and solved these using a Krylov subspace method.19 Georgiadio and Alkire20 used a X

Abstract published in AdVance ACS Abstracts, September 15, 1997.

S1089-5647(97)01936-6 CCC: $14.00

quasi-Newton method with a frontal solver to solve the convection, diffusion, and migration for the anisotropic chemical pattern etching of copper foil. Mao and White21 used a block matrix LU decomposition method to solve the one-dimensional equations for steady-state pseudo-two-dimensional boundary value problems. Balslev and Britz22 used a Newton-Raphson iterative scheme with Gauss-Jordan elimination to solve the one-dimensional transport with nonlinear kinetics at a rotating disk electrode. Britz23 also showed how kinetics and onedimensional transport to an electrode or the two-dimensional transient diffusion of a single species to a microdisk electrode can be represented as a matrix equation which may be solved using a direct LU factorization. As intimated by Britz,23 a matrix-based linear solver may be used to simulate diverse electrochemical mechanisms at a wide range of electrode geometries fully implicitly because it is not constrained by a finite difference “stencil” (that is, any sparsity pattern may be accommodated). This is the major drawback with methods such as hopscotch, ADI, SOR, SIP, and multigrid: they can only be applied to mechanisms where species may be solved sequentially by treating the term coupling the species as irreversible24 unless one of the following simplifying techniques is implemented. Kinetics may be treated explicitly25 (approximated to a function of the concentrations at the old time) in a transient simulation,26 but increasingly smaller time steps are required to simulate larger rate constants, leading to ridiculous simulation times. An explicit approximation may be iterated to self-consistency at steady-state,27 but this also leads to a huge multiplication of the CPU time. A back-to-back grid28 may be used to facilitate coupling between species, although this can result in a badly conditioned matrix, causing solver instability. A matrix method not only removes limitations from mechanistic coupling but also allows Taylor series of more than three points to be used to represent fluxes when formulating boundary conditions if necessary. A major drawback with Britz’s direct LU “brute force” approach is that all the zero elements in the matrix must be stored, imposing severe memory limitations on the size of the simulation. As noted by Gavaghan,11 two-dimensional simula© 1997 American Chemical Society

8942 J. Phys. Chem. B, Vol. 101, No. 44, 1997

Alden and Compton

tions often contain boundary singularities at electrode edges due to a change in the boundary condition between the electrode surface and its insulating surround. The presence of one or more such singularities necessitates a fine grid mesh to obtain satisfactory convergence. Although this may be overcome for symmetric problems such as the diffusion-only response of a microdisk electrode via conformal mappings such as those of Amatore and Fosset29 and Verbrugge and Baker,30 hydrodynamic electrodes such as the channel microband give rise to inherently asymmetric concentration distributions which are not conformally mapped so easily.31 Species with markedly different diffusion coefficients and the presence of homogeneous kinetics also cause conformal mappings to become less efficient. Here we introduce some sparse matrix methods: a range of solvers for sparse nonsymmetric linear systems. Unlike the direct matrix methods, these methods only require storage of, and calculation using, the nonzero elements in the coefficient matrix. Despite this, they retain the flexibility of the direct matrix methods: they may be used to solve randomly structured, randomly valued sparse matrix equations. These solvers (known as preconditioned Krylov subspace methods) have been applied successfully to finite difference and finite element simulations in quantum dynamics32,34 and other fields.35,36 Conveniently, these methods will be available in the next (Mk 18) release of the NAG FORTRAN library so they may be easily used by experimental electrochemists. We show how an electrochemical simulation for any general mechanism at any electrode geometry may be represented as a sparse matrix equation. This is further generalized to include either steady-state or transient experiments such as chronoamperometry, linear sweep voltammetry, and cyclic voltammetry. We combine preconditioned Krylov subspace methods with Newton’s method,24,37 (using an analytically derived Jacobian) for solving the nonlinear sets of equations arising from second-order mechanisms. The result is an electrochemical simulator which can be applied to a general mechanism at an electrode of arbitrary geometry. The Krylov subspace iterative solvers CGS,38 BICGSTAB(l ),39,40 and RGMRES41 have been shown42 to be highly flexible and efficient especially when used with an incomplete LU factorization43 as a preconditioner. The efficiency of these methods is compared with established simulation methods such as the multigrid method, SIP, and direct LU factorization of the linear system, for a typical electrochemical simulation.

Figure 1. Radial coordinate system used to describe diffusion to a microdisk electrode.

Figure 2. Channel microband electrode. The aspect ratio of the cell is such that the concentration distribution is effectively uniform across the electrode in the z-coordinate.

ficient. Equation 2 may be cast into a central finite difference form as

-(λrd - λrs)aj,k-1 - λzaj-1,k + (2λrd + 2λz)aj,k - λzaj+1,k - (λrd + λrs)aj,k+1 ) 0 (3) where

λz )

A+ehB

(1)

The current at a microdisk electrode, in a stationary solution where migration is eliminated, reaches a steady-state due to the convergent diffusion field,44 as depicted in Figure 1. At steadystate, the mass transport to the microdisk electrode is given by45

∂[A] ∂2[A] ∂2[A] DA ∂[A] ) DA 2 + DA 2 + )0 ∂t r ∂r ∂z ∂r

(2)

where the coordinates are as shown in Figure 1, [A] denotes the concentration of species A, and DA is its diffusion coef-

(∆z)

, λrd )

2

DA (∆r)

, and λrs )

2

DA 2∆r

(4)

and aj,k denotes the concentration of species A at a finite difference grid node with indices j ) 1...NGY in the z coordinate and k ) 1...NGX in the r coordinate. ∆r and ∆z represent the distance between adjacent finite difference nodes in the r and z coordinate, respectively. The current at a channel microband electrode (depicted in Figure 2) also reaches a steady-state due to the fresh supply of material as a result of convection.46 At steady-state, the mass transport to the electrode is given by47

∂[A] ∂2[A] ∂2[A] ∂[A] ) DA 2 + DA 2 - Vx )0 ∂t ∂x ∂x ∂y

2. Theory In this section we formulate a method with which the finite difference discretization of an electrochemical mechanism may be represented as a banded sparse matrix. We will develop the problem at steady-state and then show how it may be extended to time-dependent experiments such as chronoamperometry and cyclic voltammetry. The simplest case is considered first: a single electroactive species (A) undergoing transport-limited electrolysis:

DA

(5)

where the coordinates are defined in Figure 2 and Vx is the velocity of the solution in the x coordinate. This may be cast into a similar finite difference equation:

-(λx + λV)aj,k-1 - λyaj-1,k + (2λx + 2λy)aj,k - λyaj+1,k - (λx - λV)aj,k+1 ) 0 (6) where

λy )

DA 2

(∆y)

, λx )

DA (∆x)

, and λV )

2

Vx 2∆x

(7)

and the indices j and k and internodal distances ∆x and ∆y refer to the x and y coordinates, respectively. In both cases the equations are of an implicit five-point form, which may be expressed generally (using the coefficient labels

Electrochemical Simulations

J. Phys. Chem. B, Vol. 101, No. 44, 1997 8943

used in the NAG library implementation of the SIP finite difference method13) as

aaj,kaj,k-1 + baj,kaj-1,k + caj,kaj,k + daj,kaj+1,k + eaj,kaj,k+1 ) qaj,k (8) where qaj,k is zero and

aaj,k ) -(λrd - λrs) baj,k ) -λz caj,k ) (2λrd + 2λz)

(9)

daj,k ) -λz

mechanism with a reversible chemical step occurring at a microdisk electrode are

eaj,k ) -(λrd + λrs)

∂[A] ∂2[A] ∂2[A] DA ∂[A] ) DA 2 + DA 2 + )0 ∂t r ∂r ∂z ∂r

for the microdisk, and

aaj,k ) -(λx + λV)

∂[B] ∂2[B] ∂[B] DB ∂[B] ) DB 2 + DB 2 + ∂t r ∂r ∂z ∂r - k1[B] + k-1[C] ) 0 (14)

baj,k ) -λy caj,k ) (2λx + 2λy)

(10)

daj,k ) -λy eaj,k ) -(λx - λV)

(11)

The vector u contains the concentrations at each node on the two-dimensional grid, in the order given by a storage map equation. The row or column order is arbitrary (either j or k could be in the inner loop) but must be consistent. We therefore use the following storage mapping throughout:

ui ) aj,k where i ) k + (j - 1)NGX

∂[C] ∂2[C] ∂2[C] DC ∂[C] ) DC 2 + DC 2 + ∂t r ∂r ∂z ∂r + k1[B] - k-1[C] ) 0 These give rise to the finite difference equations:

for the channel microband. (The superscript a in the notation for q and a-e refers to species A, for clarity in the expressions below for multi species mechanisms). Equation 8 may be solved using stencil-based methods such as the SIP12 and multigrid14 solvers. This set of finite difference equations may also be represented as a matrix equation, where M is a tridiagonal matrix with two additional bands, as depicted in Figure 3:

Mu ) q

Figure 3. Matrix arising from finite difference discretization of a twodimensional mass transport equation.

(12)

Linear Homogeneous Kinetics. We now turn to homogeneous kinetics coupled to the electron transfer, using a chemically reversible ECE mechanism48 as an example: The steady-

A+efB

-(λard - λars)aj,k-1 - λaz aj-1,k + (2λard + 2λaz )aj,k - λaz aj+1,k - (λard + λars)aj,k+1 ) 0 -(λbrd - λbrs)bj,k-1 - λbz bj-1,k + (2λbrd + 2λbz - k1)bj,k - λbz bj+1,k - (λbrd + λbrs)bj,k+1 + k-1cj,k+1 ) 0 (15) -(λcrd - λcrs)cj,k-1 - λcz cj-1,k + (2λcrd + 2λcz - k-1)cj,k - λcz bj+1,k - (λcrd + λcrs)bj,k+1 + k1cj,k+1 ) 0 Species A is independent of species B and C and results in a matrix of the form shown in Figure 3. The rate of loss of species B is only a function of its own concentration, but it is regenerated from species C: this process links the two species. Similarly the rate at which species C is lost depends only on the concentration of C, but C is made from species B, also coupling the two species. Since species C is destroyed at the electrode (to form D), it cannot be related to B by the equilibrium constant, k1/k-1. The three species may all be mapped into u with an expanded storage map equation:

(13)

ui ) sj,k where i ) k + (j - 1)(NGX) + (s - 1)(NGX NGY) (16)

state diffusion-only material balance equations for an ECE

where s denotes a general species (for species A, s ) 1; B, s ) 2 etc.). The kinetic relationship leads to an extra pair of bands in the matrix, as shown in Figure 4. Homogeneous kinetic

BhC C+efD


Alden and Compton The Jacobian matrix may be obtained by differentiating the finite difference equation with respect to the concentration of each species. This is used as the basis of an iterative scheme: the Jacobian is solved by the linear solver to obtain a new approximation (z+1) of the concentration in u from the previous one (z) which is used to supply the concentration values for both the Jacobian, M and q.

M(u(z))u(z+1) ) M(u(z))u(z) - F(u(x)) ) q(u(z))

( )

where

Figure 4. Matrix arising from finite difference discretization of the material balance equations for an ECE mechanism.

relationships between species may be represented as a grid: A effect on A effect on B effect on C

affected by B

(21)

f1,1(u) f2,1(u) F(u) ) l fNGY,NGX(u)

and fj,k(u) denotes a general finite difference equation. For each reversible second-order reaction

n1L1 + n2L2 h m1R1 + m2R2

C

(22)

the kinetic equations for the four species are k1 -k1

-k1 k1

This relationship grid may used to formulate a generalized rule for a mechanism consisting of an arbitrary number of first-order steps. To index the grid, we introduce the notation A r B (which should be read as “B’s effect on A”). Thus B r C ) k-1. The grid is initially filled with zeros; then the effects of each reaction are summed into the grid. For each reversible firstorder reaction between two chemical species L and R

nL h mR

(17)

-

∂[L1] ) n1(k1[L1][L2] - k-1[R1][R2]) ∂t

-

∂[L2] ) n2(k1[L1][L2] - k-1[R1][R2]) ∂t

-

∂[R1] ) m1(-k1[L1][L2] + k-1[R1][R2]) ∂t

-

∂[R2] ) m2(-k1[L1][L2] + k-1[R1][R2]) ∂t

(23)

The homogeneous kinetics grid is updated by taking the partial derivative with respect to each concentration:

the grid is updated by

L r L W L r L + nk1 L1 r L1 W L1 r L1 + n1k1[L2]

L r R W L r R - nk-1

L1 r L2 W L1 r L2 + n1k1[L1] R r L W R r L - mk1 R r R W R r R + mk-1

L1 r R1 W L1 r R1 - n1k-1[R2] (18)

L1 r R2 W L1 r R2 - n1k-1[R1]

Finally, the grid elements are summed onto the matrix elements using the procedure

L2 r L1 W L2 r L1 + n2k1[L2]

Mi,j W Mi,j + S1 r S2

L2 r L2 W L2 r L2 + n2k1[L1]

(19)

L2 r R1 W L2 r R1 - n2k-1[R2]

where

L2 r R2 W L2 r R2 - n2k-1[R1]

i ) k + (j - 1)(NGX) + (S1 - 1)(NGX NGY) l ) k + (j - 1)(NGX) + (S2 - 1)(NGX NGY)

(20)

for S1 ) 1...NS; S2 ) 1...NS; j ) 1...NGY, and k ) 1...NGX, where NS is the number of species simulated. Nonlinear Homogeneous Kinetics. To extend the treatment to second-order kinetics, we must introduce a nonlinear method to treat the products of concentrations before we apply the linear solver. Since the nonlinearity is only due to quadratic terms, Newton’s method provides a simple and efficient treatment. This has been applied to electrochemical simulations in both one dimension by Rudolph37 using the FIFD algorithm and two dimensions by Alden24 using a multigrid method. For a more detailed description of the method applied to systems of finite difference equations, the reader is directed to refs 24 and 37.

R1 r L1 W R1 r L1 - m1k1[L2] R1 r L2 W R1 r L2 - m1k1[L1] R1 r R1 W R1 r R1 + m1k-1[R2] R1 r R2 W R1 r R2 + m1k-1[R1] R2 r L1 W R2 r L1 - m2k1[L2] R2 r L2 W R2 r L2 - m2k1[L1] R2 r R1 W R2 r R1 + m2k-1[R2] R2 r R2 W R2 r R2 + m2k-1[R1]

(24)



and q is updated by

qLj,k1 W qLj,k1 + n1(k1[L1][L2] - k-1[L1][L2]) qLj,k2 W qLj,k2 + n2(k1[L1][L2] - k-1[L1][L2]) qRj,k1 W qRj,k1 + m1(-k1[L1][L2] + k-1[L1][L2]) qRj,k2 W qRj,k2 + m2(-k1[L1][L2] + k-1[L1][L2])

(25)

Heterogeneous Kinetics. Heterogeneous kinetics couple species by inter-relating the concentrations at the electrode surface. This may be “worked around” using back-to-back grids28 with stencil methods such multigrid. The two main problems with back-to-back grids are that they make the matrix off-diagonally dominant at the boundaries causing solver instability14 and that each species may only be involved in one couple. As is evident from the treatment of homogeneous kinetics, the sparse matrix approach allows species to be coupled easily, and the boundary condition at the electrode surface is no exception. For species that are not electroactive, the boundary condition is a simple Neumann (no-flux) case. For each couple, the Butler-Volmer equations48 may be applied. Consider a single couple (written as a reduction):

A+ehB

(26)

The Butler-Volmer equations relate the forward and backward rate constants for electron transfer, kf and kb, to the applied potential (E):

F A + kb kb ; ka2 ) kf + DRAkb + FA kf + DRAkb + FA

kb1 )

FB + kf kf ; kb2 ) (33) kf + DRBkb + FB kf + DRBkb + FB

and the diffusion coefficient ratios are defined: DRA ) DA/DB and DRB ) DB/DA. This may be substituted into the general finite difference equation (8) at j ) 1 to eliminate the surface concentration, since the boundary value is not actually simulated. The finite difference equation at j ) 1 then becomes

(27)

∂a ) kfa0,k - kbb0,k ∂y

(28)

∂b DB ) kbb0,k - kfa0,k ∂y

(29)

Equations 29 and 29 may be converted into a pair of (firstorder) finite difference equations:

FA(a1,k - a0,k) ) kfa0,k - kbb0,k

(30)

FB(b1,k - b0,k) ) kbb0,k - kfa0,k

(31)

where FA ) DA/∆y and FB ) DB/∆y. To invoke the boundary condition, eqs 30 and 31 must be solved simultaneously for the surface concentration (at j ) 0) of each species as a function of the concentrations of species at the node above the electrode surface (at j ) 1). This results in an expression for each species:

a0,k ) ka1a1,k + ka2b1,k b0,k ) kb1b1,k + kb2a1,k

ka1 )

+ da1,ka2,k + ea1,ka1,k+1 ) qa1,k (34)

where the dimensionless potential, θ, ) nF(E - E°)/RT; E° is the standard potential of the couple, R and β are transfer coefficients (β ) 1 - R), and k0 is the standard heterogeneous rate constant. The flux equations at the electrode surface are

DA

where

aa1,ka1,k-1 + ba1,k(ka1a1,k + ka2b1,k) + ca1,ka1,k

kf ) k0 exp(-Rθ) kb ) k0 exp(βθ)

Figure 5. Matrix arising from finite difference discretization of the material balance equations for a couple (A + e h B).

(32)

Part of this transformation may be accomplished by the following sequence

ca1,k W ca1,k + ba1,kka1 ba1,k W 0

(35)

leaving the term ba1,kka2, which does not fit into the five-point stencil unless a back-to-back grid is used. This term may be easily incorporated into the matrix formulation in the same manner as for homogeneous kinetics, giving rise to a partial “band” relating species B to A, and another relating A to B, at the electrode surface as shown in Figure 5. Now let us generalize this to a procedure for a set of chemical species, some involved in couples, some electroinactive. Each couple may be represented as a reduction:

OX + e h RED

(36)

The Butler-Volmer equations for each species may be summarized into a matrix equation:

H1 ) GH0

(37)

By computing the inverse of matrix G, the surface concentrations H0 may be expressed in terms of the nodes above the surface H1. This may be achieved using a matrix inversion library subroutine such as NAG’s F01AAF based on an LU decomposition using Crout’s method.


Alden and Compton

The heterogeneous relations may be summarized as another grid. The two redox couples (A/B and C/D) in an ECE reaction serve as an example: affected by B C -kb1 kb1 kf2 -kf2

A kf1 -kf1

effect on A effect on B effect on C effect on D

point approximation as used in the previous section. This may be incorporated into our scheme with relative ease. Consider a three-point (second-order) expansion:10,50

|

-3a0,k + 4a1,k - a2,k ∂a ≈ ∂y y)0 2

D

-kb2 kb2

where kf1 and kb1 are the Butler-Volmer rate constants for the first couple (A + e h B), and kf2 and kb2 apply to the second couple. The general rule for assembling matrix G from the heterogeneous relationship grid is as follows: 1. Begin with G as the identity matrix. Species that are not involved in any couples just have unity on the diagonal; this corresponds to a Neumann boundary condition. 2. Sum the effects of each couple onto the grid:

The Butler-Volmer equations for each couple become

1 3 FA 2a1,k - a2,k - a0,k ) kfa0,k - kbb0,k 2 2

( F (2b B

[OX] r [RED] W [OX] r [RED] - kb [RED] r [OX] W [RED] r [OX] - kf [RED] r [RED] W [RED] r [RED] + kb

(38)

3. Transfer each grid element into the G matrix by

GS1,S2 )

[

S1 r S2 FS1

(39)

]

where FS1 ) DS1/∆y. For example, the two couples in an ECE reaction give rise to a G matrix of

[]

a1,k b1,k c1,k ) d1,k

FA + kf1 kb1 FA FA - kf1 FB + kb1 FB FB 0 0

0 0

0

0

0

0

FC + kf2 - kb2 FC FC - kf2 FD + kb2 FD FD

[]

1 - b2,k 2

[ ]

0,k

b 0,k

- kfa0,k

(45)

3 F + kf 2 A FA

-kb FA

-kf FB

3 F + kb 2 B FB

[ ]

1 2a1,k - a2,k 2 and H1 ) 1 2b1,k - b2,k 2

(46)

Note that a Taylor expansion with a maximum of a three points may be incorporated into the five-point stencil of most “standard” finite difference methods for two-dimensional simulations. Higher order expansions may be used to evaluate derivatives (e.g. in the flux calculation) but not to specify the boundaries. Since matrix methods are not restricted by such stencils, an n-point Taylor series may be used to specify any boundary. An n-point expansion may be represented by a series of n Taylor coefficients (T0...Tn-1); thus for a three-point expansion T0 ) -3/2, T1 ) 2, and T2 ) -1/2. The modifications to the rules for assembling G are as follows: 1. Begin with G having diagonal entries of T0. 2. The grid is formed as before. 3. Each grid element is now transferred to the G matrix using the expression

a0,k b0,k c0,k d0,k

(40)

1,k

) 3 - b ))k b 2

For a single couple G and H1 are given by

G)

[OX] r [OX] W [OX] r [OX] + kf

(44)

GS1,S2 ) -T0

S1 r S2 FS1

(47)

Once G-1 has been found, the boundary condition may be applied to each point above the electrode surface in the Taylor series using the update formula

Mi,l W Mi,l + TjbSj,k1 GS-1 1,S2

The boundary condition may be implemented generally by the procedure

(48)

where

Mi,l W Mi,l + bSj,k1 GS-1 1,S2

(41)

where

i ) k + (j - 1)(NGX) + (S1 - 1)(NGXNGY) l ) k + (j - 1)(NGX) + (S2 - 1)(NGXNGY)

i ) k + (S1 - 1)(NGX NGY) l ) k + (S2 - 1)(NGX NGY)

(42)

bS1,k ) 0

S ) 1...NS; k ) 1...NGX

for S1 ) 1...NS; S2 ) 1...NS; k ) 1...NGX; and j ) 1...n-1 This is then followed by

bS1,k ) 0 for S ) 1...NS; k ) 1...NGX

for S1 ) 1...NS; S2 ) 1...NS; and k ) 1...NGX. This is followed by

(43)

n-Point Flux Calculation. The efficiency of some simulations may be significantly improved if the flux to the electrode surface is represented as an n-point Taylor expansion of the concentration gradient2,11,50 rather than the (first-order) two-

(49)

(50)

n-Point Neumann Boundary Conditions. Neumann boundary conditions elsewhere in the simulation may also be implemented as an n-point approximation. As noted by Bieniasz50 and Gavaghan,11 the five-point central finite difference scheme is only accurate to a second-order approximation, so there is little to be gained from using more than three points to formulate the Neumann boundaries. The three-point ap-



Figure 6. Effect of LFILL on the (CGS) CPU time for a ||R||2 tolerance of 1 × 10-3.

proximation may be incorporated directly into the basic fivepoint stencil. For example, a Neumann boundary condition at NGY may be implemented by

4 1 aNGY,k ) aNGY-1,k - aNGY-2,k 3 3

∂[A] ∂2[A] ∂2[A] DA ∂[A] ) DA 2 + DA 2 + ∂t r ∂r ∂z ∂r

(52)

which gives rise to the finite difference equation t+1 t+1 (λrd - λrs)aj,k-1 + λzaj-1,k + (1 - 2λrd - 2λz)at+1 j,k t+1 t+1 + λzaj+1,k + (λrd + λrs)aj,k+1 ) atj,k (53)

where

DA∆t 2

(∆z) ,

λrd )

DA∆t 2

(∆r) ,

and λrs )

DA∆t 2∆r

(54)

and time is quantized into steps of ∆t so that the index t+1 represents the current time (to be solved for) and t represents the previous time step (for which a solution has been generated). Similarly for the mass transport to a channel microband electrode, eq 5 becomes 2

λy )

DA∆t DA∆t Vx∆t , λx ) , and λV ) 2 2 2∆x (∆y) (∆x)

(57)

(51)

However, should they be required, higher order finite difference schemes and Neumann boundary conditions may be incorporated into the sparse matrix. Extension to Time-Dependent Simulations. Under nonsteady-state conditions, eq 2 describing the mass transport to a microdisk electrode becomes

λz )

where

2

∂[A] ∂ [A] ∂ [A] ∂[A] ) DA 2 + DA 2 - Vx ∂t ∂x ∂x ∂y

(55)

This may be cast into a finite difference equation: t+1 t+1 (λx + λV)aj,k-1 + λyaj-1,k + (1 - 2λx - 2λy)at+1 j,k t+1 t+1 + λyaj+1,k + (λx - λV)aj,k+1 ) atj,k (56)

Both eqs 53 and 56 may be fitted into the five-point template equation (8) with qaj,k ) atj,k. By examining the relationship between eqs 3 and 53 or eqs 6 and 56, and generalizing, the relationship between the steadystate and transient problems may be expressed as the following rules: 1. For the steady-state case qaj,k ) 0, for the transient qaj,k ) atj,k 2. caj,k(transient) ) 1 + caj,k(steady-state) 3. λ(transient) ) λ(steady-state)∆t 4. homogeneous kinetic terms: S1 r S2(transient) ) S1 r S2(steady-state)∆t The transient simulation involves using the solution from the previous time step (t) as the right-hand side, q, in the equation for the new time step (t+1). Since the iterative solvers require a starting approximation, the solution at the previous time may also be used for this purpose. Preconditioned Krylov Subspace (F11) Methods. These methods are able to solve large sparse nonsymmetric linear systems of algebraic equations (11). Solution occurs in two stages: the first involves the generation of a preconditioner matrix P. The second stage is to solve the preconditioned system iteratively:

PNu ) q where N ) P-1M

(58)

A common approach to preconditioning is to derive P via an incomplete LU (ILU) factorization52 of the matrix M: a truncated form of Gaussian elimination. A parameter known as the “fill level” (LFILL) controls how much of the LU factorization is retained and how much is discarded. With a fill level of zero, the only decomposed elements that are retained are at points corresponding to nonzero elements in the original matrix. As the level of fill is increased, more elements of the decomposition are retained, providing a better preconditioning matrix but for a higher cost in memory and CPU time.


Alden and Compton

Figure 7. Effect of MILU on the CPU time of the three Krylov subspace solvers for a ||R||2 tolerance of 1 × 10-3.

TABLE 1: Parameters for the Model Channel Microband Simulation (See Figure 2)

i-1

30 µm 0.4 cm 0.04 cm 0.6 cm 0.01 cm3 s-1 1 × 10-6 mol cm-3 2 × 10-5 cm2 s-1 0.5 0.5 0.2 N

xe w 2h d Vf [A]bulk D Ra βa φb PSTRATc

ui ) R0 + R1 + ... + Ri-1 )

(I - M)jR0 ∑ j)0

(60)

which are in increasing Krylov subspaces

ui ∈ Ki where Ki ) {R0, R1, ..., Ri-1}

a R and β are the amounts of space (in number of electrode lengths) allowed upstream and downstream, respectively, of the electrode for the propagation of axial diffusion. b φ is the fraction of the channel simulated in the y coordinate. c PSTRAT is the pivoting strategy used by the ILU preconditioner. The value N corresponds to no pivoting (suitable for a diagonally dominant matrix).

It is possible to use a range of other methods as preconditioners; the NAG library also contains subroutines such as F11DEF which may be used for Jacobi preconditioning or SSOR preconditioning of the Krylov subspace methods. Recently Meese53 showed how SIP could be extended to act as a similar preconditioner to ILU, giving efficiency improvements. Once the preconditioner matrix, P, has been constructed, one of the three iterative Krylov subspace solvers can be applied. These work by splitting the matrix M to generate an iterative scheme for which each successive approximation of the solution is an element of increasing Krylov subspaces. The simplest splitting M ) I - (I - M) gives rise to the Richardson iteration:

ui ) (I - M)ui-1 + q

which generates solutions

(59)

(61)

where Ri is the residual at step i; Ri ) q - Mui. Note that in this section, i is used to represent the iteration number and j is used as a generic index. The general splitting

M ) P - Q ) P - (P - M)

(62)

can be rewritten as the standard splitting

N ) I - (I - N)

(63)

for the preconditioned matrix

N ) P-1M

(64)

In this way the preconditioner is incorporated into the Krylov scheme. The Krylov subspace projection methods fall into three classes:54 1. Construct ui for which the residual is orthonormal to the current subspace, namely, q - Mui ⊥ Ki(M;R0). This is known as the Ritz-Galerkin approach. 2. Identify ui for which the residual two-norm, ||q - Mui||2, is minimized over Ki(M;R0). This is, not surprisingly, known as the minimum residual approach.

TABLE 2: Boundary Conditions for the Model Channel Microband Simulation region of simulation space

definition in real space

finite difference nodes

upstream downstream upper retaining wall lower wall containing electrode before electrode at electrode after electrode

x)0 x ) (R + β + 1)xe y ) 2φh y)0 0 e x < Rxe Rxe e x e (R + 1)xe (R + 1)xe < x e (R + β + 1)xe

k)0 k ) NGX j ) NGY j)0 k ) 1 w RKE k ) RKE + 1 w (R + 1)KE k ) (R + 1)KE + 1 w NGX

boundary condition [A] ) [A]bulk ∂[A]/∂x ) 0 ∂[A]/∂y ) 0 ∂[A]/∂y ) 0 [A] ) 0 ∂[A]/∂y ) 0



Figure 8. Effect of the polynomial order, l , on the BICGSTAB(l ) CPU time for a ||R||2 tolerance of (a) 1 × 10-3 and (b) 1 × 10-6.

3. Find ui so that the residual is orthonormal to some other suitable k-dimensional subspace. This is the Petrov-Galerkin approach. The three methods (CGS, BICGSTAB(l ), and RGMRES) used in this paper are all advanced hybrids of the three basic classes. For a detailed introduction to Krylov subspace methods see ref 54. Computing. The simulation code was written in C++ calling FORTRAN subroutines for the mathematical libraries. NAG’s beta code was used for the F11 routines. Callable IDL was also used (called via ANSI C wrapper routines) to provide real-time graphics for simulations. The code was compiled (optimized for the mips2 architecture) and executed on a Silicon Graphics Indigo2 with an R4400 processor and 192 Mb of RAM. CPU times were calculated using the NAG library routine X05BAF. 3. Results and Discussion In this section we compare the efficiencies of the F11 methods to established simulation methods such as the multigrid method,

TABLE 3: Convergence Criteria for Multigrid, SIP, and the F11 NAG Library Subroutines method multigrid SIPb

convergence criterion ||R||2 max j,k

| |

Rj,k if Cj,k * 0 Cj,k

max|Rj,k| if Cj,k ) 0 j,k

F11a

||R||∞ e F(||q||∞ + ||M||∞||u||∞)

a Where F is the tolerance value specified by the variable TOL in the NAG library. b This is only for CONRES (normalized residual). CONCHN (change in normalized residual) was set to a high value so it was always met. See the NAG library documentation for more information.

SIP, and direct LU decomposition. For our model problem we simulate the transport-limited current at a 30 µm long (xe) channel microband electrode in Cartesian coordinates, with the parameters shown in Table 1. A five-point Taylor expansion was used to compute the flux, and a three-point Taylor


Alden and Compton

TABLE 4: Equivalent Convergence Thresholds for Multigrid, SIP, and the F11 NAG Library Subroutines multigrid ACC ) 1 × 10-3 NGX and NGY

multigrid

F11

33 65 97 129 161 193 225

3.31 × 10 1.21 × 10-4 3.69 × 10-4 8.68 × 10-4 1.76 × 10-4 3.78 × 10-4 6.14 × 10-4

1.01 × 10 1.73 × 10-9 2.42 × 10-9 3.79 × 10-9 5.82 × 10-10 8.86 × 10-10 1.06 × 10-9

-4

multigrid ACC ) 1 × 10-6 SIP 3.25 × 10 4.38 × 10-9 6.04 × 10-9 9.25 × 10-9 1.37 × 10-9 2.06 × 10-9 2.45 × 10-9

-8

-8

multigrid

F11

4.12 × 10 1.23 × 10-7 3.39 × 10-7 7.37 × 10-7 1.34 × 10-7 3.06 × 10-7 5.41 × 10-7

1.99 × 10 1.93 × 10-12 2.37 × 10-12 3.36 × 10-12 3.99 × 10-13 6.95 × 10-13 8.99 × 10-13

-7

SIP 4.57 × 10-11 4.90 × 10-12 5.77 × 10-12 8.00 × 10-12 9.61 × 10-13 1.66 × 10-12 2.13 × 10-12

-11

TABLE 5: Estimated CPU Times for Various Problem Sizes CPU time 104 nodesa

a

105 nodesa

106 nodesb

problem size tolerance, ||R||2

1 × 10-3

1 × 10-6

1 × 10-3

1 × 10-6

1 × 10-3

multigrid SIP BICGSTAB(l ) with MILU CGS with MILU RGMRES(10) with MILU BICGSTAB(l ): LFILL ) 5 CGS: LFILL ) 5 RGMRES(10): LFILL ) 5 BICGSTAB(l ) CGS RGMRES(10) DIRECT LU

0.5 s 2.75 s 5.2 s 5.2 s 7.3 s 6.5 s 7.3 s 8.6 s 8.5 s 18.2 s 10.0 s

0.65 s 4.3 s 6.9 s 5.8 s 8.5 s 10.6 s 9.9 s 12.4 s 12.9 s 18.4 s 18.2 s

7.5 s 56.0 s 2.2 min 2.2 min 3.4 min 4.1 min 6.1 min 6.3 min 7.1 min 9.1 min 15.0 min

9.6 s 1.4 min 3.5 min 3.2 min 5.4 min 5.4 min 6.9 min 8.2 min 8.8 min 10.8 min 20.1 min

1.8 min 18.9 min 55.7 min 54.3 min 1.6 h 2.6 h 4.4 h 5.0 h 5.9 h 4.6 h 22 h

2.6 min

3.6 h

1 × 10-6

2.3 min 26.2 min 1.7 h 1.8 h 3.4 h 2.7 h 4.8 h 5.4 h 6.0 h 6.3 h 22 h 12.3 days

b

Values were interpolated along the regression line. Value were extrapolated from the regression line.

Figure 9. Effect of the restart subspace, m, on the RGMRES CPU time for a ||R||2 tolerance of 1 × 10-6.

expansion was used to apply second-order Neumann boundary conditions. The boundary conditions for this problem are shown in Table 2. The ILU preconditioner and Krylov subspace solvers have a number of adjustable parameters which may be “tuned” for each type of problem. Fill may be introduced into the ILU preconditioner to improve the preconditioning matrix and thus reduce the number of iterations performed by the solver at the expense of a greater storage requirement. Figure 6 shows the effect of varying the level of fill in the ILU preconditioner with the CGS solver. The optimum value of LFILL is 5 for this problem. Larger values increase the time taken by the preconditioner and the solver per iteration; smaller values increase the number of iterations performed by the solver. In all cases the extra memory required for the fill was less than the memory used with no fill.

The NAG library subroutine allows the ILU preconditioner to be modified to preserve the row sums in the matrix; the option is known as modified ILU (MILU). Instead of simply discarding fill elements within a row, they are subtracted from the diagonal element. This has the effect of accelerating the solution at the cost of instability for some systems.36,42 The system of equations arising from discretization of the channel flow problem is well conditioned, so MILU may be used without any problems. The effect of using MILU is shown in Figure 7: a dramatic improvement in efficiency is observed. Also note that the optimal fill level with MILU is zero. The Krylov subspace methods BICGSTAB(l ) and RGMRES both have a single adjustable parameter which may also be “tuned”. Figure 8a,b shows the effect of varying the polynomial order, l , in the BICGSTAB(l ) solver for multigrid tolerances



Figure 10. Plot of log10(CPU time) vs log10(N) for various solvers for a ||R||2 tolerance of (a) 1 × 10-3 and (b) 1 × 10-6.

TABLE 6: Memory Requirement for a Medium (104 Node) and Very Large (106 Node) Problem memory requirement (Mb) method

104 nodes

106 nodes

multigrid SIP BICGSTAB(4) CGS RGMRES(10) DIRECT LU

1.11 0.76 2.89 2.44 2.90 1526

109 76.3 289 244 290 1.5 × 107

of 1 × 10-3 and 1 × 10-6, respectively. For the less converged solution (a) the optimum l value is 1. For (b) the optimum value seems to be about 3-4, although larger and smaller values make little difference to the CPU time. Figure 9 shows the effect of varying the size of the restart subspace, m, in the RGMRES solver for a tolerance level of 1 × 10-3. The optimum value is approximately 10 regardless of the tolerance level. The CPU time increases significantly if m is much smaller. The optimum values and behavior are very similar to that found by Salvini and Shaw for a range of problems.42

To compare CPU times of various iterative linear methods, it is necessary to specify equivalent iteration convergence criteria. These are defined differently in many of the NAG library routines, as shown in Table 3. To compare like with like, the multigrid method was used to solve the test problem for a range of meshes with a residual two-norm threshold of 1 × 10-3 and 1 × 10-6. The former is suitable for most electrochemical simulations, but the latter was used to provide a comparison for extremely converged solutions. The thresholds for the other algorithms were computed from the residuals of these solutions, as shown in Table 4. A value of 30 was used for the SIP acceleration parameter and five mesh levels were used for the multigrid solver (that is, the number of nodes was chosen so that NGX - 1 and NGY - 1 were divisible by 25). For larger problems than those specified in Table 4 it would be advisable to use a higher SIP acceleration parameter and number of grids in the multigrid scheme to ensure similar performance. Figure 10 shows the decadic log of the CPU time as a function of the decadic log of the number of nodes for (a) a multigrid tolerance of 1 × 10-3


Alden and Compton

TABLE 7. CPU Time and Memory Requirements as a Function of the Number of Nodes for Various Methods CPU time solver method direct LU SOR (optimal ω) preconditioned Krylov subspace: CGS or BICGSTAB(4) preconditioned Krylov subspace with MILU SIP multigrid a

Poisson’s equationa 2

O(N ) O(N3/2) O(N3/2) O(N5/4) not investigated O(N)

channel flow 1.9

O(N ) not investigated O(N1.6) O(N1.5) O(N1.3) O(N1.15)

memory requirement O(N2) O(N) O(N) O(N) O(N) O(N)

The data for Poisson’s equation was kindly provided by Gareth Shaw at NAG.

Figure 11. Plot of memory requirement (in Mb) vs N for various solvers.

and (b) a multigrid tolerance of 1 × 10-6. The optimal values of l for the BICGSTAB(l ) (l ) 1 for (a) and l ) 4 for (b)) and m ) 10 for RGMRES were used. Zero fill is contrasted with an optimal fill level of 5 and with the effect of MILU. The data used to generate Figure 10 at higher CPU times (104 nodes upward) were linearly regressed to generate Table 5, which shows estimated CPU times for 104, 105, and 106 node problems. Both SIP and the Krylov subspace methods lie between the multigrid method in which CPU time depends on the total number of nodes, N, as O(N1.15) and the direct LU method in which CPU time depends on N as O(N1.9). Comparing the three Krylov subspace methods, BICGSTAB is seen to be the fastest for the channel problem. In virtually all cases CGS and BICGSTAB perform significantly better than RGMRES. RGMRES has been found to be less efficient than CGS or BICGSTAB for other problems.35,36,42 Fill can reduce the CPU time by a factor of approximately 2. MILU has a more dramatic effect on the CPU time; the dependence on the number of nodes is reduced from O(N1.6) to O(N1.5). The very large (up to 106 node) problems which have been simulated by SIP12 and Multigrid24 would require approximately 3 or 30 times the CPU time, respectively, using the preconditioned Krylov subspace methods. However it should be noted that in the simulations cited simple Cartesian or one-dimensional expanding grids were used and all current evaluations, Neumann boundary conditions, and electrode flux boundary conditions were first order. Using an improved conformal mapping and higher order boundary conditions, it should be possible to reduce the number of nodes required to implement the Krylov subspace methods with CPU times comparable to the multigrid method; it is only necessary to reduce the number of nodes by a factor of four in each coordinate to achieve this. It is also possible to

use a multigrid preconditioner55,56 with the Krylov subspace methods, which would combine the advantages (speed and flexibility) of both methods. Figure 11 and Table 6 show the memory requirement (in Mb) of the various methods as a function of the number of nodes. The Krylov subspace methods require a amount of storage space that is linearly proportional to the number of nodes, a great improvement over the direct LU method, which has a quadratic dependence. The largest problem that could be simulated using direct LU factorization was a 72 × 72 node mesh using a Silicon Graphics Indigo2 with 192 Mb of RAM. Simulating a mesh of several hundred nodes in either direction presented no problem with the Krylov subspace methods even when some fill was used. The multigrid method uses approximately half, and SIP one-third, of the memory of the Krylov subspace methods. 4. Conclusions A strategy for a general electrochemical simulator has been presented. Unlike many commercial simulators currently available, this is not algorithmically constrained to one spatial dimension. Both transient experiments, including chronoamperometry, linear sweep and cyclic voltammetry, and steady-state voltammetry, may be simulated. Migration effects could also be incorporated into this sparse matrix approach by concatenating the vector of unknown concentrations with the coupled mesh of unknown potentials for each species and adding an extra band into the coefficient matrix to relate the variables. The preconditioned Krylov subspace methods available in the NAG library have been shown to be suited to solving the large sparse nonsymmetric matrixes arising from application of the finite difference method to general electrochemical



problems. Of the three solvers available in the NAG library, either BICGSTAB or CGS is recommended, depending on the desired tolerance, since RGMRES seems to be less efficient. However in terms of robustness, RGMRES > BICGSTAB > CGS,54 so the former may prove useful for difficult problems. CGS has the advantage that there are no adjustable parameters to “tune”, although BICGSTAB is more efficient when optimally tuned. The recommended values are l ) 4 for BICGSTAB(l ) and m ) 10 for RGMRES. A small amount of fill (LFILL ) 2-5) will generally accelerate solution (except with MILU). It should be noted that nonoptimal values of adjustable parameters are far less detrimental than with large SIP and especially SOR simulations.

SIP

The memory and CPU dependencies on the number of nodes (N) for the direct LU decomposition, the preconditioned Krylov subspace solvers, SOR, SIP, and multigrid methods are summarized in Table 7. If the problem allows MILU to be used (as in the case of the steady-state Cartesian simulation of the channel microband electrode), then the preconditioned Krylov subspace methods compare well with finite difference solvers such as SOR and SIP, although the high efficiency of the multigrid method should be noted. For ill-conditioned problems where MILU cannot be used, it is likely that finite difference solvers such as SIP and multigrid would diverge.

preconditioner

The stencil-based finite difference solvers are unsuitable for simulating more complex mechanisms where the material balance equations are coupled; here matrix methods are necessary. The F11 solvers compare very favorably to the direct LU factorization of the matrix proposed by Britz or even stateof-the-art direct solvers as shown by Salvini and Shaw.42

matrix method

implicit method

explicit method

hopscotch

ADI

SOR

SSOR

LU factorization

ILU

fill

Glossary stencil

condition of a matrix

The sparsity pattern that a particular method can accommodate (sometimes referred to as “finite difference molecule”). A method that treats a (linear) system as a matrix equation and can therefore accommodate any arbitrary sparsity pattern. The unknown values at a given node and the ones around it are solved from the previous value at that node (and possibly those around it) and the boundary conditions. The unknown value at a given node is solved from a combination of the previous values at that node and the ones around it and the boundary conditions. A relaxation method for solving finite difference equations that is fully explicit overall but offers greater stability than a simple explicit (Jacobi) method. Alternating direction implicit method: a relaxation method for solving finite difference equations related to the one-dimensional Crank-Nicholson method in that it contains a mixture of implicit terms in one dimension and explicit terms. Successive over-relaxation: an iterative method for solving finite difference equations, whose performance depends critically on an adjustable parameter, ω. Symmetric successive over-relaxation: a variant of SOR (essentially SOR in the natural ordering then SOR in the reverse ordering). This is commonly used as a preconditioner as it preserves the symmetry in the underlying system.

LFILL MILU

Krylov subspace method

CGS BICGSTAB(l ) RGMRES W [S1] r [S2]

Strongly implicit procedure: an iterative method for solving finite difference equations. More efficient than SOR and less critically dependent on its acceleration parameter. The sensitivity of the solution to small changes in input data. A quantitative measure is the condition number of the linear system, κ(M) ) ||M|| ||M1|| (if two norms are used, this is known as the spectral condition number). A well-conditioned system has κ(M) close to unity. An ill-conditioned system has a large κ(M), which may cause problems with numerical methods using finite precision arithmetic and lead to slow convergence of iterative methods. An auxiliary matrix in an iterative method that approximates either the coefficient matrix, M, or its inverse (M-1). The preconditioner is applied in every step of the iterative method. It is used to improve the condition of the coefficient matrix to accelerate solution by an iterative method. Solution of the linear equations system q ) Mu using Gauss-Jordan elimination to factor the matrix M into upper (U) and lower (L) triangular matrixes, which can be used to find M-1. A number of algorithms may be used for this1 including those of Doolittle, Crout, Banachiewicz, and Cholesky. Incomplete LU factorization: the LU factorization is truncated so that only elements corresponding to the original sparsity pattern of M are retained, plus extras depending on the level of fill. Commonly used as a preconditioner. A position that is zero in the matrix M, but not in the exact factorization of M. Level of fill in the ILU preconditioner. Modified ILU in order to preserve row sums in M, often gives better performance with finite difference methods. An iterative method based on splitting the coefficient matrix, M. The approximate solution vectors after each iteration form a KryloV sequence. The subspace spanned by these vectors is known as KryloV subspace. The new iterate is generated from a projection onto this subspace. Conjugate gradient squared Krylov subspace method. Biconjugate gradient stabilized Krylov subspace method (based on BICG). Restarted generalized minimal residual Krylov subspace method (based on GMRES). Assignment operator: the left operand is replaced with the right operand. Kinetic relation operator: the effect of species S2 on species S1.

Acknowledgment. We thank the Numerical Algorithms Group for allowing us to use their beta code and especially Gareth Shaw at NAG for his advice and assistance. We thank EPSRC for supercomputing time (Grant No. GR/L/36413). We also thank EPRSC for a studentship and Keble College, Oxford, for a Senior Scholarship for J.A.A. References and Notes (1) Fox, L. An Introduction to Numerical Linear Algebra; Clarendon Press: Oxford, 1964.

8954 J. Phys. Chem. B, Vol. 101, No. 44, 1997 (2) Britz, D. Digital Simulation in Electrochemistry, 2nd ed.; Springer-Verlag: Berlin, 1988. (3) Speiser, B. to Electroanalytical Chemistry; Bard, A. J., Ed.; Marcel-Dekker: New York, 1997; Vol. 19. (4) Rudolph, M. J. Electroanal. Chem. 1991, 314, 13. (5) Rudolph, M., Reddy, D. P.; Feldberg, S. W. Anal. Chem. 1994, 66, 589. (6) Bieniasz, L. K. Comput. Chem. 1997, 21, 1. (7) Flanagan, J. B.; Marcoux, L. J. Phys. Chem. 1973, 77, 1051. (8) Shoup, D.; Szabo, A. J. Electroanal. Chem. 1984, 160, 1. (9) Feldberg, S. W. J. Electroanal. Chem. 1987, 222, 101. (10) Heinze, J. J. Electroanal. Chem. 1981, 124, 73. (11) Gavaghan, D. J. J. Electroanal. Chem. 1997, 420, 147. (12) Compton, R. G.; Dryfe, R. A. W.; Wellington, R. G.; Hirst, J. J. Electroanal. Chem. 1995, 383, 13. (13) NAG Fortran Library Mark 17, Section D03EBF, p 1, Numerical Algorithms Group, Oxford. The Numerical Algorithms Group have an information desk which may be contacted by email: [email protected]. (14) Alden, J. A.; Compton, R. G. J. Electroanal. Chem. 1996, 415, 1. (15) NAG Fortran Library Mark 17, Section D03EDF, p 1, Numerical Algorithms Group, Oxford. (16) Friedrichs, M. S; Friesner, R. A.; Bard, A. J. J. Electroanal. Chem. 1989, 258, 243. (17) Kavanaugh, T. C.; Friedrichs, M. S.; Friesner, R. A.; Bard, A. J. J. Electroanal. Chem. 1990, 283, 1. (18) Bard, A. J.; Denuault, G.; Friesner, R. A.; Dornblaser, B. C.; Tuckerman, L. S. Anal. Chem. 1991, 63, 1282. (19) Friesner, R. A.; Tuckerman, L. S.; Dornblaser, B. C.; Russo, T. V. J. Sci. Comput. 1989, 4, 327. (20) Georgiadio, M.; Alkire, R. C.; J. Electrochem. Soc. 1994, 141, 679. (21) Mao, Z.; White, R. E. J. Electrochem. Soc. 1994, 141, 1. (22) Balslev, H.; Britz, D. Acta Chem. Scand. 1992, 46, 949. (23) Britz, D. J. Electroanal. Chem. 1996, 406, 15. (24) Alden, J. A.; Compton, R. G. J. Phys. Chem. B, in press. (25) Alden, J. A.; Hutchinson, F.; Compton, R. G. J. Phys. Chem. B 1997, 101, 949. (26) Micheal, A. C.; Wightman, R. M.; Amatore, C. A. J. Electroanal. Chem. 1989, 26, 33. (27) Compton, R. G.; Pilkington, M. B. G. J. Chem. Soc., Faraday. Trans. 1, 1989, 85, 2255. (28) Bidwell, M. J.; Alden, J. A.; Compton, R. G. J. Electroanal. Chem. 1996, 417, 119.

Alden and Compton (29) Amatore, C. A.; Fosset, B. J. Electroanal. Chem. 1992, 328, 21. (30) Verbrugge, M. W.; Baker, D. R. J. Phys. Chem. 1992, 96, 4572. (31) Pastore, P.; Magno, F.; Lavagnin, I.; Amatore, C. J. Electroanal. Chem. 1991, 301, 1. (32) Yu, H. G.; Smith, S. C. J. Chem. Soc., Faraday Trans. 1997, 93, 861. (33) Peskin, U.; Miller, W. H. Edlund, A. J. Chem. Phys. 1995, 103, 10030. (34) Tao, G. H.; Wyatt, R. E. Chem. Phys. Lett. 1995, 239, 207. (35) Tsai, C. C.; Liu, T. J. J. Non-Newtonian Fluid Mech. 1995, 60, 155. (36) Yang, D. Y.; Chen, G. S; Chou, H. P. Ann. Nucl. Energy 1993, 20, 9. (37) Rudolph, M. J. Electroanal. Chem. 1992, 338, 85. (38) Sonneveld, P. SIAM J. Sci. Statist. Comput. 1989, 10, 36. (39) Van der Vorst, H. SIAM J. Sci. Statist. Comput. 1989, 13, 631. (40) Sleijpen, G. L. G.; Fokkema, D. R. Electron Trans. Numer. Anal. 1993, 1, 11. (41) Saad, Y.; Schultz, M. SIAM J. Sci. Statist. Comput. 1986, 7, 856. (42) Salvini, S.; Shaw, G. NAG Technical Report TR2/96. (43) Meijerink, J.; Van der Vorst, H. Math. Comput. 1977, 31, 148. (44) Aoki, K. Electroanalysis 1993, 5, 627. (45) Saito, Y. ReV. Polarogr. 1968, 15, 177. (46) Compton, R. G.; Unwin, P. R. J. Electroanal. Chem. 1986, 205, 1. (47) Compton, R. G.; Fisher, A. C.; Wellington, R. G.; Dobson, P. J.; Leigh, P. A. J. Phys. Chem. 1993, 97, 10410. (48) Testa, A. C.; Reinmuth, W. Anal. Chem. 1961, 33, 1320. (49) Bard, A. J.; Faulkner, L. R. Electrochemical Methods; Wiley: New York, 1980. (50) Bieniasz, L. K.; Osterby, O.; Britz, D. Comput. Chem. 1995, 19, 351. (51) Randles, J. E. B. J. Chem. Soc., Faraday. Trans. 1948, 44, 322. (52) Meijerink, J.; Van der Vorst, H. J. Comput. Phys. 1991, 44, 134. (53) Meese, E Copper Mountain Conference of IteratiVe Methods, Proceedings Vol 1, 1996, University of Colorado. Also see http:// www.termo.unit.no/mtf/people/eammtf/cm abstr eng.html. (54) Van der Vorst, H. Lecture Notes on IteratiVe Methods, http:// www.math.ruu.nl/people/vorst/cgnotes/cgnotes.html. (55) Axelsson, O.; Vassilevski, P. S. Numer. Math. 1989, 56, 157. (56) Axelsson, O.; Vassilevski, P. S. SIAM J. Numer. Anal. 1990, 57, 1569.

A General Method for Electrochemical Simulations. 1. Formulation of

Recommend Documents