Cyclic Voltammetric Current Functions Determined with a Prescribed

Anal. Chem. 2008, 80, 9659–9665

Cyclic Voltammetric Current Functions Determined with a Prescribed Accuracy by the Adaptive Huber Method for Abel Integral Equations Lesław K. Bieniasz* Institute of Physical Chemistry of the Polish Academy of Sciences, Department of Complex Systems and Chemical Processing of Information, ul. Niezapominajek 8, 30-239 Cracow, Poland, and Institute of Teleinformatics, Faculty of Electrical and Computer Engineering, Cracow University of Technology, ul. Warszawska 24, 31-155 Cracow, Poland Modern electroanalytical applications of cyclic voltammetry require that theoretical current functions be obtainable automatically, efficiently, and with a prescribed accuracy, by computer simulation algorithms. One of the classical simulation approaches relies on formulating and solving relevant integral equations. Numerical solution methods, used for this purpose so far, are nonautomatic and do not provide information about accuracy of the results. The adaptive variant of the Huber method, developed by the present author, can generate theoretical cyclic voltammograms automatically, with a given target accuracy, and more efficiently than the formerly studied patch-adaptive direct simulation method. This is demonstrated using examples of cyclic voltammograms described by the firstkind Abel integral equations. The method is therefore a candidate for automatic integral equation solvers, needed for building a new generation of problem solving environments for electroanalytical chemistry and for widely understood automation of electroanalytical investigations. Out of a variety of contemporary transient electrochemical techniques, cyclic voltammetry (CV) is one of the most frequently used, and most robust, experimental methods suitable for electroanalytical purposes.1 However, for the successful use of the method one needs to know precisely the theoretical current functions associated with a particular electrochemical reaction system, and these are typically not available in a closed form. Therefore, some kind of computer modeling,2 or digital simulation,3,4 is necessary. One of the oldest approaches to such simulations, relevant for electrochemical systems described by kinetic partial differential equations (PDEs) of the reaction-diffusion type, relies on converting the PDEs into integral equations (IEs), by means of the Laplace transform, and numerically solving the IEs. For a * Corresponding author. Phone: (+48 12) 639 52 12. Fax: (+48 12) 425 19 23. E-mail: [email protected]. (1) See, e.g.: Speiser, B. In Encyclopedia of Electrochemistry; Bard, A. J., Stratmann, M., Unwin, P. R., Eds.; Wiley: Weinheim, Germany, 2003; Vol. 3, Chapter 2.1, p 81. (2) Mocak, J. In Encyclopedia of Analytical Science; Worsfold, P. J., Townshend, A., Poole, C. F., Eds.; Elsevier: Oxford, U.K., 2005; p 208. (3) Feldberg, S. W. In Historical Perspectives on the Evolution of Electrochemical Tools; Leddy, J., Birss, V., Vany´sek, P., Eds.; The Electrochemical Society: Pennington, NJ, 2004; p 189. (4) Britz, D. Digital Simulation in Electrochemistry, 3rd ed.; Springer: Berlin, Germany, 2005. 10.1021/ac801412f CCC: $40.75  2008 American Chemical Society Published on Web 11/12/2008

guide, readers are referred to the known literature,5-8 focused on the applications of the method to spatially one-dimensional PDEs. The IE approach to simulations played an important role in the early days of the development of the CV theory, enabling, for the first time, the solution of a number of elementary kinetic models, probably best exemplified by the seminal Nicholson and Shain paper.5 The main weakness of the IE approach is that it is restricted to linear PDEs, which precludes its application to models involving nonlinear homogeneous reaction terms arising, for example, in the case of second- and higher-order reactions or for Michaelis-Menten kinetics. This is one of the reasons why digital simulations by means of direct PDE solution methods, such as finite-difference, finite-element, etc., methods3,4 are currently more widely used and developed. However, in cases when it is feasible, the IE approach can still be useful, either as a method of obtaining additional insights into the theory or as an alternative way of computing theoretical CV currents. Simulationists may also want to compare numerical solutions obtained by various methods, to verify their correctness. The derivations of the IEs sometimes facilitate an identification of the small number of lumped parameters on which the current functions effectively depend. The lumped parameters are analytical expressions, typically involving several independent model parameters occurring in the PDE system. The knowledge of the lumped parameters obviates the need to study the effects of each of the independent model parameters separately and improves our understanding of the model dependences. For example, in the case of a Nernstian charge transfer, the IE reveals, in a natural way,5 that the cyclic voltammetric current is a function of E - E1/2 rather that E - E0, where E is electrode potential, E0 is a formal potential of the charge-transfer reaction, and E1/2 is a related half-wave potential. This suggests the use of Estart - E1/2 as a lumped model parameter, instead of Estart - E0, where Estart is the potential at which the CV scan begins. This result cannot be easily obtained just from the inspection of the PDEs and initial/boundary conditions, or from the normalization of variables involved, in the direct simulation approach. The IE approach has also a potential for being computationally less expensive than the direct simulation, because (5) Nicholson, R. S.; Shain, I. Anal. Chem. 1964, 36, 706. (6) Nicholson, R. S.; Olmstead, M. L. In Computers in Chemistry and Instrumentation; Mattson, J. S., Mark, H. B., Jr., MacDonald H. C., Jr., Eds.; Marcel Dekker: New York, 1972; Vol. 2, p 119. (7) Mirkin, M. V.; Nilov, A. P. Comput. Chem. 1991, 15, 55. (8) Mircˇeski, V. J. Electroanal. Chem. 2003, 545, 29.

Analytical Chemistry, Vol. 80, No. 24, December 15, 2008

9659

in the IE approach one computes the current function(s) only, without determining spatio-temporal concentration profiles. The relative performance of the two kinds of simulations will depend, however, on the particular numerical methods used. In connection with the above facts it has to be noted that the current level of advancement of the simulation methods used to solve IEs in electrochemistry is rather low. Elementary, fixedgrid product-integration methods,9 such as the Huber method,10 have been used so far.5-8 This refers also to some of the existing software for electrochemical simulations, notably the ELSIM program.11 These methods are nonautomatic, i.e., they have no built-in control of the accuracy of the solutions obtained. Hence, the burden of assessing the quality of the solutions obtained remains with the method users, and it requires well-designed and tedious trial-and-error tests and result comparisons. This situation contrasts with the current development of adaptive methods for the direct solution of electrochemical kinetic PDEs,12-15 where the ability to automatically compute model responses with a prescribed accuracy is the essential feature of the methods. Without having analogous automatic methods for the solution of the IEs, it is difficult to think of entering a new phase of computational electrochemistry,16 and building an advanced generation of problem solving environments (software) for electroanalytical modeling, from which the users would expect a high degree of automation.17,18 Trends toward automation are also pervading the modern electroanalytical equipment and methodology.19-21 Simulation techniques have to follow these trends. Owing to its aforementioned limitations, the IE approach cannot serve as a unique or basic simulation method in a general problem solving environment for electrochemistry. However, this type of software is expected17,18 to supply the users not only with an automatic operation of the numerical procedures but also with a comprehensive collection of diverse methods, from which a method most optimal for solving a particular problem can be selected automatically or by the user. Therefore, the availability of an adaptive IE method could contribute to the flexibility and power of such software. Having in mind the above requirements, in the previous studies of the present author22,23 an adaptive product-integration method for solving the first-kind Abel IEs has been elaborated. The purpose of the present communication is to demonstrate that the method can be used for obtaining theoretical CV current functions with a given target accuracy. This will be shown using examples of electrochemical IEs of the Abel type. Further work is planned to extend the method onto other types of IEs pertinent to the theory of CV, including all the types previously solvable by the ELSIM program.11

provide final equations only. The method can be currently applied to electrochemical IEs that take the form of the first-kind Abel equations:

∫ K(t, τ)χ(τ) dτ ) F(t) t

0

in which χ(t) is the unknown CV current function dependent on dimensionless time t (often denoted by at in the literature5), F(t) is a known function, and K(t, τ) ) H(t, τ)(t - τ)-1/2 is the known kernel of the integral transformation. Two most often encountered functions H(t, τ) are H(t, τ) ) 1 and H(t, τ) ) exp[-λ(t - τ)] with λ > 0. We assume that F(0) ) 0, since if the opposite is true, the solution of eq 1 can be decomposed into a sum χ(t) ) φ(t) + ψ(t), where φ(t) satisfies

(9) Baker, C. T. H. The Numerical Treatment of Integral Equations; Clarendon Press: Oxford, U.K., 1977. (10) Huber, A. Monatsschr. Math. Phys. 1939, 47, 240. (11) Bieniasz, L. K. Comput. Chem. 1992, 16, 11. (12) Bieniasz, L. K. Electrochim. Acta 2007, 52, 3929. (13) Gavaghan, D. J.; Gillow, K.; Süli, E. Langmuir 2006, 22, 10666. (14) Henley, I. E.; Fisher, A. C.; Compton, R. G.; Banks, C. E. J. Phys. Chem. B 2005, 109, 7843.

9660


∫ K(t, τ)φ(τ) dτ ) F(0)

(2)

∫ K(t, τ)ψ(τ) dτ ) F(t) - F(0)

(3)

t

0

whereas ψ(t) obeys t

0

where the right-hand side is zero at t ) 0. Consequently, φ(t) may contain a singularity at t ) 0; for example, in the case of H(t, τ) ) 1 eq 2 yields the Cottrell-type expression for the electric current associated with the potential step to the starting potential: φ(t) ) F(0)π-1t-1/2

(4)

whereas ψ(t) is not singular. In electrochemical Abel equations F(0) is often very close to zero, but not exactly zero. Component φ(t) is then neglected, and for the purpose of the numerical calculations of the component ψ(t), eq 1 is preferably transformed to the form 3. Such a transformation has been applied in this work in examples 1 and 3 described below. The present adaptive method22,23 combines the Huber-type discretization10 with an automatic integration step selection based on controlling local error estimates. Thus, the current function is replaced by a piecewise linear spline function, spanned on the generally nonuniform grid of nodes tn (n ) 0, 1, ..., N) with step sizes hn ) tn - tn-1 (n ) 1, ..., N). The approximate nodal values χn of the current function are determined, one after another, from the formulas

χn ) METHOD The derivations of the adaptive method are described in detail in separate articles,22,23 so they will not be repeated here. We

(1)

{

2b1/2F(t1/2) - b1F(t1) for n)0 2a1/2 - a1 - 1 anχn-1 - bngn for n ) 1, ..., N

(5)

where t1/2 ) (t0 + t1)/2 with t0 ) 0,

gn )

{

n-1

F(tn) -

∑ k)1

[(

F(tn) for n ) 1 χk-1 - tk-1

)

χk - χk-1 Q0,n,k-1,k + hk

]

χk - χk-1 Q1,n,k-1,k for n ) 2, ..., N (6) hk

{ {

{

t1/2Q0,1/2,0,1/2 - Q1,1/2,0,1/2 for n ) 1/2 t0Q0,1/2,0,1/2 - Q1,1/2,0,1/2 an ) tnQ0,n,n-1,n - Q1,n,n-1,n for n ) 1, ..., N tn-1Q0,n,n-1,n - Q1,n,n-1,n

est(δn) ) (7)

h1/2 for n ) 1/2 t0Q0,1/2,0,1/2 - Q1,1/2,0,1/2 bn ) hn for n ) 1, ..., N tn-1Q0,n,n-1,n - Q1,n,n-1,n

(8)

where γn ) hn-1/hn is the local grid ratio; χˆ 2 is a provisory solution obtained assuming t2 ) 2h1. The provisory solution serves only for approximating the second derivative at t1 and is later discarded. Coefficients cn and d1/2 are additional method coefficients22,23 defined by

(9)

1 2 cn ) bnhn-2[(tn-1 + hntn-1)Q0,n,n-1,n - (2tn-1 + hn)Q1,n,n-1,n + 2 for n ) 1, ... , N (11) Q2,n,n-1,n]

with

Qm,n,l,k )

∫

tk

tl

K(tn, τ)τm dτ

Hence, Qm,n,l,k, an, and bn are analytically determined method coefficients dependent on the particular function H(t, τ) and step sizes hn. Tables in refs 22 and 23 provide expressions for these coefficients, relevant for functions H(t, τ) ) 1 and H(t, τ) ) exp[-λ(t - τ)]. Readers may wish to consult ref 22 in the first place, as it provides the derivations most important for electrochemical IEs. The second paper23 describes an addition to the method, related to the automatic determination of the starting solution values and their errors. The starting values are normally equal to zero in electrochemistry. The approximate solutions χn differ from the true solutions χ(tn) by errors δn ) χ(tn) - χn. Estimates est(δn) of these errors are obtained at the following assumptions. First, the errors are regarded as local, with the exception of the first two errors δ0 and δ1. This means that while calculating the error of χn for n > 1 all previous nodal values χk for k ) 0, 1, ..., n - 1, and the current function integrals over intervals [0, tn-1], are considered as exact. Owing to the typically fairly strong error damping on passing from one discrete time level to a next one, such local errors approximate the true global errors quite well, as will be shown below. The first two errors δ0 and δ1 are determined simultaneously so that they are global errors, i.e., they involve complete contributions from error sources at both time levels. Second, the local error estimates are assumed equal to the first, dominant terms of the Taylor expansions of the exact local errors. Third, second derivatives of the current function, occurring in the dominant error expansion terms, are replaced by three-point finite-difference quotients, using computed discrete solutions. Thus, the equations for the error estimates are (15) Ludwig, K.; Speiser, B. J. Electroanal. Chem. 2006, 588, 74. (16) Bieniasz, L. K. In Modern Aspects of Electrochemistry; Conway, B. E., White, R. E., Eds.; Kluwer/Plenum: New York, 2002; Vol. 35, p 135 and references therein. (17) Gallopoulos, E.; Houstis, E.; Rice, J. R. Future Research Directions in Problem Solving Environments for Computational Science; Technical Report CSD TR-92-032; Computer Science Department, Purdue University: West Lafayette, IN, 1992. (18) Houstis, E. N.; Rice, J. R. Math. Comput. Simul. 2000, 54, 243. (19) Bond, A. M.; Sˇvestka, M. Collect. Czech. Chem. Commun. 1993, 58, 2769. (20) Richards, E.; Bessant, C.; Saini, S. Sens. Actuators, B 2003, 88, 149. (21) Erichsen, T.; Reiter, S.; Ryabova, V.; Bonsen, E. M.; Schuhmann, W.; Ma¨rkle, W.; Tittel, C.; Jung, G.; Speiser, B. Rev. Sci. Instrum. 2005, 76, 062204. (22) Bieniasz, L. K. Computing 2008, 83, 25. (23) Bieniasz L. K. Computing, DOI 10.1007/s00607-008-0020-9.

c1 - 2d1/2 for n)0 2a1/2 - a1 - 1 for n)1 a1 est(δ0) + (χ0 - 2χ1 + χˆ 2)c1 2 2 2 χ - χ + χ c γn(1 + γn) n-2 γn n-1 (1 + γn) n n for n ) 2, ..., N (10) (χ0 - 2χ1 + χˆ 2)

[

]

and d1/2 )

h1Q1,1/2,0,1/2 - Q2,1/2,0,1/2 4h1Q1,1/2,0,1/2

(12)

Equation 10 was derived22,23 assuming that the solution χ(t) was sufficiently regular at t > 0, for all Taylor expansions used in the error derivations to exist. However, in the CV simulations considered in the present work we have to account for the possible nonsmoothness of the χ(t) function at the potential sweep switching moments (hereafter denoted by ts). As the second derivatives d2χ(ts)/dt2 do not exist, eq 10 does not give satisfactory predictions in the case of time-stepping over ts, and it has to be modified. It can be shown, in the way similar to that in Bieniasz,23 that if tn-1 is forced to coincide with ts, then the valid error estimate at the next time level is est(δn) ) (χn-1 - 2χn + χˆ n+1)cn

(13)

where χˆ n+1 is a provisory solution obtained assuming hn+1 ) hn. In other words, in such cases the second derivative is evaluated at tn rather than at tn-1, where it does not exist. Given the error estimates, the control theoretic algorithm of Gustafsson24 is used to adjust the step sizes so that the requested accuracy criteria are met. The criteria are S max{|est(δ0)|, |est(δ1)|} e tol

for n ) 1

(14)

|est(δn)| e tol

for n > 1 and tn-1 * ts

(15)

S|est(δn)| e tol

for n > 1 and tn-1 ) ts

(16)

where tol is a prescribed absolute error tolerance, and S is a heuristic safety factor, chosen to account for the discrepancy between the true and estimated errors, which may occur at t ) 0 and t ) ts, due to the current function irregularities.23 The value S ) 10 usually proves sufficient and has been used in this work. On the first time step, h1 ) hstart is initially attempted, where hstart is a predefined step value, expected to be optimal (24) Gustafsson, K. ACM Trans. Math. Software 1994, 20, 496.


9661

for a given IE, and the criterion 14 is checked, possibly resulting in h1 modifications. A maximum, IE-dependent limit hmax for the steps is also introduced, to avoid excessively large steps to be taken. EXAMPLES We demonstrate the operation of the adaptive method using three examples of eq 1, chosen because they possess semianalytical solutions that can be used for result verification. Example 1: CV current function for the reversible (Nernstian) charge-transfer reaction A + e- h B at a planar macroelectrode. The original IE,5 corrected for the initial singularity according to eq 3, is dτ ∫ (tχ(τ) - τ) t

0

1/2

)

1 1 1 + exp[u - S(t, ts)] 1 + exp(u)

(17)

where

S(t, ts) )

{

t for t e ts 2ts - t for t > ts

(18)

and u is a dimensionless parameter related to the difference between the starting potential and the half-wave potential of the charge-transfer reaction (cf. the discussion in the introduction). The u parameter must be taken sufficiently large for the effect of u on χ(t) to be much smaller than the assumed tol parameter. A number of authors presented semianalytical solutions of eq 17, in the form of quadratures,25-28 truncated series,29-35 or other approximate expressions.36-38 Out of these, the series summation proposed by Myland and Oldham32 is probably the most convenient for computing reference solutions, and it has been adopted here, because the summation algorithm has been published, and the series solution can be easily obtained for any t value. However, the truncated Myland and Oldham series has a limited accuracy. By comparing the Myland and Oldham series summation for the forward CV peak, with the apparently most accurate results available in the literature, taken from Table 1 in Mocak and Bond,35 we have found that a sum of 1000 terms of the Myland and Oldham series differs from the Mocak and Bond results by no more than 3.5 × 10-8, when u ) 20 and t ∈ [0, 30]. Unfortunately, Mocak and Bond35 used an undisclosed summation acceleration algorithm embedded in the MATHEMATICA program, which precludes its duplication in a stand-alone program. Summing up more terms of the Myland and Oldham series does not bring a visible improvement in accuracy, when the calculations are done in extended (25) (26) (27) (28) (29) (30) (31) (32) (33) (34) (35) (36) (37) (38)

Sˇevcˇik, A. Collect. Czech. Chem. Commun. 1948, 13, 349. Ramamurthy, A. C.; Rangarajan, S. K. Electrochim. Acta 1981, 26, 111. Lether, F. G.; Shivpuri, V. Comput. Chem. 1990, 14, 177. Natarajan, A.; Mohankumar, N. Comput. Chem. 1997, 21, 315. Reinmuth, W. H. Anal. Chem. 1961, 33, 1793. De Vries, V. T.; Van Dalen, E. J. Electroanal. Chem. 1963, 6, 490. Oldham, K. B. J. Electroanal. Chem. 1979, 105, 373. Myland, J. C.; Oldham, K. B. J. Electroanal. Chem. 1983, 153, 43. MacLeod, A. J. Appl. Math. Comput. 1993, 57, 305. Mocak, J. Electrochem. Commun. 2002, 4, 803. Mocak, J.; Bond, A. M. J. Electroanal. Chem. 2004, 561, 191. Lether, F. G.; Evans, O. M. Appl. Math. Comput. 1985, 16, 253. Sivakumar, S.; Basha, C. A. Russ. J. Electrochem. 2005, 41, 421. Montella, C. J. Electroanal. Chem. 2008, 614, 121.

9662


Figure 1. Adaptively simulated current functions (O) compared with semianalytical solutions (solid curves) in examples 1 (a), 2 (b), and 3 (c). Simulation parameters are u ) 20, ts ) 30, hstart ) hmax ) 1, tol ) 10-4 (in example 1); ts ) 2, hstart ) hmax ) 0.1, tol ) 10-4 (in example 2); λ ) 1012, u ) 20, ts ) 30, hstart ) hmax ) 1, tol ) 2 × 10-9 (in example 3).

Figure 3. Variations of the adaptively adjusted local step sizes hn in the course of the integration, obtained for examples 1 (a), 2 (b), and 3 (c). Simulation parameters are as in Figure 1. Figure 2. Estimates est(δn) of the absolute local errors of the current functions (-O-), compared with the true absolute errors (-b-), obtained for examples 1 (a), 2 (b), and 3 (c). Simulation parameters are as in Figure 1.

precision (“C++” (long double variables). Therefore, the sums of 1000 terms have been computed here, and the verifications of the accuracy of the adaptive solutions have been restricted to tol g 10-6, although the adaptive algorithm works without failures caused by machine error interference22 for tol g 10-7.

Example 2: CV current function for the reversible deposition A + e- h Bdep of metals on a solid planar macroelectrode, assuming a constant activity of Bdep. The current function can be shown to obey the IE:

∫ (t - τ) t

0

-1/2

χ(τ) dτ ) 1 - exp[-S(t, ts)]

(19)

(39) Berzins, T.; Delahay, P. J. Am. Chem. Soc. 1953, 75, 555.


9663

Figure 4. Largest absolute values of the true absolute errors of the adaptively simulated current functions, as functions of the error tolerance parameter tol, in examples 1 (-b-) and 2 (-O-). Other simulation parameters are u ) 20, ts ) 30, hstart ) hmax ) 1 (in example 1); ts ) 2, hstart ) hmax ) 0.1 (in example 2).

Figure 6. Largest absolute values of the true absolute errors of the adaptively simulated current functions, as functions of the computational time (ct), in examples 1 (-b-) and 2 (-O-). Other simulation parameters are as in Figure 4. The computational times are averages over 500 identical runs of the program.

diffusion coefficients of A and B, the original IE,5,41,42 corrected for the initial singularity, is 1 - τ)] 1 χ(τ) dτ ) ∫ exp[-λ(t 1 + exp[u - S(t, t )] 1 + exp(u) (t - τ) t

1/2

0

s

(21) where λ has a meaning of a dimensionless rate constant of the homogeneous reaction. In the limit of λ f ∞ (extremely fast homogeneous reaction) a steady state is approached, for which an analytical solution exists.43 The steady-state solution for eq 21 becomes χ(t) )

1 ( πλ ) { 1 + exp[u1- S(t, t )] - 1 + exp(u) } 1/2

(22)

s

Figure 5. Largest absolute values of the true absolute errors of the adaptively simulated current functions, as functions of the total number N of integration steps needed, in examples 1 (-b-) and 2 (-O-). Other simulation parameters are as in Figure 4.

from which the semianalytical solution can be derived by using tables of Laplace transforms:40 χ(t) )

{

-1

2π

2π-1 daw(t1/2) for t e ts ) - π-1/2 exp(-ts){2π-1/2 daw[(t - ts)1/2]-

1/2

daw(t

erex[(t - ts)1/2] + exp(t - ts)} for t > ts (20) where daw(x) is the Dawson integral and erex(x) ) exp(x ) erfc(x). Example 3: CV current function for the catalytic reaction mechanism A + e- h B, B f A involving a reversible charge transfer at a planar macroelectrode and a (pseudo)first-order irreversible homogeneous reaction. At the assumption of equal 40

2

(40) Abramowitz, M.; Stegun, I. A. Handbook of Mathematical Functions; Dover Publications: New York, 1972.

9664


COMPUTATIONAL DETAILS Computational experiments were performed on an IBMcompatible personal computer having an Intel Pentium D processor operating at 3 GHz, and having a 2 GB operational memory. The numerical code was written in C++ using mostly extended precision (long double variables) and compiled as a 32-bit console application, using Borland C++ Builder 6.0 compiler. The code was run under MS Windows XP Professional. Function erex(x) was calculated with a relative error of about 10-16, by using a C translation of the packet calerf of Cody.44 Dawson integral daw(x) was calculated with a similar relative error by using the algorithm of Cody et al.45 RESULTS AND DISCUSSION Figures 1-3 present typical results obtained. As can be seen in Figure 1, a good qualitative agreement between the adaptive and analytical solutions can be easily achieved. The grid nodes (41) Save´ant, J. M.; Vianello, E. Advances in Polarography; Pergamon Press: London, 1960; p 367. (42) Bieniasz, L. K. J. Electroanal. Chem. 1984, 170, 77. (43) Miller, S. L. J. Am. Chem. Soc. 1952, 74, 4130. (44) Cody, W. J. Math. Comput. 1969, 23, 631. (45) Cody, W. J.; Paciorek, K. A.; Thacher, H. C., Jr. Math. Comput. 1970, 24, 171.

clearly concentrate in critical regions requiring an increased resolution, i.e., in regions of a high curvature of the current function plots, and in regions adjacent to the points of nonsmoothness of the current functions. Figure 3 displays more quantitatively the local variations of the adapted grid spacing, showing that larger variations (over about 6 orders of magnitude) occur in the case of nonsmoothness in example 2 than in examples 1 and 3 (grid spacing variations over 1 order of magnitude). The need for adaptive adjustments of the grid is thus strongly documented. Figure 2 demonstrates that the local error estimates 10 and 13 approximate well the true global errors in examples 1 and 2. Only sporadically larger discrepancies between est(δn) and δn are observed in regions of rapid temporal changes of the errors (see, for example, the spike of the true error at t ≈ 43.37 in Figure 2a). The surprisingly small errors achieved in example 3 (cf., Figure 2c) result from special properties of the function H(t, τ) ) exp[-λ(t - τ)]. According to eq 22, the current function and its derivatives increase proportionally to λ1/2, with increasing λ, so that in the limit of large λ (note that Figures 1c, 2c, and 3c correspond to about the largest physically meaningful λ) one intuitively expects a difficulty with obtaining absolute errors comparable in magnitude with the errors in examples 1 and 2, where the current function values are of the order of unity or smaller. However, this is a wrong expectation. The coefficient cn for H(t, τ) ) exp[-λ(t - τ)] equals22

cn )

(3/4 - λhn/2)π1/2 erf[(λhn)1/2] - (3/2)(λhn)1/2 exp(-λhn) λhn{(1 - 2λhn)π1/2 erf[(λhn)1/2] - 2(λhn)1/2 exp(-λhn)} (23)

so that for large λ it behaves like cn ∼ λ-1, implying that the error estimates 10 and 13 tend to zero with increasing λ, like est(δn) ∼ λ-1/2. Hence, paradoxically, much smaller absolute errors are achievable for large λ than for small λ. Another interesting result seen in Figure 2c is the difference between the true and estimated absolute errors, occurring in the regions of rapid changes of the current function. The absolute difference there reaches a maximum of about 7 × 10-8, which corresponds to the relative current function error of about 3 × 10-13. The difference does not vanish when tol is decreased or when a sequence of fixed, uniform grids, with a decreasing step size is applied instead of the adaptive grid. This result suggests that the steady-state solution (eq 22) differs slightly from the true transient current function, even for λ as large as 1012, so that the adaptive solution is, in fact, more accurate than eq 22. To the knowledge of the present author, such subtle theoretical effects have not been revealed in any CV simulation studies of the catalytic mechanism published so far. The adaptive Huber method thus appears quite powerful in revealing the peculiarities of electrochemical kinetics. Of course, the effect is too small to be detected experimentally, as the typical accuracy of experimental results is much worse (relative error about 1%).

Parts a and b of Figure 2 confirm that a particular target accuracy tol is achieved in examples 1 and 2 (results for example 3 are less informative, since eq 22 is not a sufficiently accurate reference solution, as was discussed above). Figure 4 further documents this feature, for tol varying over several orders of magnitude. As can be seen, the largest absolute error of the simulated current function increases nearly linearly with the tol assumed, and its values are close to tol. Such a proportionality of the error, to the tolerance parameter, is generally required from adaptive solvers.46 Figure 5 reveals that a change of the number N of integration steps by 1 order of magnitude roughly corresponds to the change of the errors by 2 orders of magnitude. This result is consistent with the theoretical second-order accuracy, predicted for the method,22 since an average step size is inversely proportional to N. Finally, Figure 6 presents efficiency plots for the adaptive Huber method. From the slope of the nearly linear logarithmic plots we find that the computational time is roughly inversely proportional to the error achieved. This result is consistent with the second accuracy order and with the fact that the number of arithmetic operations in the adaptive algorithm is proportional to N2, since eq 5 is essentially a linear algebraic equation system with a triangular matrix. By comparing Figure 6 with Figure 1 in ref 12 we note that computational times needed for adaptive solutions of eq 1 are much smaller (roughly at least by 1 order of magnitude) than the computational times characteristic of the patch-adaptive method, when directly solving kinetic PDEs with a comparable accuracy. Hence, the adaptive Huber method provides the solutions more efficiently. We stress that this result refers to the comparison of adaptive methods. It is possible that by using suitably designed fixed grids (for example, exponentially expanding spatial grids), a conventional direct simulation will be faster than the present Huber method. However, such fixed-grid calculations do not provide error estimates and are nonautomatic: the grids must be predefined by a simulationist, based on the a priori knowledge of the solutions. In contrast, the adaptive space grids are generated by the patch-adaptive algorithm, which requires some extra cost, but the a priori knowledge of the concentration profiles is not required. CONCLUSIONS We conclude that the adaptive Huber method is a useful simulation tool for electroanalytical chemistry, which can be used for obtaining CV current functions efficiently and with a prescribed accuracy. The method is therefore a good candidate for automatic IE solvers, needed for building a new generation of problem solving environments for electroanalytical chemistry and for the widely understood automation of electroanalytical investigations. Received for review July 9, 2008. Accepted September 26, 2008. AC801412F (46) Shampine, L. F. Lect. Notes Math. 1989, 1386, 118.


9665

Cyclic Voltammetric Current Functions Determined with a Prescribed

Recommend Documents