A reliable automatic simulation of singular electroanalytical transients, by the adaptive Huber method for Volterra integral equations

Journal of Electroanalytical Chemistry 799 (2017) 40–52

Contents lists available at ScienceDirect

Journal of Electroanalytical Chemistry journal homepage: www.elsevier.com/locate/jelechem

A reliable automatic simulation of singular electroanalytical transients, by the adaptive Huber method for Volterra integral equations

MARK

L.K. Bieniasz Faculty of Physics, Mathematics, and Computer Science, Cracow University of Technology, ul. Warszawska 24, Cracow 31-155, Poland

A R T I C L E I N F O

A B S T R A C T

Keywords: Volterra integral equations Weakly singular kernels Weakly singular solutions Adaptive methods Product-integration Computational electrochemistry

In a number of transient electroanalytical experiments (such as, in particular, the potential step chronoamperometry) one expects current-time responses having a singularity at time t = 0. A reliable simulation of such singular responses is a difficult task. Conventional simulation techniques fail to provide accurate results close to the singularity. A recent extension [L. K. Bieniasz, J. Comput. Appl. Math. 323 (2017) 136] of the adaptive Huber method for solving Volterra integral equations (IEs) is shown to overcome this difficulty in cases when the current behaves like t −1/2 for t → 0. The method requires an availability of highly accurate approximants for computing certain integrals of the kernel functions occurring in the IEs. Relevant approximants are elaborated for several most important kernel functions specific of one-dimensional diffusion, and one kernel function for two-dimensional diffusion. The resulting algorithm is tested on two examples of single IEs describing chronoamperometry and cyclic voltammetry for a single reversible charge transfer, and one example of an IE system describing chronoamperometry for an ErevErev mechanism. Singular transients are simulated automatically with a prescribed accuracy, even for t arbitrarily close to the singularity.

1. Introduction In a variety of controlled-potential transient electroanalytical experiments [1] one expects theoretically a singularity of the currenttime response i(t) at the initial time moment t = 0. Such temporal singularities typically occur in the presence of reversible (Nernstian) heterogeneous charge transfer reactions, when there is a discontinuity at t = 0, between initial and boundary conditions accompanying partial differential equations that describe a given experiment. Hence, the singularities are most pronounced in the case of potential step chronoamperometric experiments [1]. In other electroanalytical experiments, such as, in particular, the popular linear potential sweep or cyclic voltammetric experiments [1], the singularities are also frequently present formally, even though they tend to be ignored by modellers, in cases when they are unimportant for the analysis of experimental data. Under conditions of purely diffusional transport in the absence of homogeneous reactions, the singular current-time responses behave asymptotically like i(t) ∼ t −1/2 when t → 0. This behaviour was revealed, probably for the first time, by Cottrell, in his seminal analysis of chronoamperometry [2]. Later on, the behaviour was proven theoretically in a general way for electrodes of arbitrary geometry [3,4], and confirmed in dozens of studies devoted to particular electrode

E-mail address: [email protected]. URL: http://www.cyf-kr.edu.pl/~nbbienia. http://dx.doi.org/10.1016/j.jelechem.2017.05.027 Received 18 April 2017; Received in revised form 15 May 2017; Accepted 16 May 2017 Available online 19 May 2017 1572-6657/ © 2017 Elsevier B.V. All rights reserved.

arrangements. There is also a number of theoretical papers (see, for example, Refs. [5–9]) showing an analogous behaviour for various homogeneous reaction-diffusion systems (although a general proof does not appear available). Further instances of the Cottrell-type singular transients were discovered in systems with migration-diffusion transport [10–13]. A reliable digital simulation [14] of singular current-time responses is a difficult matter. Conventional simulation methods by finite differences, finite elements, or other techniques for the direct numerical solution of partial differential equations, are widely used for this purpose. However, it is probably not commonly realised that it is mathematically not legal to use such methods at the singularity at t = 0. This is because most of these methods require, for convergence, an appropriate regularity of the solutions, which is absent at t = 0. As a consequence, such simulations are grossly inaccurate close to t = 0. Readers interested in the magnitudes of such errors can consult Fig. 2 in Ref. [15] or Fig. 4 in Ref. [16], from which it is seen that errors of the order of 10%, 100%, and even 1000%, are well possible. The concrete values of the initial errors depend on the simulation method, the selection of discretisation grids, and the geometry of the diffusion field (but the errors occur under all kinds of diffusion conditions). The fact that numerical results obtained in such simulations for t > 0 are not entirely useless, results from the relatively fast damping of the initial

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

above calculation of approximate solutions and their error estimates requires an availability of analytical formulae for the moment integrals of the kernel functions Kκ (t , τ ):

errors. However, the damping requires several (possibly many) discrete time steppings to be performed. Therefore, the use of nonuniform temporal and/or spatial grids, or adaptive grids (with temporal grid nodes concentrated near t = 0), can improve the simulation of the singular transients at t > 0 (see, for example, Refs. [17–19]), although it cannot reduce the errors at nodes closest to t = 0. In certain cases one can eliminate the singularities from numerical simulations by using analytical solutions for simplified systems [20,21], if only such analytical solutions are available. In the present paper we shall demonstrate that the above difficulties can be avoided by an appropriate use of the integral equation (IE) method [22]. We shall describe a new simulation approach to singular transients, employing a specially dedicated recent extension [23] of the adaptive IE solution technique described in Refs. [24–42]. The technique is based on the Huber method [43] for one-dimensional Volterra IEs (with either nonsingular or weakly singular kernels), and it provides automatically solutions possessing a desired accuracy. The approach described here applies to one-dimensional diffusion transport without homogeneous reactions (with the exception of one example of twodimensional diffusion). Further generalisation of the approach, to onedimensional diffusion coupled with first order homogeneous reactions, is possible, but it requires an additional work that exceeds the scope of this study. The latter generalisation is planned to be investigated and described separately. There have not been many former attempts to model singular transients by means of the IEs. Some examples, corresponding to the one-dimensional diffusion transport without homogeneous reactions, are available in Refs. [44–47]. In all these examples the authors tried to solve the IEs analytically, rather than numerically, possibly due in part to the lack of a suitable numerical method. However, analytical solutions may not always be obtainable. As we shall see, the present numerical method is a viable, powerful, and more general alternative to analytical solutions.

t

∫

with m = 0, 1, 2. Alternatively, in cases when such formulae do not exist, highly accurate approximants for the moment integrals must be developed (throughout this paper we use the phrase “highly accurate” to denote quantities, the relative error of which is comparable to, or even lower than the error of machine representation of double precision variables, which is about 10 −16 according to the IEEE 754 standard [48]). This is required if we want the convergence of the method to be dependent primarily on quadrature errors and not to be affected by an inaccurate representation of the moment integrals. Relevant approximants have been elaborated for several kernel functions often encountered in electroanalytical modelling [24–42]. The local linearisation of the solutions, and error estimation based on truncated Taylor expansions, imply that the continuity of U(t) for t ≥ 0 and the existence of at least two derivatives of U(t) with respect to t, are required by the method described in Refs. [24–42] (although the method was found to operate satisfactorily also when the solution is continuous but not differentiable at t = 0). However, the method of Refs. [24–42] could not calculate singular solutions. This deficiency was eliminated in the recent work [23], by assuming that U(t) can be decomposed into the sum:

∼ U (t ) = U (t ) + U (t ). (4) ∼ In Eq. (4) U (t ) = [c1 t −1/2, …, cM t −1/2]T represents singular solution components, with c = [c1, …, cM] T denoting a vector of unknown coefficients, and U (t ) represents nonsingular solution components (possessing at least two derivatives for all t ≥ 0, with a possible exception of t = 0, where U (t ) must be bounded). The decomposition (4) of U(t) implies an analogous decomposition of Y (t): ∼ Y (t ) = Y (t ) + Y (t ).

As was described in Refs. [24–42], the adaptive Huber method is designed for solving generally nonlinear systems of Volterra IEs, that can be written in the form:

t

(1)

Yj (t ) =

0

∫

Kκ (t , τ ) Uμ (τ ) dτ , (6)

0

∼ whereas Y (t ) can be expressed as the vector ∼ ∼ ∼ Y (t ) = [cμ (1) Z κ (1) (t ), …, cμ (I ) Z κ (I ) (t )]T ,

(7) ∼ in which Z κ (t ) (for κ = 1, …, L) are integrals characteristic of the various kernels:

t

∫

(5)

The elements of Y (t ) (for j = 1, …, I) are given by the formulae analogous to Eq. (2):

In Eq. (1) F(⋅) = [F1(⋅), …, FM(⋅)]T is a vector of M functions representing the individual IEs, 0 denotes the zero vector, U(t) = [U1(t), …, UM(t)]T is a vector of M unknown functions of t, and Y (t) = [Y 1(t), …, Y T I(t)] is a vector of I integrals:

Y j (t ) =

(3)

0

2. The adaptive Huber method for singular transients

F (t , U (t ), Y (t )) = 0 .

Kκ (t , τ ) τ m dτ

Kκ (t , τ ) Uμ (τ ) dτ (2)

t

with j = 1, …, I. In Eq. (2) Kκ (t , τ ) with κ = κ(j) = 1, …, L denotes one of L possible kernels, and Uμ(t) with μ = μ(j) = 1, …, M is one of the unknowns U(t) (the one associated with the jth integral). According to the method, solutions U(t) are approximated by linear splines over a grid of dynamically selected discrete nodes tn (n = 0, 1, …) along the t axis, with local grid step sizes hn = tn − tn −1. Consequently, integrals Y (t) are approximated by product-integration trapezium quadratures, with coefficients resulting from the integration of the splines. In this way one obtains systems of nonlinear algebraic equations, which are solved numerically for the approximations Un to U (tn) at every successive node tn. Apart from determining approximate solutions, the method calculates estimates of their local errors. The estimates are used for deciding whether a discrete solution obtained possesses a requested accuracy, or has to be recalculated using a reduced step size hn. They also serve for predicting subsequent step sizes. In the spirit of the so-called product-integration methods for IEs, the

∼ Z κ (t ) =

∫ 0

Kκ (t , τ ) τ −1/2 dτ. (8)

Under a number of rather straightforward assumptions [23] concerning the functions F(⋅), one can determine coefficients c, and discrete nonsingular solution components un (together with their error estimates) numerically from the IEs. One important assumption is that the explicit dependence of F(⋅) on U(t) must vanish in the limit of t → 0. In the electrochemical context this assumption can be satisfied by requiring that quasi-reversible and irreversible heterogenous reactions do not contribute in any way to the IEs when t → 0, which is consistent with the wisdom that for such reactions one does not expect temporal singularities of U(t). Coefficients c can be determined almost exactly, simultaneously with the initial discrete solutions u0 and u1 (throughout this paper, we use the phrase “almost exact” to denote quantities, the error of which consists only of machine errors resulting from the finite precision arithmetic on computers). Once c is known, the determination 41

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

kfun(x) ≡kplp(x) [22,37], for which

of subsequent discrete values u2 , u3 , … amounts to solving a modified IE system:

1/2  (s ) = 1−tanh (s ) . ϕ s1/2

∼ ∼ F (t , U (t )+U (t ), Y (t )+Y (t )) = 0, (9) ∼ ∼ in which both U (t ) and Y (t ) can be considered to be known almost exactly, provided that analytical expressions for (or highly accurate ∼ approximants to) Z κ (t ) are available. As U (t ) is nonsingular, the numerical solution of Eq. (9) proceeds as in Refs. [24–42]. In addition, owing to the negligible error of c, the error estimates for un are practically equivalent to the error estimates of un. In view of the above discussion, in the case of singular solutions highly accurate procedures for computing integrals (8) must be supplied to the method, in addition to the procedures for comput ing moment integrals (3). For a few simple kernels relevant to ∼ electrochemistry, analytical expressions for Z κ (t ) exist, and are reported in Ref. [23]. This refers to Kκ (t , τ ) = 1, Kκ (t , τ ) = (t −τ )−1/2 , Kκ (t , τ ) = exp[−λ (t − τ )], Kκ (t , τ ) = (t − τ )−1/2 exp[− λ (t − τ )], Kκ Kκ (t , τ ) = erex{[λ (t −τ )]1/2 }, (t , τ ) = erfc{[λ (t −τ )]1/2 }, Kκ and (t , τ ) = daw{[λ (t −τ )]1/2 }, where λ ≥ 0 is a parameter, erex(z) = exp

The sign in Eq. (10) is negative, and parameter λ is the same as in ∼ case C. An approximation to Z κ (t ) is needed. (E) For one-dimensional diffusion inside a sphere (e.g. inside a mercury drop electrode), kfun(x) ≡kins(x) [22,39], for which

 (s ) = ϕ

(z )erfc(z), and daw(z ) =

exp(−z 2 )

∫

exp(ξ 2 )

(15)

where K0 and K1 denote modified Bessel functions of the second kind and order 0 and 1, respectively. The sign in Eq. (10) is negative, and ∼ parameter λ is the same as in case B. An approximation to Z κ (t ) is needed. (G) For one-dimensional diffusion inside a cylinder (e.g. inside a tubular electrode without electrolyte flow), kfun(x) ≡kinc(x) [22,40], for which

dξ is the Dawson integral.

∼ Some of these analytical expressions for Z κ (t ) involve special functions (such as erfc(z), erex(z), daw(z), or Bessel functions), the calculation of which requires appropriate numerical libraries. However, for a number of other important kernels it is necessary to develop highly accurate approximants, since no library routines are available. This effort was undertaken in the present work, and the resulting approximants are described in the next section.

1/2 ⎡ ⎤  (s ) = 1 ⎢ I 0 (s ) −1⎥ , ϕ 1/2 1/2 s ⎣ I1(s ) ⎦

(16)

where I0 and I1 denote modified Bessel functions of the first kind and order 0 and 1, respectively. The sign in Eq. (10) is positive, and ∼ parameter λ is the same as in case B. An approximation to Z κ (t ) is needed. (H) Diffusion to a microband electrode is not one-, but twodimensional. However, when diffusion coefficients of the members of a redox couple are identical one can use [42] Eq. (10) with kfun(x) ≡kband(x), for which

∼ 3. Approximants for integrals Z κ (t ) Most of the kernels associated with one-dimensional pure diffusion are of convolution type, and can be written in the form [22]: (10)

where λ ≥ 0 is a certain parameter, kfun(x) is a certain function of x = λ(t − τ)1/2, and the sign ( + or −) depends on kfun(x). There  (s ) of the Laplace variable s, such that kfun(x) usually exists a function ϕ  (s ), taken can be defined in terms of the inverse Laplace transform of ϕ 2 between the domains of θ = x , and s:

⎧  (s ) = s−1/2 + 2 ⎨π ϕ ⎩ ⎪

⎪

∞

∑

m) 2 [A (2 0 (s )]

m =0

−1 Fek′2m (0, −s ) ⎫ ⎬ , Fek2m (0, −s ) ⎭ ⎪

⎪

(17)

where Fek2m is one of the so-called radial, modified Mathieu functions of the third kind, Fek′2m denotes the derivative of Fek2m with respect to m) denotes special expansion coefficients its first argument, and A (2 0 occurring in the theory of the Mathieu functions (see Ref. [42] for details and references). In the limit of large s (which corresponds to small x), the following asymptotic formula is more practical for  (s ) values [42]: calculating ϕ

 (s )}. kfun(x ) = kfun(θ1/2 ) = L −1 {ϕ

(11) ∼ Eq. (11) can be helpful in obtaining approximants to Z κ (t ), as is shown below. In the present study we focus on the following cases A-H of the electrode|electrolyte geometry, and related functions kfun(x): (A) In the simplest case of one-dimensional diffusion in a semiinfinite region adjacent to a planar macroelectrode, kfun(x) ≡ 0, and ∼ parameter λ is not needed. In this case Z κ (t ) is obtainable analytically [23]. (B) For one-dimensional diffusion in a semi-infinite region outside a spherical electrode, kfun(x) ≡erex(x), the sign in Eq. (10) is negative, and parameter λ = D1/2/r0, where D is a diffusion coefficient and r0 is ∼ the electrode radius [22]. Also in this case Z κ (t ) is obtainable analytically [23]. (C) For one-dimensional diffusion inside a finite electrolyte layer adjacent to a planar electrode, with an impermeable external boundary, kfun(x) ≡kpli(x) [22,37], for which

1/2  (s ) = coth (s )−1 . ϕ s1/2

(14)

1/2 ⎤ ⎡  (s ) = 1 ⎢1− K 0 (s ) ⎥ , ϕ s1/2 ⎣ K1(s1/2 ) ⎦

0

Kκ (t , τ ) = [π (t −τ )]−1/2 ± λ kfun[λ (t −τ )1/2],

1 1 − 1/2 . s1/2 coth (s1/2 )−1 s

The sign in Eq. (10) is positive, and parameter λ is the same as in case B. ∼ An approximation to Z κ (t ) is needed. (F) For one-dimensional diffusion in a semi-infinite region outside a cylindrical wire electrode, kfun(x) ≡kcylw(x) [22,33,34], for which

z

2

(13)

⎡ ⎤−1 1  (s ) = s−1/2 − ⎢s1/2+ 1 − ϕ s−3/4 exp(−4s1/2 ) ⎥ . 1/2 ⎣ ⎦ 4 32(2π )

(18)

The sign in Eq. (10) is negative, and parameter λ = 4D /a, where a is ∼ the width of the microband. An approximation to Z κ (t ) is needed. ∼ As the kernel components in Eq. (10) are additive, and integral Z κ (t ) −1/2 for the component [π(t − τ)] is analytically obtainable [23], we can ∼ focus on looking for approximants to integrals Z κ (t ) corresponding to the kernel components 1/2

Kκ (t , τ ) = kfun[λ (t −τ )1/2]

(19)

with kfun(x) ≡ kpli(x), kplp(x), kins(x), kcylw(x), kinc(x), and kband (x). By taking into account Eqs. (8) and (11) it is straightforward to show that for such kernels

(12)

t

The sign in Eq. (10) is positive, and parameter λ = D /l, where l is the ∼ thickness of the electrolyte layer. An approximation to Z κ (t ) is needed. (D) For one-dimensional diffusion inside a finite electrolyte layer adjacent to a planar electrode, with a permeable external boundary, 1/2

∼ Z κ (t ) =

∫ 0

kfun[λ (t −τ )1/2 ] τ −1/2 dτ =

2 kfun sqrt (λt1/2 ), λ

(20)

where the function denoted by kfunsqrt(x), with x = θ

1/2

42

= λt

1/2

, can be

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

kfunsqrt(x) values was employed, to determine local polynomial approximations having the desired accuracy (see Subsections 3.1–3.6 for specific formulae). The present author's “C++” translation of the discrete minimax code of Schmitt [54] was used for the fitting. The lowx, large-x, and intermediate-x approximations were combined together in the form of C++ procedures serving for an accurate and efficient computation of the various functions kfunsqrt(x). The procedures were later combined with the code of the adaptive Huber method. As the amount of coefficients needed by all these procedures is considerable, we do not present these coefficients in the present paper, but we provide them in the Supplementary Material, together with the C++ code of the procedures. Intervals of x in which particular sets of coefficients are to be used, are also given. This part of the calculations, as well as testing simulations described below in Section 5, were accomplished using the Bloodshed/Orwell Dev-C++ 5.7.1 environment [55,56], equipped with the TDM GCC 4.8.1 compiler (a 32 bit release). The resulting programs were run as sequential, 32-bit console applications on a laptop computer with an Intel Centrino 2 processor (2.4 GHz), operating under Windows Vista. All calculations were done in extended precision (18–19 significant digits). Special function erfc, present in some of the equations in Subsections 3.1–3.6, was computed by using the present author translation (into C++) of the relevant routine from the publicly available (from netlib [57]) package SPECFUN of Cody [58]. Out of the approximants developed, the least accurate is the one for kbandsqrt(x), with the largest relative error of about 1.62 × 10 −16. The dependences on x, of the relative errors of the present approximants, are very similar to those previously obtained for a variety of approximants to kernels, their moment integrals, and other functions [33,37-40,42,59,60]. In the following Subsections 3.1–3.6 we provide detailed formulae adopted in highly accurate approximants to functions kfunsqrt(x).

expressed by any of the following alternative three formulae:

kfun sqrt (x ) =

kfun sqrt (x ) =

1 2

θ

∫

kfun[(θ − ϑ)1/2]ϑ

−1/2dϑ,

(21)

0

⎧ϕ  (s ) ⎫ π 1/2 L −1 ⎨ 1/2 ⎬, 2 ⎩ s ⎭

(22)

or x2

kfun sqrt (x ) =

∫

kfun[(x 2− ξ 2 )1/2] dξ. (23)

0

Hence, we need approximants to kfunsqrt(x). In order to construct such approximants, we used Eq. (22) to tabulate functions kfunsqrt(x) over large intervals of x ≥ 0, and with an accuracy of at least 19 significant digits. The GWR[...] procedure of Valkó and Abate [49–51], for the numerical Laplace transform inversion by the Gaver-Wynn-rho method, was run under the control of MATHEMATICA [52]. For computing the Mathieu functions, the package “Mathieu”, of Blose et al. [53] was chosen. These calculations were done on an Hewlett-Packard Apollo 8000 Gen9 cluster involving 53568 computational cores based on Intel Xeon E5-2680v3 processors, under Linux CentOS 7 operating system. Fig. 1 compares plots of the tabulated functions kplisqrt(x), kplpsqrt(x), kinssqrt(x), kcylwsqrt(x), kincsqrt(x), and kbandsqrt(x). As can be seen, functions kplisqrt(x), kinssqrt(x), and kincsqrt(x) grow unboundedly with increasing x. In contrast, functions kplpsqrt(x), kcylwsqrt(x), and kbandsqrt(x) grow up to a common limiting value (equal π1/2/2, see equations in Subsections 3.2, 3.4 and 3.6 below), when x →∞. Another observation, expected from the properties of the functions kfun(x) [22], is that the difference between kplisqrt(x) and kplpsqrt(x) vanishes for x → 0. Similarly, the difference between kcylwsqrt(x) and kincsqrt(x) vanishes for x → 0. For all functions kfunsqrt(x) investigated, series expansions or other analytical approximations can be derived (see Subsections 3.1–3.6), that are valid in the limits of either small, or large x. Comparison with the tabulated kfunsqrt(x) values gave domains of x where these approximations are sufficiently accurate. Relative errors not exceeding about 10 −16 were demanded from the approximations, in agreement with the adopted meaning of “high accuracy” (cf. Section 2). In the intermediate intervals of x, in which neither the low-x nor large-x approximations were satisfactory, minimax fitting to the tabulated

3.1. Approximations to kplisqrt(x) Function kpli(x) has two alternative representations [22,37]:

kpli(x ) =

∞

2 π 1/2x

∑ ν =0

⎡ ⎛ ν+1 ⎞2 ⎤ ⎟ ⎥, exp ⎢−⎜ ⎣ ⎝ x ⎠⎦

(24)

or

kpli(x ) = 1 −

∞

1 π 1/2x

+ 2 ∑ exp [−(ν+1)2π 2x 2]. ν =0

(25)

Setting Eq. (24) into Eq. (21) gives:

2.0

∞

kinssqrt

kplisqrt

kplisqrt (x ) = π 1/2

ν =0

kincsqrt

kfunsqrt(x )

kplisqrt (x ) = x −

π 1/2 2 + 2 π

∞

∑ ν =0

daw[(ν+1) πx ] . ν+1

(27)

By truncating the series (26) to several terms, one obtains a good approximation to kplisqrt(x), valid for small x. Eq. (27) is not practical for calculations, since a very large number of terms would have to be summed. However, for sufficiently large x one can modify Eq. (27) by replacing daw(z) by its standard asymptotic expansion:

kplpsqrt kcylwsqrt

0.5

(26)

whereas setting Eq. (25) into Eq. (21) gives:

1.5

1.0

⎛ ν+1 ⎞ ⎟, x ⎠

∑ erfc ⎜⎝

kbandsqrt

∞

daw(z ) ≈

∑ σ =0

(2σ −1) !! 1 . 2σ +1 z 2σ +1

(28)

This gives (after changing the summation order):

0.0

-3

-2

-1

0

1

2

3

4

kplisqrt (x ) ≈ x −

log(x )

π 1/2 + 2

with coefficients

Fig. 1. Plots of various functions kfunsqrt(x).

43

∞

∑ σ =0

aσ , x 2σ +1

(29)

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

aσ =

∞

(2σ −1) !! 2σ π 2σ +2

∑

(ν+1)−2(σ +1) =

ν =0

(2σ −1) !! ζ (2σ + 2), 2σ π 2σ +2

kinssqrt (x ) ≈ (30)

∞

kinssqrt (x ) ≈

∑

aσ x σ ,

∑

aσ x σ ,

(41)

σ =0

where a0 = 0, a1 = 1, a2 = π /2, a3 = 2/3, a4 = π /4, etc. By truncating the series (41) to several terms, one obtains a good approximation to kinssqrt(x), valid for small x. Numerical values of the first few coefficients aσ of the series (41) are listed in Table A.9. Another expression for kins(x), this time valid for all x, is [39]: 1/2

σmax

(31)

σ =0

(40)

or, after expanding erex(−x) into the standard power series:

∞

where ζ (α ) = ∑ν =1 ν−α (for α = 2, 3, …) is the Riemann Zeta function. By truncating the series (29) to several terms, one obtains a good approximation to kplisqrt(x), valid for large x. Numerical values of the first few coefficients aσ of the series (29) are listed in Table A.4 (in Supplementary Material). In the intermediate intervals of x we use the polynomial approximations

kplisqrt (x ) ≈

π 1/2 [erex(−x )−1], 2

1/2

∞

1 + 2 ∑ exp (−bν2 x 2 ), π 1/2x ν =1

kins(x ) = 3 −

with a suitably chosen σmax. Numerical values of the coefficients aσ of the polynomial (31) are listed in Tables A.1–A.3.

(42)

where bν (for ν = 1, 2, …) are positive roots of the equation tan(b) = b. By setting Eq. (42) into Eq. (21) one obtains:

3.2. Approximations to kplpsqrt(x)

∞

Function kplp(x) has two alternative representations [22,37]:

kplp(x ) =

⎡ ⎛ ν+1 ⎞2 ⎤ ⎟ ⎥, (−1)ν exp ⎢−⎜ ⎣ ⎝ x ⎠⎦

∞

2 π 1/2x

∑ ν =0

π 1/2 daw(bν x ) +2∑ . 2 bν ν =1

kinssqrt (x ) = 3x −

Combination of Eqs. (43) and (28) yields: (32)

π 1/2 + 2

kinssqrt (x ) ≈ 3x −

or ∞

kplp(x ) =

1 −2∑ π 1/2x ν =0

⎡ ⎛ 1 ⎞2 ⎤ exp ⎢−⎜ν+ ⎟ π 2x 2⎥. ⎝ ⎠ 2 ⎣ ⎦

aσ =

⎛ ν+1 ⎞ ⎟, x ⎠

∑ (−1)ν erfc ⎜⎝

kplpsqrt (x ) = π 1/2

ν =0

π 1/2 2 − 2 π

1 daw ⎡⎣ (ν+ 2 ) πx ⎤⎦

∞

∑

1

.

ν+ 2

ν =0

(35)

(2σ −1) !! 2σ

∞

∑ σ =0

−

bν−2(σ +1).

ν =1

x

(36)

with coefficients

aσ = =

(2σ − 1) ! ! 2σ π 2σ +2 (2σ − 1) ! ! 2σ π 2σ +2

∞

1 −2(σ +1)

( )

∑ν =0 ν+ 2

(22σ +2 −1) ζ (2σ +2).

(37)

∑

aσ x −σ ,

37 118 z 5− 197071875 z 7+…. 3031875

(46)

∑

aσ x σ ,

(47)

(for smaller x), and σmax

kinssqrt (x ) ≈ 3x +

∑

aσ x −σ ,

σ =0

(48)

(for larger x), with a suitably chosen σmax. Numerical values of the coefficients aσ of the polynomials in Eqs. (47) and (48) are listed in Tables A.10–A.36.

(38)

with a suitably chosen σmax. Numerical values of the coefficients aσ of the polynomial (38) are listed in Tables A.5–A.7.

3.4. Approximations to kcylwsqrt(x)

3.3. Approximations to kinssqrt(x)

Function kcylw(x) has the following asymptotic series expansion for x → 0 [33]:

Function kins(x) has an approximation valid for small x [39]:

kins(x ) ≈ exp(x 2 )erfc(−x ).

2

σ =0

σmax σ =0

1

σmax

kinssqrt (x ) ≈

By truncating the series (36) to several terms, one obtains a good approximation to kplpsqrt(x), valid for large x. Numerical values of the first few coefficients aσ of the series (36) are listed in Table A.8. In the intermediate intervals of x we use the polynomial approximations

kplpsqrt (x ) ≈

(45)

Hence, s1 = 1/10, s2 = 1/350, s3 = 1/7875, s4 = 37/6063750, s5 = 59/ 197071875, etc. (see also Ref. [62], p. 269). By truncating the series (44) to several terms, one obtains a good approximation to kinssqrt(x), valid for large x. Numerical values of the first few coefficients aσ of the series (44) are listed in Table A.37. In the intermediate intervals of x we use the polynomial approximations

aσ

, 2σ +1

(44)

∞

∑

1

π 1/2 − 2

aσ , x 2σ +1

f (z )≈3z−2 − 5 z−1− 175 z− 7875 z 3

By truncating the series (34) to several terms, one obtains a good approximation to kplpsqrt(x), valid for small x. Combining Eq. (35) with Eq. (28) yields, in turn:

kplpsqrt (x ) ≈

σ =0

The infinite sums of the inverse even powers of bν, present in Eq. (45), are known exactly. According to Pichereau [61], by using the calculus ∞ of residues, it can be proven that the sum sα = ∑ν =1 bν−2α (for α = 1, 2, 3, …) is equal to −1/2 times the coefficient standing at the successive powers of z (beginning with the z −1 power) in the Laurent series expansion for the function f(z) = [tan(z) − z] −1 + z −1, around z = 0. The expansion (as obtained using MATHEMATICA [52]) is:

(34)

whereas setting Eq. (33) into Eq. (21) gives:

kplpsqrt (x ) =

∞

∑

with coefficients (33)

Setting Eq. (32) into Eq. (21) gives: ∞

(43)

kcylw(x ) ≈

(39)

Setting Eq. (39) into Eq. (21) gives:

1 3 3 21 3 27 4 − x + x2 − x + x −… 2 4π 1/2 8 32π 1/2 64

Setting Eq. (49) into Eq. (21) gives: 44

(49)

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

∞

∞

∑

kcylwsqrt (x ) ≈

aσ x σ ,

(50)

σ =0

where a0 = 0, a1 = 1/2, a2 = −3π /16, a3 = 1/4, a4 = −63π /512, etc. By truncating the series (50) to several terms, one obtains a good approximation to kcylwsqrt(x), valid for small x. Numerical values of the first few coefficients aσ of the series (50) are listed in Table A.38. For large x, an alternative approximation is obtainable from Eqs. (15) and (22), by noting that for s → 0 (which corresponds to large x) there is K1(s1/2) ≈ s −1/2. This gives the inverse Laplace transform for kcylwsqrt(x): 1/2

1/2

⎧ 1 K (s1/2 ) ⎫ kcylwsqrt (x ) ≈ L −1 ⎨ − 0 1/2 ⎬, ⎩s ⎭ s 2

π 1/2 ⎧ ⎨1− 2 ⎩

aσ =

∑

⎭

(61)

The infinite sums of the inverse even powers of bν, present in Eq. (61), are known exactly. According to the properties of the Rayleigh special ∞ function [63–65], if we denote sα = ∑ν =1 bν−2α (for α = 1, 2, 3, …), then

(52)

s1 = 1/8, s2 = 1/192, and sα = α + 1 ∑γ =1 sγ sα − γ for α = 3, 4, 5, …. By truncating the series (60) to several terms, one obtains a good approximation to kincsqrt(x), valid for large x. Numerical values of the first few coefficients aσ of the series (60) are listed in Table A88. In the intermediate intervals of x we use the polynomial approximations

α −1

σmax

kincsqrt (x ) ≈

∑

aσ x σ ,

(62)

σ =0

(for smaller x), and σmax

(53)

∑

aσ x −σ ,

(63)

σ =0

(for larger x), with a suitably chosen σmax. Numerical values of the coefficients aσ of the polynomials in Eqs. (62) and (63) are listed in Tables A.69–A.87.

σmax

aσ x σ ,

bν −2(σ +1).

(51)

In the intermediate intervals of x we use the polynomial approximations

∑

(60)

ν =1

1

exp[−(8x 2 )−1][ ln(2)− γE+ln(8x 2 )] ⎫ ⎬. 2π 1/2x

(54)

σ =0

σ =0

aσ , x 2σ +1

∞

(2σ −1) !! 2σ

kincsqrt (x ) ≈ 2x +

kcylwsqrt (x ) ≈

∞

∑

with coefficients

In the interval of x where Eq. (52) is a valid approximation, the factor K0 [(8x2) −1] can be replaced by the first term of the asymptotic expansion of K0(z) for z → 0, which is K0(z) ≈ ln(2) − γE − ln(z), with γE denoting the Euler gamma constant. Consequently, we obtain the final formula for the large-x approximation:

kcylwsqrt (x ) ≈

π 1/2 + 2

kincsqrt (x ) ≈ 2x −

which can be derived analytically:

π 1/2 ⎧ exp[−(8x 2 )−1]K 0 [(8x 2 )−1] ⎫ ⎨1− ⎬. ⎭ 2 ⎩ 2π 1/2x

(59)

Combination of Eqs. (59) and (28) yields:

π 1/2

kcylwsqrt (x ) ≈

π 1/2 daw(bν x ) +2∑ . 2 bν ν =1

kincsqrt (x ) = 2x −

(for smaller x), and

3.6. Approximations to kbandsqrt(x)

σmax

∑

kcylwsqrt (x ) ≈

aσ x −σ ,

Function kband(x) has an approximation valid for small x [42]:

(55)

σ =0

⎛x⎞ 1 erex ⎜ ⎟ . ⎝4⎠ 4

(for larger x), with a suitably chosen σmax. Numerical values of the coefficients aσ of the polynomials in Eqs. (54) and (55) are listed in Tables A.39–A.67.

kband(x ) ≈

3.5. Approximations to kincsqrt(x)

kband sqrt (x ) ≈

Function kinc(x) has the following asymptotic series expansion for x → 0 [40]:

or, after expanding erex(x/4) into the standard power series:

kinc(x ) ≈

3 3 21 3 27 4 1 + x + x2 + x + x +… 2 4π 1/2 8 32π 1/2 64

Setting Eq. (64) into Eq. (21) gives:

aσ x σ ,

(56)

1 π 1/2x

kband sqrt (x ) ≈

ν =1

(66) 1/2

π 1/2 −1 ⎧ 1 4K 0 (s1/2 )I 0 (s1/2 ) ⎫ ⎬. L ⎨ − ⎩s ⎭ πs1/2 2

(67)

For s → 0 (which corresponds to large x) there is I0(s inversion of the Laplace transform then gives:

) ≈ 1. Analytical

1/2

∞

+ 2 ∑ exp (−bν2 x 2 ),

aσ x σ ,

where a0 = 0, a1 = 1/4, a2 = −π /32, a3 = 1/96, a4 = −π /1024, etc. By truncating the series (66) to several terms, one obtains a good approximation to kbandsqrt(x), valid for small x. Numerical values of the first few coefficients aσ of the series (66) are listed in Table A.89. For large x, an alternative approximation is obtainable from Eqs. (17) and (22), and from the asymptotic formula (58) in Ref. [42], which give together:

where a0 = 0, a1 = 1/2, a2 = 3π1/2/16, a3 = 1/4, a4 = 63π1/2/512, etc. By truncating the series (57) to several terms, one obtains a good approximation to kincsqrt(x), valid for small x. Numerical values of the first few coefficients aσ of the series (57) are listed in Table A.68. Another expression for kinc(x), valid for all x, is [40]:

kinc(x ) = 2 −

∑

1/2

(57)

σ =0

(65)

σ =0

∞

∑

⎛ x ⎞⎤ π 1/2 ⎡ ⎢1-erex ⎜⎝ ⎟⎠ ⎥⎦ , 2 ⎣ 4

∞

kband sqrt (x ) ≈

(the absolute values of the expansion coefficients are identical as in Eq. (49), but all coefficients are now positive). Setting Eq. (56) into Eq. (21) gives:

kincsqrt (x ) ≈

(64)

kband sqrt (x ) ≈

(58)

π 1/2 ⎧ 2 exp[−(8x 2 )−1]K 0 [(8x 2 )−1] ⎫ ⎨1− ⎬. ⎭ π 3/2x 2 ⎩

(68)

Similar to Eq. (52), the factor K0 [(8x2) −1] can be replaced by the first term of its asymptotic expansion, giving the final result for the large-x approximation:

where bν (for ν = 1, 2, …) are positive roots of the equation J1(b) = 0, with J1 denoting the Bessel function of the first kind and order 1. By setting Eq. (58) into Eq. (21) one obtains: 45

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

π 1/2 ⎧ 2 exp[−(8x 2 )−1][ ln(2)−γE+ln(8x 2 )] ⎫ ⎨1− ⎬. 2 ⎩ π 3/2x ⎭

kband sqrt (x ) ≈

i (t ) = (πt )−1/2 − λ kplp(λ t 1/2 ). (69)

In case D the solution is:

In the intermediate intervals of x we use the polynomial approximations

i (t ) = (πt )−1/2 + λ kpli(λ t 1/2 ).

∑

aσ x σ ,

(70)

σ =0

(for smaller x), and

i (t ) = (πt )−1/2 + λ icylw(λ t 1/2 ),

σmax

kband sqrt (x ) ≈

∑

aσ x −σ ,

(for larger x), with a suitably chosen σmax. Numerical values of the coefficients aσ of the polynomials in Eqs. (70) and (71) are listed in Tables A.90–A.113.

i (t ) =

In order to test the validity of the adaptive Huber method [24–42] with extension [23], for simulating singular electrochemical transients, we choose the following examples 1–3. 4.1. Example 1: Potential step chronoamperometry for an E reaction

(72)

iins (x ) = L −1 {s−1/2 [1−coth (s1/2 )+s−1/2]},

(where only species O is initially present at concentration c ), subject to a cathodic potential step at t = 0 (an anodic step for a reverse reaction is described by identical equations, after appropriate sign changes). Assuming pure diffusion in the absence of homogeneous reactions, the concentration c† of species O (at the electrode) is related to its flux component J⊥ perpendicular to the electrode, by the formula (5.6) in Ref. [22]:

∫

with the Laplace transform between the domains of θ = x and s. For large s, ε = exp( −2s1/2) is small, so that ∞

1 − coth (s1/2 ) = −2

ε = −2 ∑ ε σ 1−ε σ =1

(83)

is a rapidly convergent series. Hence,

Kκ (t , τ ) [ ±

J ⊥]

∞ ⎧1 exp(−2σs1/2 ) ⎫ ⎬. iins (x ) = L −1 ⎨ −2 ∑ s s1/2 ⎩ ⎭ σ =1

dτ. (73)

0

We test all cases A–H (see Section 3) of the kernel Kκ (t , τ ) given by Eq. (10). The minus sign before J⊥ in Eq. (73) must be taken in cases E and G, the plus sign must be taken in remaining cases. Under limiting current conditions c† = 0 in Eq. (73) which then becomes an IE for J⊥, or for the related Faradaic current. It is convenient to express this IE using dimensionless variables t = t / tstep , 1/2 τ = τ / tstep , and λ = λtstep , where tstep is the duration of the experiment. can be defined as The dimensionless current i (t ) i (t ) = ± i (t )/[nFAc★ (D /tstep )1/2 ], where i(t) is the dimensional current, A is the electrode area, F is the Faraday constant, and D is the diffusion coefficient of O. The plus sign before i(t) must be taken in cases E and G, the minus sign in remaining cases. Consequently, by combining Eqs. (10) and (73) one obtains the first kind Volterra IE for the dimensionless current response: t

∫

(82) 2

t

+

(81)

where

★

=

(80)

i (t ) = (πt )−1/2 − λ iins (λ t 1/2 ),

In this example we consider the charge transfer reaction:

c★

⎡ ⎛ ⎞2 ⎤ λ λ iband ⎢ ⎜⎜ ⎟⎟ t ⎥ , ⎢⎣ ⎝ 4 ⎠ ⎥⎦ 4

where the expression for iband (t ) is defined by Eq. (7) in Ref. [60] (note that different dimensionless variables were used in Ref. [60], hence the re-scaling of t by (λ /4)2 is necessary). Similar approximants to the solutions for cases E and G have not been published thus far. Therefore, in the present study such approximants have been elaborated, and are described below. In case E, by applying the Laplace transform to Eq. (74) with kfun(x) ≡kins(x), one obtains:

4. Examples

O + ne− ⇋ R

(79)

where icylw(x) is the function defined by Eq. (11) in Ref. [59]. In case H we use

(71)

σ =0

c†

(78)

Highly accurate approximants to the solutions for cases F and H were elaborated in Refs. [59] and [60], respectively. These approximants have been used here. In particular, in case F we use:

σmax

kband sqrt (x ) ≈

(77)

0

∫ 0

⎪

⎪

⎪

(84)

Inverting the series analytically, term by term, gives:

⎡∞ ⎛ σ2 ⎞⎤ iins (x ) = 1 − 2 ⎢ ∑ exp ⎜− 2 ⎟ ⎥ (π 1/2x )−1. ⎝ x ⎠ ⎥⎦ ⎢⎣ σ =1

(85)

By truncating the series (85) to several terms, one obtains a good approximation to iins(x), valid for small x. For small s, in turn, s −1/2 [1 − coth(s1/2) + s −1/2] ≈ s −1/2, so that, after inversion, one obtains an approximation valid for large x:

iins (x ) ≈ (π 1/2x )−1.

(86)

Eqs. (85) and (86) suffice to compute iins(x) with a high accuracy, for any x ≥ 0. The relevant C++ code for calculating iins(x) is provided in the Supplementary Material. In case G one finds, in a similar way, that

t

[π (t −τ )]−1/2 i (τ ) dτ ± λ

⎪

kfun[λ (t −τ )1/2 ] i (τ ) dτ − 1 = 0.

i (t ) = (πt )−1/2 − λ iinc (λ t 1/2 ),

(74)

(87)

where

The advantage of this example is that analytical solutions, or highly accurate approximations to the solutions of Eq. (74) are available, and can be used as a reference for quantitative determinations of simulation errors. In case A, the well known solution [2] is:

⎧ 1 ⎡ I (s1/2 ) ⎤ ⎫ iinc (x ) = L −1 ⎨ 1/2 ⎢1− 1 1/2 ⎥ ⎬. ⎩ s ⎣ I 0 (s ) ⎦ ⎭

i (t ) = (πt )−1/2 .

(75)

For large s, from asymptotic expansions of the Bessel functions one obtains:

(76)

s1/2

In case B, the solution is:

i (t ) =

(πt )−1/2

+ λ.

1 ⎡ I1(s1/2 ) ⎤ 1 1 1 25 −5/2 s + …. ⎢1− ⎥ ≈ s−1 + s−3/2 + s−2 + ⎣ I 0 (s1/2 ) ⎦ 2 8 8 128

In case C the solution is:

Term by term Laplace inversion gives: 46

(88)

(89)

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz ∞

iinc (x ) ≈

∑

diffusion coefficient ratio), S (t , ts ) is the sawtooth function:

aσ x σ ,

(90)

σ =0

⎧t S (t , ts ) = ⎨ ⎩ 2ts−t

where a0 = 1/2, a1 = 1/(4π1/2), a2 = 1/8, a3 = 25/(96π1/2), a4 = 13/ 64, etc. By truncating the series (90) to several terms, one obtains a good approximation to iinc(x), valid for small x. Numerical values of the first few coefficients aσ of the series (90) are listed in Table A.114. For small s, in turn, s −1/2 [1 −I1(s1/2)/I0(s1/2)] ≈ s −1/2, so that, after inversion, one obtains an approximation valid for large x:

iinc (x ) ≈ (π 1/2x )−1.

(91)

σmax

∑

aσ x σ ,

(92)

σ =0

(for smaller x), and σmax

iinc (x ) ≈

∑

aσ x −σ ,

(93)

σ =0

(for larger x), with a suitably chosen σmax. Numerical values of the coefficients aσ of the polynomials in Eqs. (92) and (93) are listed in Tables A.115–A.150. The relevant C++ code for calculating iinc(x) is also provided in the Supplementary Material. 4.2. Example 2: Cyclic voltammetry for an Erev reaction

4.3. Example 3: Potential step chronoamperometry for an ErevErev mechanism

In this example we consider again reaction (72), assuming its reversibility. We simulate cyclic voltammetry, paying special attention to the initial singularity of the current. As was noted in Section 1, such a singularity is usually ignored by the modellers (for example, large discrete time steps are used at t = 0 to overstep the singularity). However, this is not a preferred practice, since it introduces large and uncontrolled numerical errors into simulations. These errors can be large, since a theoretically infinite current is simulated by a merely finite current jump. They are uncontrolled, in the sense that it is not possible to estimate and modify their magnitudes by following standard formulae for the dependences of the truncation or discretisation errors (of the conventional simulation methods [14]) on integration steps. The formulae assume an appropriate regularity of the solution (usually a few terms of the Taylor expansion must exist), which is not the case at the singularity. In former tests of the adaptive Huber method, cyclic voltammetry for reaction (72) was simulated by subtracting the singular component of the current, already at the stage of formulating the IEs [29,31,32,34,36-38,40,42]. In this way only the nonsingular component of the current was simulated. The subtraction is no longer necessary when using the present extension of the method (which facilitates the formulation of the IEs, especially for complicated reaction mechanisms). By using standard normalisations of variables [66], the IE describing cyclic voltammetry for the reaction (72) is: t

∫

Examples 1 and 2 are indispensable for testing the adaptive Huber method, but they are described by single linear IEs. As was demonstrated in Ref. [23], the method is capable of handling much more complicated, especially nonlinear IEs and IE systems possessing singular solutions. Unfortunately finding similarly complicated examples of IEs describing electrochemical transient experiments, that would also have analytical (or other highly accurate) reference solutions, is difficult. We therefore restrict further tests to a linear IE system describing potential step chronoamperometry for an ErevErev mechanism at a planar electrode in a semiinfinite electrolyte domain. We consider a cathodic potential step for two reversible charge transfers:

A1 + n1e− ⇋ A2

A2 + n2

∫

0

e−

(97) (98)

⇋ A3

E10

kfun[λ (t −τ )1/2 ] χ (τ ) dτ − f (t ) = 0.

t

(γ12+ε1) ∫

0

(94)

0

t ψ1 (τ ) [π (t − τ )]1/2

dτ +γ12 ∫

t

In Eq. (94) χ (t ) is the dimensionless current function (in case A known as the Randles-Ševčik function), t = t (nFv )/(RT ) is the dimensionless time (often denoted by at in the literature [66]), τ is the dimensionless integration variable, and λ = λ [RT /(nFv )]1/2 , with v denoting the rate of the potential sweep. Function f (t ) is defined as

f (t ) = {1+exp [w−S (t , ts )]}−1 ,

E20 ,

and respectively. We assume generally having formal potentials different diffusion coefficients D1, D2 and D3 of species A1, A2 and A3, out of which only A1 is initially present. Similar to example 1, we introduce dimensionless variables t = t / tstep , τ = τ / tstep , and dimensionless fluxes ψ1 (t ) and ψ3 (t ) of species A1 and A3 at the electrode, defined as ψ1 (t ) = −J1⊥ /[c★ (D1/ tstep )1/2], ψ3 (t ) = −J3⊥ /[c★ (D1/ tstep )1/2 ]. By applying the usual rules [22] of converting the diffusional initial boundary value problems into IEs, we obtain the following system of first kind Volterra IEs to be solved for the unknowns ψ1 (t ) and ψ3 (t ):

t

(t − τ )−1/2 χ (τ ) dτ ± π 1/2λ

(96)

with ts denoting a dimensionless switching time. The singularity subtraction used in Refs. [29,31,32,34,36-38,40,42] consisted in replacing f (t ) by f (t ) − f (0) in Eq. (94). Since for t ≈ 0 the function f (t ) can be approximated by f(0), the comparison of Eqs. (74) and (94) reveals that the initial part of χ (t ) should be well approximated by π −1/2f (0) g (t ), where the formulae for g (t ) are formally identical to the formulae for i (t ) in example 1 (cf. Eqs. (75)–(81), and (87)), taking into account somewhat differently defined t and λ . This initial singular part of the current function decays quickly (on the time scale of a voltammetric experiment), so that at later times there remains only a χ (t ) response to the potential sweep perturbation. For the latter response independently obtained accurate reference values are not immediately available. We therefore verify long-time adaptive solutions only in cases A and H. In case A we compare adaptive solutions with the predictions of the series expansion algorithm of Myland and Oldham [67], which provides χ (t ) values accurate to about ± 5 × 10 −6. In case H we make use of the recently elaborated approximant [42] employing the QUADPACK code with the absolute error tolerance parameter equal zero, and the relative error tolerance parameter equal 10 −7.

In the intermediate intervals of x we use the polynomial approximations

iinc (x ) ≈

for t ≤ts , for t > ts

ε2 ∫ 0

ψ1 (τ ) [π (t − τ )]1/2

0 t

dτ +(γ23+ε2 ) ∫ 0

ψ 3 (τ ) [π (t − τ )]1/2 ψ 3 (τ ) [π (t − τ )]1/2

⎫ dτ −ε1=0 ⎪ ⎪ ⎬. ⎪ dτ =0 ⎪ ⎭

(99)

In Eq. (99) γ12 = (D1/D2)1/2, γ23 = (D2/D3)1/2, ε1 = exp(n1w1), ε2 = exp (n2w2), with w1 = F (E10−E )/(RT ) and w2 = F (E20−E )/(RT ), where E is the potential step value. The total dimensionless current can be defined as

(95)

where w is a parameter dependent on the difference between the starting potential and the formal potential of reaction (72) (in the case of different diffusion coefficients of O and R it also depends on the

i (t ) = n1ψ1 (t ) − n2 ψ3 (t ). System (99) can be solved analytically with the result 47

(100)

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

⎫ (πt )−1/2 ⎪ ⎬, ε ε ψ3 (t )=− ε (γ + ε1 )2+ γ γ (πt )−1/2 ⎪ ⎭ 1 23 2 12 23 ψ1 (t )= ε

ε1 (γ23 + ε 2)

1 (γ23 + ε 2 ) + γ12 γ23

-3

(101)

from which

log( max|error|)

-4

n1ε1 (γ23+ε2 )+n2 ε1 ε2 i (t ) = (πt )−1/2 . ε1 (γ23+ε2 )+γ12 γ23

(102)

Although in example 3 only planar diffusion is considered, it should be stressed that there is no difficulty to take into account other electrode geometries, by including appropriate functions kfun(x) into the IE system (99). We do not do this here only because of difficulties with obtaining accurate reference solutions in such situations.

-5 -6 -7

5. Results and discussion

-8

Simulation results obtained for cases A–H in example 1 can be divided into three groups. The first group contains results for cases A and B, for which current transients i (t ) are almost exact, independently of the values of the model parameter λ . No grid adaptation is necessary. This happens because according to Eqs. (75) and (76), the nonsingular contribution to the solution i (t ) (cf. the discussion of Eq. (4) in Section 2) is then a constant (and hence linear) function of t , and in such cases the Huber method yields theoretically exact solutions. We do not show such results, because more interesting are analogous results for example 3 (see later). The second group contains results for cases C–G, for which the nonsingular contribution to i (t ) is a non-negligible nonlinear function of t . Consequently, the results are contaminated by discretisation errors, and a grid adaptation is needed to control these errors. The grid adaptation is found to operate very satisfactorily, similar to all previously studied example IEs, to which the adaptive Huber method was applied [23–42]. Fig. 2 shows typical simulated transients corresponding to cases C–G, assuming a moderately large value of λ = 1. In general, the grid nodes are concentrated in places requiring an increased resolution of the simulated current, although in some instances we observe a relatively low density of the grid nodes in the initial parts of the transients (see the next paragraph for explanation). The errors of the simulated values are always smaller than or compar-

-8

– – i(t) 4

2

0 0.6

0.8

-4

-3

able to the requested error tolerance tol, within the entire interval of t , including t arbitrarily close to the singularity point. This is seen in Fig. 3, which shows how the maximum absolute errors of the transients depend on tol in cases C-G. Such a successful error control would not be achievable by conventional digital simulation methods based on directly solving partial differential equations (cf. Section 1). Fig. 3 demonstrates that it is possible to obtain a requested error level, down to about 10 −8. Lower errors may also be achievable, but the computational cost increases with decreasing tol. The temporal grid adaptation usually begins with reducing the starting integration step hstart supplied to the method, down to a value h1 at which the initial solution values satisfy the requested accuracy criteria [23–27]. In the present simulations of singular transients this implies that with decreasing tol the first accepted discrete solution approaches the singularity (i.e. the solution value grows in magnitude). Obtaining a controlled error magnitude in such circumstances appears to be a considerable success. The third group of results involves results that show an intermediate behaviour between the former two groups. When λ = 1 such a behaviour is observed for case H. At relatively small t , or when the error tolerance tol is not very small, the transients are simulated very accurately (compared to tol), which indicates that the nonsingular solution component is negligible in comparison with the singular component in Eq. (4). Only after a sufficiently long time t , and/or when tol is chosen very small, the grid adaptation comes into action. This behaviour is illustrated by Fig. 4, which shows results obtained in case H assuming a rather demanding tol = 10 −10. As can be seen, there is a need for grid adaptation only for t greater than about 0.2. Note that this is an opposite behaviour to that observed in conventional simulations by directly solving partial differential equations, mentioned in Section 1, where the errors are biggest initially, and decrease with time. As a consequence, in case H it was possible to obtain much more accurate solutions than in cases C-G, without excessive computational costs. In particular, for λ = 1, hstart = 0.01, and hmax = 0.1, errors as low as 10 −12 could easily be obtained, and the proportionality of the errors

6

0.4

-5

Fig. 3. Dependence of the decimal logarithm of the maximum true absolute error max|error| on the decimal logarithm of the error tolerance parameter tol, observed in example 1, cases C (▴ ), D (▾ ), E(♦ ), F(□ ), and G(∘), assuming λ = 1. Symbols are joined by solid lines. Method parameters are: starting step hstart = 0.01, maximum step hmax = 0.1. The values of tol are: 10 −3, 5 × 10 −4, 2 × 10 −4, 10 −4, 5 × 10 −5, 2 × 10 −5, 10 −5, 5 × 10 −6, 2 × 10 −6, 10 −6, 5 × 10 −7, 2 × 10 −7, 10 −7, 5 × 10 −8, 2 × 10 −8, and 10 −8 (from right to left). The dashed line represents the equality: max|error| = tol. Computational times (per simulated transient) needed by these simulations vary between about 0.01 s and 294 s (case C), 0.01 s and 285 s (case D), 0.01 s and 194 s (case E), 0.006 s and 62 s (case F), and 0.01 s and 363 s (case G).

8

0.2

-6

log(tol)

10

0

-7

1

– t Fig. 2. Adaptive solutions of Eq. (74), obtained in cases C (▴ ), D (▾ ), E(♦ ), F(□ ), and G(∘), assuming λ = 1, and method parameters: error tolerance tol = 10 −4, starting step hstart = 0.01, maximum step hmax = 0.1. Solid lines denote reference solutions calculated from Eqs. (77)–(79), (81) and (87).

48

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

10

10

a

a 8

6

6

– – i(t)

– – i(t)

8

4

2

2

0

0 0

0.2

0.4

– t

0.6

0.8

0

1 3

0.2

0.4

0.2

0.4

– t

0.6

0.8

1

0.6

0.8

1

b

b

0

2.5 2 – – log( i ( t ))

-0.2 error × 1010

4

-0.4

1.5 1

-0.6 0.5

-0.8 0

-1

-0.5

0

0.2

0.4

– t

0.6

0.8

1

0

– t

Fig. 5. Effect of parameter λ on the adaptive solution of Eq. (74) in cases C (a), and H (b), obtained with method parameters: error tolerance tol = 10 −4, starting step hstart = 0.01, maximum step hmax = 0.1. Notations in subfigure (a): λ = 0.1 (■), 1 (•), 2 (▴ ), 5 (▾ ), and 10 (♦ ). Notations in subfigure (b): λ = 0.1 (■), 1 (•), 10 (▴ ), 102 (▾ ), and 103 (♦ ). Solid lines denote reference solutions calculated from Eqs. (77) and (80).

Fig. 4. Adaptive solution of Eq. (74) (a), and its absolute errors (b), obtained in case H, assuming λ = 1, and method parameters: error tolerance tol = 10 −10, starting step hstart = 0.01, maximum step hmax = 0.1. Notations in subfigure “a”: (▵ ) - adaptive solutions; solid line - reference solution calculated from Eq. (80). Notations in subfigure “b”: ( ×, joined by dashed lines) - estimated errors; (▵ , joined by solid lines) - true errors. Computational time needed by this simulation equals about 0.2 s.

one can see that higher λ implies a greater number of grid nodes. Hence, depending on the value of λ , simulation results can belong to the second or the third group, in all cases C–H. Simulation results obtained for example 2 are shown in Fig. 6. The adaptive Huber method provides reliable solutions both in the long time scale (cf. Fig. 6a) and in the short time scale, just after the initial singularity (cf. Fig. 6b). At sufficiently small times t , the adaptive solutions corresponding to the various cases A–H are visually not distinguishable, reflecting mostly the t −1/2 -like decay of χ (t ) (at least when λ = 1, as is assumed in Fig. 6), in agreement with the short-time predictions indicated in Subsection 4.2. The differences between the voltammograms corresponding to the various cases A–H are clearly seen in the long time scale. Voltammograms obtained in cases A and H match well the long-time reference solutions indicated in Subsection 4.2. Since in example 3 the nonsingular components of the unknowns ψ1 (t ) and ψ3 (t ) are theoretically equal zero (cf. Eq. (101)), we should obtain almost exact simulation results, similar to cases A and B in example 1. Fig. 7 depicts some of the results obtained, clearly

to the tol parameter was fairly well obeyed for 10 −12 ≤ tol ≤ 10 −7 (with computational times, per simulated transient, varying between about 0.002 s and 187 s), whereas for tol > 10 −7 the errors stayed at the level of 10 −7. Obviously, it would not be possible to achieve such low errors without having highly accurate approximants for functions: kband(x), its moment integrals, and kbandsqrt(x). Depending on the value of λ , results obtained for other cases (C–G) may also enter the third group, especially when λ is small, and the associated contribution of the nonsingular solution component in Eq. (4) is small, too. In fact, a relatively sparse t-grid is observed in some cases also in Fig. 2, as was noted. Owing to the limited role played by the grid adaptation in the first and third groups of the results, one must use appropriately small values of hstart, and of the maximum step size hmax, if solution values densely tabulated along the t axis are of interest. The computational cost increases with increasing λ , independently of whether the solution decreases (as is in cases C, E, G) or increases (as is in cases D, F, H) with increasing λ . This is illustrated by Fig. 5, where 49

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

10

0.8

a

a 8

0.4

6

– χ(t )

– – i(t)

0.6

0.2

4

2

0

-0.2

0 0

10

20

30 – t

40

50

60

0

1.4

0.2

0.4

– t

0.6

0.8

1

2

b

b

1.2

1 1

0 error × 1011

– χ(t )

0.8 0.6 0.4

-1 -2

0.2

-3 0 0

0.2

0.4 0.6 – t × 105

0.8

1

-16

Fig. 6. Long-time (a) and short-time (b) behaviour of the adaptive solutions of Eq. (94), obtained in cases A (■), B (•), C (▴ ), D (▾ ), E(♦ ), F (□ ), G (∘), and H (▵ ), assuming λ = 1, w = 9, ts = 30 , 0 ≤ t ≤ 60 , and method parameters: error tolerance tol = 10 −4, starting step hstart = 10 −9, maximum step hmax = 1. Solid lines denote reference solutions indicated in Subsection 4.2.

-14

-12

-10

-8 -6 – log( t )

-4

-2

0

Fig. 7. Effect of parameter w2 on the simulated currents in example 3 (a), and corresponding true absolute errors of the current (b). The values of w2 are: 8 (■), 2 (•), 1 (▴ ), 0 (▾ ), −1 (♦ ), −2 (□ ), and −8 (∘). The remaining model parameter values are: n1 = n2 = 1, w1 = 8, γ12 = γ23 = 1. Method parameters are: error tolerance tol = 10 −4, starting step hstart = 10 −16, maximum step hmax = 0.02. Solid lines in subfigure “a” denote analytical predictions of Eq. (102). In subfigure “b” solid lines only join the data symbols.

demonstrating the perfect agreement between numerical and analytical solutions corresponding to various combinations of model parameters w1, w2, γ12 and γ23. In view of the fact that for IE systems the determination of coefficients c (cf. Eq. (4)) requires somewhat more elaborate computations than in examples 1 and 2 (see Ref. [23] for details), it is interesting to investigate how big the errors of the simulated currents are, depending on the step sizes used. By choosing a very small starting step hstart (such as hstart = 10 −16 in Fig. 7) we force several initial integration steps hn to be also comparably small, until they increase up to the maximum step size hmax requested. As can be seen in Fig. 7b, this causes the absolute errors to be larger initially (reaching about ± 4 × 10 −11) than at the subsequent steps (where the errors do not exceed about 10 −13). In contrast, relative errors increase somewhat with t , between about 10 −18 at low t , up to about 10 −13 at t approaching unity. This is again an opposite behaviour to that typical for conventional simulations based on direct solution of partial differential equations. Even the largest errors observed in Fig. 7b are much smaller than the smallest values of the error tolerances tol that are practical for simulations by the adaptive Huber method [24–42], which are usually

not smaller than 10 −8. It is also unlikely that errors as small as ± 4 × 10 −11 may be needed in practical applications of the simulation results to experimental data analysis. Hence, the accuracy achieved in example 3 should be sufficient for most of the applications of the method. 6. Conclusions The results of the simulations, performed in this work, demonstrate that by using the recent extension [23] of the adaptive Huber method [24–42] one can simulate singular transient electrochemical responses reliably and with a prescribed accuracy. Errors are well controlled, arbitrarily close to the singularity. Hence, the method is free from the usual problem associated with the conventional simulation methods (by directly solving partial differential equations), which are grossly inaccurate in the neighbourhood of the temporal singularity. Furthermore, the extended adaptive Huber method operates automati50

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

[22] L.K. Bieniasz, Modelling Electroanalytical Experiments by the Integral Equation Method, Springer, Berlin, 2015. [23] L.K. Bieniasz, An adaptive Huber method for nonlinear systems of Volterra integral equations with weakly singular kernels and solutions, J. Comput. Appl. Math. 323 (2017) 136–146. [24] L.K. Bieniasz, An adaptive Huber method with local error control, for the numerical solution of the first kind Abel integral equations, Computing 83 (2008) 25–39. [25] L.K. Bieniasz, Initialisation of the adaptive Huber method for solving the first kind Abel integral equation, Computing 83 (2008) 163–174. [26] L.K. Bieniasz, An adaptive Huber method for weakly singular second kind Volterra integral equations with non-linear dependencies between unknowns and their integrals, Computing 87 (2010) 35–54. [27] L.K. Bieniasz, An adaptive Huber method for non-linear systems of weakly singular second kind Volterra integral equations, Appl. Math. Comput. 217 (2011) 5622–5631. [28] L.K. Bieniasz, Extension of the adaptive Huber method for Volterra integral equations arising in electroanalytical chemistry, to convolution kernels exp [ −α(t − τ)]erex {[β(t − τ)]1/2} and exp[ −α(t − τ)]daw {[β(t − τ)]1/2}, J. Comput. Meth. Sci. Eng. 11 (2011) 323–338. [29] L.K. Bieniasz, Cyclic voltammetric current functions determined with a prescribed accuracy by the adaptive Huber method for Abel integral equations, Anal. Chem. 80 (2008) 9659–9665. [30] L.K. Bieniasz, Automatic simulation of cyclic voltammograms by the adaptive Huber method for weakly singular second kind Volterra integral equations, Electrochim. Acta 55 (2010) 721–728. [31] L.K. Bieniasz, Automatic simulation of cyclic voltammograms by the adaptive Huber method for systems of weakly singular Volterra integral equations, J. Electroanal. Chem. 642 (2010) 127–134. [32] L.K. Bieniasz, Extension of the adaptive Huber method for solving integral equations occurring in electroanalysis, onto kernel function representing fractional diffusion, Electroanalysis 23 (2011) 1506–1511. [33] L.K. Bieniasz, A highly accurate, inexpensive procedure for computing integral transformation kernel and its moment integrals for cylindrical wire electrodes, J. Electroanal. Chem. 661 (2011) 280–286. [34] L.K. Bieniasz, Automatic simulation of electrochemical transients at cylindrical wire electrodes, by the adaptive Huber method for Volterra integral equations, J. Electroanal. Chem. 662 (2011) 371–378. [35] L.K. Bieniasz, Automatic simulation of electrochemical transients by the adaptive Huber method for Volterra integral equations involving kernel terms exp [ −α(t − τ)]erex {[β(t − τ)]1/2} and exp[ −α(t − τ)]daw {[β(t − τ)]1/2}, J. Math. Chem. 50 (2012) 765–781. [36] L.K. Bieniasz, Automatic solution of integral equations pertinent to diffusion with first order homogeneous reactions at cylindrical wire electrodes, J. Electroanal. Chem. 674 (2012) 38–47. [37] L.K. Bieniasz, Automatic simulation of electrochemical transients, assuming finite diffusion space at planar interfaces, by the adaptive Huber method for Volterra integral equations, J. Electroanal. Chem. 684 (2012) 20–31. [38] L.K. Bieniasz, Automatic solution of the Singh and Dutt integral equations for channel or tubular electrodes, by the adaptive Huber method, J. Electroanal. Chem. 693 (2013) 95–104. [39] L.K. Bieniasz, Automatic solution of integral equations describing electrochemical transients under conditions of internal spherical diffusion, J. Electroanal. Chem. 694 (2013) 104–113. [40] L.K. Bieniasz, Automatic solution of integral equations describing electrochemical transients under conditions of internal cylindrical diffusion, J. Electroanal. Chem. 700 (2013) 30–39. [41] L.K. Bieniasz, Automatic solution of integral equations describing electrochemical transients at dropping mercury electrodes, J. Electroanal. Chem. 705 (2013) 44–51. [42] L.K. Bieniasz, A new theory, automatic computation of reversible cyclic voltammograms at a microband electrode, J. Electroanal. Chem. 767 (2016) 123–133. [43] A. Huber, Eine Näherungsmethode zur Auflösung Volterrascher Integralgleichungen, Monatsschr, Math. Phys. 47 (1939) 240–246. [44] W.T. De Vries, Potential-step electrolysis followed by linear-sweep voltammetry, J. Electroanal. Chem. 14 (1967) 75–81. [45] W.T. De Vries, Potential-step electrolysis followed by linear-sweep voltammetry at a plane mercury-film electrode, J. Electroanal. Chem. 16 (1968) 295–312. [46] S.Y. Park, K. Aoki, K. Tokuda, H. Matsuda, Determination of diffusion coefficients of metals dissolved in mercury by double potential step chronoamperometry at hanging mercury drop electrodes, Bull. Chem. Soc. Jpn. 56 (1983) 2133–2137. [47] A. Neudeck, L. Dunsch, Microstructured electrode materials in UV-visible spectroelectrochemistry, J. Electroanal. Chem. 386 (1995) 135–148. [48] IEEE, IEEE 754: Standard for binary floating point arithmetic, http://grouper.ieee. org/groups/754, AAccessed 8 AprilApril, 2017. [49] P.P. Valkó, J. Abate, Comparison of sequence accelerators for the Gaver method of numerical Laplace transform inversion, Comput. Math. Appl. 48 (2004) 629–636. [50] J. Abate, P.P. Valkó, Multi-precision Laplace transform inversion, Intern. J. Numer. Meth. Eng. 60 (2004) 979–993. [51] GWR, http://library.wolfram.com/infocenter/MathSource/4738, Accessed 8 April, 2017. [52] MATHEMATICA, Wolfram Res., Inc., Champaigne, IL, http://www.wolfram.com, Accessed 8 April, 2017. [53] E.N. Blose, B. Ghimire, N. Graham, J. Stratton-Smith, Edge corrections to electromagnetic Casimir energies from general-purpose Mathieu-function routines, Phys. Rev. A 91 (2015) 012501. [54] H. Schmitt, Algorithm 409. Discrete Chebychev curve fit, Commun. ACM 14 (1971) 355–356.

cally, by adjusting the discrete temporal grids to the local variations of the electrochemical responses, to meet accuracy criteria. The present study is restricted to kernel functions and IEs representing pure diffusion. Further work is needed to obtain approximants to integrals ∼ Z κ (t ) associated with kernel functions arising for homogeneous reaction-diffusion systems, and for systems involving other modes of transport, e.g. convection. This should allow a simulation of singular transients arising for such systems. Acknowledgements This research was supported in part by PLGrid Infrastructure. The access to the HP Apollo 8000 Gen9 cluster, and to MATHEMATICA 10, was provided by the Cracow Academic Computing Center CYFRONET AGH. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.jelechem.2017.05.027. References [1] A.J. Bard, L.R. Faulkner, Electrochemical Methods, Fundamentals and Applications, Wiley, New York, 2001. [2] F.G. Cottrell, Der Reststrom bei galvanischer Polarisation, betrachtet als ein Diffusionsproblem, Z. Phys. Chem. 42 (1902) 385–431. [3] C.G. Phillips, K.M. Jansons, The short-time transient of diffusion outside a conducting body, Proc. R. Soc. Lond. A 428 (1990) 431–449. [4] K.B. Oldham, The short-time chronoamperometric behaviour of an electrode of arbitrary shape, J. Electroanal. Chem. 297 (1991) 317–348. [5] J. Koutecký, R. Brdička, Fundamental equation for the electrolytic current when depending on the formation rate of the depolariser jointly with diffusion and its polarographic verification, Collect. Cecoslov. Chem. Commun. 12 (1947) 337–355. [6] P. Delahay, G.L. Stiehl, Theory of catalytic polarographic currents, J. Am. Chem. Soc. 74 (1952) 3500–3505. [7] H.B. Herman, H.N. Blount, Chronocoulometric measurement of chemical reaction rates. The ECE mechanism at plane and spherical electrodes, J. Phys. Chem. 75 (1969) 1406–1413. [8] M. Lovrić, I. Ruži, Extension of an analytical solution for polarographic current influenced by first-order coupled chemical reaction. Planar diffusion model, J. Electroanal. Chem. 146 (1983) 253–261. [9] L. Bieniasz, Influence of diffusion coefficient ratio DO/DR on potential-step chronoamperometric and linear voltammetric current at stationary planar electrodes in the case of a pseudo-first-order EC catalytic reaction scheme, J. Electroanal. Chem. 170 (1984) 77–87. [10] R. Lange, K. Doblhofer, The transient response of electrodes coated with membrane-type polymer films under conditions of diffusion and migration of the redox ions, J. Electroanal. Chem. 237 (1987) 13–26. [11] J.C. Myland, K.B. Oldham, Limiting currents in potentiostatic voltammetry without supporting electrolyte, Electrochem. Commun. 1 (1999) 467–471. [12] S. Cohn, K. Pfabe, J. Redepenning, A similarity solution to a problem in nonlinear ion transport with a nonlocal condition, Math. Models Methods Appl. Sci. 9 (1999) 445–461. [13] L.K. Bieniasz, Analytical formulae for chronoamperometry of a charge neutralisation process under conditions of linear migration and diffusion, Electrochem. Commun. 4 (2002) 917–921. [14] D. Britz, J. Strutwolf, Digital Simulation in Electrochemistry, Springer, Berlin, 2016. [15] J. Strutwolf, W.W. Schoeller, Digital simulation of potential step experiments using the extrapolation method, Electroanalysis 9 (1997) 1403–1408. [16] J. Strutwolf, D. Britz, High order spatial discretisations in electrochemical digital simulation. 2. Combination with the extrapolation algorithm, Comput. Chem. 25 (2001) 205–214. [17] D.J. Gavaghan, An exponentially expanding mesh ideally suited to the fast and efficient simulation of diffusion processes at microdisc electrodes. 2. Application to chronoamperometry, J. Electroanal. Chem. 456 (1998) 13–23. [18] L.K. Bieniasz, Use of dynamically adaptive grid techniques for the solution of electrochemical kinetic equations. Part 6. Testing of the finite-difference patchadaptive strategy on example models with solution difficulties at the electrodes, in one-dimensional space geometry, J. Electroanal. Chem. 481 (2000) 134–151. [19] K. Harriman, D.J. Gavaghan, E. Süli, Adaptive finite element simulation of chronoamperometry at microdisc electrodes, Electrochem. Commun. 5 (2003) 519–529. [20] L.K. Bieniasz, Use of potential-step formulae to reduce computational time in the simulation of linear voltammetry by orthogonal collocation, J. Electroanal. Chem. 208 (1986) 165–168. [21] L.K. Bieniasz, A singularity correction procedure for digital simulation of potentialstep chronoamperometric transients in one-dimensional homogeneous reactiondiffusion systems, Electrochim. Acta 50 (2005) 3253–3261.

51

Journal of Electroanalytical Chemistry 799 (2017) 40–52

L.K. Bieniasz

Cambridge, 2009. [63] J.W. Rayleigh, Note on the numerical calculation of the roots of fluctuating functions, Proc. London Math. Soc. 5 (1874) 112–194. [64] M.K. Kerimov, Overview of some new results concerning the theory and applications of the Rayleigh special function, Comput. Math. Math. Phys. 48 (2008) 1454–1507. [65] J.L. deLyra, On the sums of inverse even powers of zeros of regular BBessel functions, arXXiv:1305.0228v3 [math-ph] 13 Feb 2014, https://arxiv.org/pdf/ 1305.0228.pdf, Accessed 8 April, 2017. [66] R.S. Nicholson, I. Shain, Theory of stationary electrode polarography. Single scan and cyclic methods applied to reversible, irreversible, and kinetic systems, Anal. Chem. 36 (1964) 706–723. [67] J.C. Myland, K.B. Oldham, An analytical expression for the current-voltage relationship during reversible cyclic voltammetry, J. Electroanal. Chem. 153 (1983) 43–54.

[55] Bloodshed Software, Dev-C++, http://www.bloodshed.net/devcpp.html, Accessed 8 April, 2017. [56] Orwell, Dev-C++ Blog, http://orwelldevcpp.blogspot.com, Accessed 8 April, 2017. [57] http://www.netlib.org, Accessed 8 April,2017. [58] W.J. Cody, Algorithm. 715: SPECFUN - A portable FORTRAN package of special function routines and test drivers, ACM Trans. Math. Softw. 19 (1993) 22–32. [59] L.K. Bieniasz, A highly accurate, iinexpensive procedure for computing theoretical chronoamperometric current at cylindrical wire electrodes, Electrochim. Acta 56 (2011) 6982–6988. [60] L.K. Bieniasz, Highly accurate, inexpensive procedures for computing theoretical chronoamperometric currents at single straight electrode edges and at single microband electrodes, J. Electroanal. Chem. 760 (2016) 71–79. [61] A. Pichereau, Sur l’équation tan(x)=x, x réel, http://alain.pichereau.pagespersoorange.fr/equation_tan(x)=x.html,AAccessed 8 April, 2017. [62] P. Flajolet, R. Sedgewick, Analytic Combinatorics, Cambridge University Press,

52