- Email: [email protected]
An efficient numerical method of solving the Abel integral equation for cyclic voltammetry
CompurersChm. Vol. 16, No. 4, pp. 311-317, 1992 Printed in Great Britain. AU rights rmervcd
0097-8485/92 s5.00 + 0.00 1992 Pqanl*n press Ltd
tipyright Q
METHOD OF AN EFFICIENT NUMERICAL SOLVING THE ABEL INTEGRAL EQUATION FOR CYCLIC VOLTAMMETRY LESLAW K. ~IJEMIASZ* titute
of Pbysioal Chemistry of the Polish Academy of Soienas, Molten Salts Laboratory, ul. Zagrody 13. 30-318 Cracow, Poland (Received 23 December i991; tr revised form IO AprU 1992)
Abstract-An e5cient, apparently numerically stable and su5cieatly accurate algorithm of solving the Abel integral equation for cyclic voltamtuetry has been proposed. The method is analogous to that of
replacing the kernel of a Fredholm integral equation by a degenerate one. An approximate expression for the kernel (c -z)-m has been suggested which allows the determination of current function values in the range C E [O.OOlSS, 30.1 of the nondimensioaal time parameter C. Depending on the number of integration steps, the algorithm can be faster by l-3 orders of magnitude than the traditional Huber method. The xne.moryrequiremen ts arc also much lower. The method can he extended to other cyclic vohmmetric integral equations involving ({ - z)-ln kernels. A similar approach ~811be wed for an approximate semi-integration, enabling e.g. a fast global analysis of cyolic voltammograms.
INTRODUCTION Cyclic voltammetty (CV) is one of the most powerful and frequently used transient techniques in electroanalytical chemistry (Bard & Faulkner, 1980). Unfortunately, the discussion of experimental CV data is difficult, because for most electrochemical reaction mechanisms an analytical description of voltammograms is not possible and numerical calculations are necessary. One of the techniques used in this context consists of converting kinetic partial differential equations into corresponding integral equations and solving the latter by computer (Macdonald. 1977). This method of theoretica modelllng is less general than a direct solution of the kinetic equations by e.g. a &rite difference approach (Feldberg, 1969; Britz, 1988). For example it is essentialIy not applicable to problems involving homogeneous reactions of order higher than 1, except in a few cases (Nicholson, 1965). However, it usually enables the expression of the influence of certain parameters (such as potential scan rate, bulk concentrations or diffusion coefficients) in a semi-analytical way. This leads to a reduction of numerical calculations and facilitates applications of the theory. Therefore, the method is still employed (see for example Aoki 8c Kato, 1988; Singh et al., 1991; Mirkin & Nilov, 1991; Bieniasx, 1991), after the period of its widest use in the 1960s and 1970s (Nicholson & Olmstead, 1972 and Refs cited therein).
l
Temporary address: Kcmisk Institut, Aarlun Unlversitat, Laugelandsgade 140, DK-8000 Aarhus C. Denmark.
The CV integral equations often involve convolution integrals: r WC - Z)X(Z) d& s0 where the kernel K({ - z) takes the form: K(I - 2) = (I -z)-1’2,
(I)
(2)
C is a non-dimenaiond time parameter and x(c) denotes a non-dimensional CV current function (see e.g. Macdonald, 1977 for standard definitions). As a reprcacntative example the Abel equation: I
:(~-z)~~~~x(~)~={~+B~~II[-~S(~,C.I~~-’(3)
can be considered, with W
I‘) =
c { 2C. - c
for C Q C, for I > G
(4)
where c, is a non-dimensional switching time. This equation corresponds to a reversible redox reaction: Ox+ne-* Red, which takes place at a planar electrode in semi-infinite diffusion conditions. The parameter 13determines the position of a voltammogram along the F axis. A detailed derivation of equation (3) has been presented by Nicholson and Shain (1964). Numerical methods, traditionally used to solve the CV integral equations have heen outlined by Nicholson and Olmstead (1972). They helong to the class of product-integration methods (Baker, 1978) and consist of the following. The considered range of f: is divided into N_ small subintervals OGlGLX of Iength A each (integration steps). Using this grid, integrals (1) are replaced by finite sums resulting
311
312
LE~~AW K. BIENIASZ
from certain quadrature formulas. One looks for x(f) values at points [ =NA, where N=O,l,..., N,. For any N 2 1 the finite sum for integral (1) contains a term involving x(NA) and term(s) involving previously calculated x(c) values. Since the value of the integral (1) is determined by the integral equation, this allows all the x(NA) values to be calculated one after another. However, for every new N the terms involving previously calculated x(C) values need to be recalculated, because the kernel function depends on c_ As a result, when N increases, the computation of x(NA) requires progressively more and more operations and the entire computational time is proportional to N&/2. The integral (1) can also be considered as the so-called semi-integral of the ~(6) function (Oldham, 1981). For the calculation of semi-integrals similar algorithms are employed (see e.g. Bond et al., 1985). The above methods have been implemented in the recently developed user-friendly PC program “ELSIM”, intended for electrochemical kinetic simulations (Bieniasz, 1992). Although they possess good numerical properties (Baker, 1978; Nicholson k Olmstead, 1972), the above inefficiency with respect to the computational time disadvantageously affects the performance of the “ELSIM” program since personal computers are still relatively slow. Therefore, a search for alternative, faster methods is highly desirable. Faster algorithms for semi-integration can also be of use for the so-called global analysis of cyclic voltammograms (Bond ef al., 1985; Zoski et al., 1991). The method of replacing the kernel function by a degenerate one has been used to solve integral equations in various areas of science and technology (see e.g. Krasnov et al., 1971; Kopchenova & Maron, 1981). In this method one looks for a set of M + 1 functions ui(c), wi(z) enabling the following approximate expression for a kernel function K(c, z) of interest: K(c, z, u 5
Ui(C)Witz)*
should be proportional to N_ (and M). Consequently, it should be much lower than for the traditional algorithms described above. Of course, the success of this method depends on finding a suitable approximation (5) for the kernel. No applications of such an approach to the CV integral equations can be found in the electrochemical literature. Therefore, in the present work the possibility of solving the CV equations in this way is examined. The discussion is focused on equation (3), which allows basic properties of the method to be examined and clearly demonstrated, without any influence of additional factors. The uniform grid of [ values, described above, is assumed. PRINCIPLES (I) Elimination kernel
of the problem of the singularity
The main difficulty, which at fist sight appears to preclude any application of the considered method to equation (3), is that the kernel (2) is singular at z = C_Any attempts to find an expansion (5) which would retain this property and enable a direct use of equation (6) seem rather hopeless. However, within the framework of numerical integration this problem can be overcome with little effort, as shown below. Firstly, similar to the assumptions accepted in the methods described by Nicholson and Olmstead (1972), the following two most convenient approximations (7) and (8) for the current function can be considered over each subinterval z E [(k - l)A, kA], wherek=l,2,...,Nand N>l: thestepfunction approximation x(z) = #A),
Although this method is usually discussed in conjunction with Fredholm integral equations (see e.g. Baker, 1978; Krasnov et al., 1971; Kopchenova & Maron, 198 1; Kagiwada & Kalaba, 1974), one can think of a similar approach in the case of Volterra equations (most CV equations belong to that category). In that case one would have: I i K(5*z)X(x) dx = f50 ui(C) W&)X(Z) dz. (6) s0 s0 The advantage of expression (6) is that with a numerical integration procedure the integrals on the righthand side need not be recalculated for every C, because integrated functions do not depend on C, Therefore, the computational time necessary to calculate integral (1) for all C = NA (corresponding to the previously mentioned uniform grid of C values)
(7)
and the linear approximation x(z) y x(kA) + A-‘h(kA) -
(5)
i-0
of the
xKk - 1)4l(z --k4.
(8)
In the following text they are referred to as cases 1 and 2, respectively. Secondly, for every N b 1 the integral (1) in equation (3) can be separated into two components:
oc
(g - z).+x(z)
f
dz =
s(-* s
(5 -
o
zr”*x(z> dz
c
+
(t; - z)-‘“x(z)
dz.
(9)
i-A
The second component can be easily integrated: (I - z)-Wz)
dz =&(C)
+ WC -
=dxWA)+hxW--1)A)l
41 (10)
where g = 2A’j2, h = 0 for case 1 and g = (4/3)A”‘, h = l/2 for case 2.
Solving the Abel integral equation Table I. Approximate valuesof the ceftkkntsa,, 6, of formula (II), as well as cmrespmding collocation points xl, determined by the proceduredcscriw in the Appendix. For i z 11 the absolute accuraq of a, is poor, but this does not influeactthe BECUI(LC~ of fomlula (11)
i
a,
0
0.3010139437B + 00 0.3987372698E +OO 0.61120347288 + 00 0.100738741OE + 01 0.17047444063 + 01 0.2911638772E -I-01 0.4989617665E + 01 0.8561125026B +Ol 0.15405269088 + 02 0.953950656513 tO2 -0.1766463285E +04 0.99629686768 f 06 -0.14434697498 + I2 0.4391864389B + 22 - 0.4777665570E + 43 0.2962918922E +99
1 2 3 4 : 7 : 10 11 12 13 14 I5
C-A Y? = exp(--bJ)
1 0 0 1 1 2 2 3 3 4 4 4 5 s 5 6
3o.ooooO ll.2mol 3.77001 1.27001 0.43101 0.14601 0.04951 0.01681 0.00453 om305 O.aO230 0.00192 0.00172 0.00162 0.00157 0.00155
for the kernel
The approximation for the kernel needs to have the form (S), which significantly limits the number of possibilities. In addition, the most evident polynomial expansion cannot describe simultaneously rapid variations of function (2) at z ‘Y [ and its slow decay when z Q c. Instead, some combination of power functions and exponential functions seems more appropriate. Therefore, in the present work the following particular expression is proposed: (I; -2)-l”
N f
&exP(-b(l
-Z)),
i-0
exp(--b,NA)
Expression (11) can be now substituted into the first integral in the right-hand side of the equation (9). This gives: = f aiyr, i-0
kA
c
k-1
k-l x
[P,x((k
-
114
(12)
+
ax&AX
(14)
where PI=0
(15)
qi = [l - exp(-b,A)]b;’
(16)
for case 1, or p, = -{exp(-_,A) -[l q,={l-[1
- exp(-b,A)](b,A)-‘}b;’
(17)
-exp(-b,A)](b,A)-‘}b;’
(18)
for case 2. From equation (14) it follows that for any N 2 2 the values of y” and yp-’ are interrelated: Y?=
&[Y?-
’ + P,XW
- 2)A) + q&N
l)A)ls (19)
-
where Rt = exp( - b, A).
(20)
By combining equations (3), (9), (10) and (12) one obtains: x(NA) = g-’
[
1 -hx(W -l)A), (21)
F(NA) - 5 a&“’ 1-O
where F(C) represents the right-hand side of equation (3). Finally, after denoting: YN =
a,yYg
-
1
(22)
Pi= a,p,R,g-’
(23)
Q,= w,R,g-’
(24)
equations (19) and (21) can be rewritten: Y;“= Q((N
- 2)A) +
QdO
- l)A) + R, Y;- ’
(25)
for I’ = 0, 1, . . . , M and x(NA) = g-‘F(NA)
rule
(z ) dz
exp(&z)X(z) do. (13) 1 (k-1)A The last expression is valid for N 2 2. When N = 1, one simply has y F = 0. Using the approximations (7) or (8) one obtains: N--l yp = exp(-b,NA) 1 exp(b,kA) =
(11)
which is relatively simple, involves easily calculated and integrable functions and consequently leads to a convenient, recursive algorithm for the function x(c), as shown below. The procedure of finding constants ‘I~,b, is described in the Appendix. Approximate values of these coefficients are listed in Table 1. With M = 15 and coefficients from Table 1, the absolute accuracy of expression (11) is better than 0.005 for 0.00155 < < -z d 30.0. As will be shown, this range and accuracy are fully sticient for solving equation (3).
(C -z)-“2X(z)dz
cxpV&
I0 N-l
0.16666666678 0.1870952637E + 0.75567844608 + 0.24626869168 + 0.7520127382B + 0.22456284028 + 0.66519076683 + 0.1964557424E + 0.5891772I53E + 0.230959962OE + 0.43779635lOE + 0.8573015812E + 0.16448518558+ 0.3242s129753+ 0.64388113lSE+ O.l486313334E.+
The integration range z E 10, [ -A] of the first integral in the right-hand side of equation (9) does not involve the singularity point. Therefore it is convenient to apply an approximate expression (5) for the kernel (2) to this integral only and to modify appropriately the integration limits in equation (6).
(3) Integration
where
b,
(2) The approximation
313
- hx((N - 1)A) -
;
Y”.
(26)
i-0
ALGORITHM (I) Practical formulation A practical algorithm of solving equation (3) results from the above equations. At the beginning
Lesww K. BIENMSZ
314
(N = 0) an initial value of the current substituted: x(O) = 0.
function
(21)
When N = 1 ([ = A) one uses equation
I VI d Ix I,l~rl~p(-b,A)(b,g)-‘.
(10):
x(A) = g -‘F(A).
Simultaneously, all Yj are set to 0, on the basis of equations (13) and (22). Next, for N 2 2 the two-step reeurrenee procedure can be applied to calculate successively all the x(NA) values. In the fust step the integrals Y;” are updated, according to equation (25) and in the second step the corresponding x(NA) value is calculated from equation (26). It is important that even with a minimal A acceptable with the coefficients a, and b, from Table 1 (A = 0.00155), not all the integrals Yy need to be calculated. Equation (13) implies that
where
(29)
lxlmu is an expected maximal value of lx([)l. /*Determine M =
15;
a necessary
number
of
integrals
Initial
Chi
=
condition
0.0;
Chi_2
for
First
for(i Chi
=
chi);
calculated 0;
i c=
n;
= F(Delta)/g;
Calculation8
for(N
*
2;
i++)Y{i] Chi_l
0
l/
for
N ==
'/
1 f/
= 0.0;
= Chi;
Chi);
for
N c=
<
M--;
/+ output
values
printf(1V%10.5g\nu,
ff
N -=
*/
= Chi;
printf("%10.5g\n'*,
/*
~[i]
whila(-b[M]*Delta-log(b[Wl)+log(fabs(aIlq)) log(O.O005*g/Chimax))
/*
/* output
*/
N z-a 2 */
Nmax
; N++)
{ /*
Update
SUm
=
for(i
integrals
Y[i]
l/
0.0; = 0;
i c-
M;
i++)
{Y[i]-P[i]*chi_z+q[i]*~i_l+R[i]*Y[i]; I*
Set
Chi
=
Chi_2
the
solution
F(N*Delta)/g = Chi_l;
printf(*'%10.5g\nn, 1
Chi_l
sum+=y[i];)
*/ - h*Chi_l
- sum;
= Chi; chi);
(30)
A safe estimation is 1~1, < 1. It can also be assumed that every Yr for which the above upper bound on I YfI is smak than e.g. 0.0005 can be neglected [this value is much lower than the maximal error of the ~(5) values resulting from the kernel approximation (ll), as is shown below]. Taking into account the coefficients a, and b, from Table 1, inequality (30) implies that terms with low i dominate in the sum (12). If so, the necessary number M of terms to be summed up equals 12 for A = 0.00155 [cases 1 and 2), 8forA=O.W(cases1and2)and7(ease1)orS(case 2) for A = 0.01. This result is very convenient, because repeated computations involving coeflicents P,, Q, and R, derived from Usand b, values for i > 9 in Table 1 can lead to underflow and/or uncontrolled roundoff errors. This also reduces the computational effort. A relevant fragment of the code in the “C” language (Kernighan & Ritcble, 1978), designed to solve equation (3) can look as follows:
(28)
IvYI d Ix lrrmx exp(--b,A)b;‘,
(29) can be combined with equation (22)
Inequality to give:
is
/* output
*/
Solving the Abel integral equation An extended procedure of this kind has been included into the “ELSIM” program (Bieniasz, 1992) and served for the calculations described in the Results section. The values of a, and b, from Table 1 were used. (2) Error an&sir Three kinds of errors are expected in the considered algorithm: the error resulting from the approximation (1 l), errors due to the replacement of the integrals by !inite sums and machine errors. In order to estimate the effect of the kernel approximation on x(c), equation (3) can be rewritten with the kernel disturbed by a small constant &:
I
,’ KC - z)-‘fi f
&*~~I) dz = e3,
where x*(C) is a corresponding function. After denoting
(31)
disturbed CV current
6x(1) = x(T) -X’(l)
(32)
equation (31) can be combined with equation (3), rearranged and partially solved analyticalIy to give the approximation: I~~(C)Ic I%l(lX lmsx+I&ax)~C”7~-’
(33)
is a maximal value of Ifx(l)l. This where l&u implies that the error is bounded and that Is(c)1 can increase at most proportionally with Cl”. The maximal possible error can be attained at the maximal considered value {, of [, so that one obtains from formula (33): ILxl_=
“J2C~a-’rlXl_.
1 - laJ25$7c-
(34)
one Ix I, = 1, ~1~~~=0.005and I,=30 gets 14Imu= 0.0177, i.e. the maximal error of the current function, induced by the approximation (1 l), is comparable with the error of this expression (.Q. In practice, predictions (33) and (34) are overestimated and la,(C)1 is much smaller. One reason is that the error involved in the approximation (11) is not constant, as assumed in equation (31), but its value and sign vary with C -z. The errors resulting from the replacement of integrals by finite sums can always be reduced by selecting a sufficiently small A value. The integration step A = 0.01 has been checked and found quite satisfactory for the calculation of cyclic voltammograms (Nicholson & Shain, 1964). With the present algorithm any A 2 0.0015 can be used. This can be important for the extensions of the algorithm to more complicated integral equations. In a good algorithm the machine errors should not accumulate too much during the calculations (see e.g. Dahlquist & Bjiirck, 1974; Baker, 1978). A formal analysis of this property of the considered algorithm on the basis of equations (25) and (26) appears rather complex. The resistance of tbe algorithm to errors of this kind was tested experimentally (see Results).
Assuming
315 RESULTS
The operational testing of the proposed algorithm was performed in the following way. First, equation (3) was solved using Huber’s method (Huber, 1939; Nicholson & Ohm&ad, 1972), implemented in the “ELSIM” program (Bieniasz, 1992), with A = 0.004 and iV_ = 7500. This calculation of the CV curve took 39.21 min, when using an IBM compatible PC with an Intel 80386 processor operating at 33 MHz and an 80387 math coprocessor. The value: ln(0) ~6.5, recommended by Nicholson and Shain (1964) was used and {,- - 15 was assumed. The obtained solution was considered accurate. Next, the present algorithm was applied (both cases 1 and 2), using the same integration step. These calculations took about 0.13 min each. Double precision variables were used in all the calculations. The results were compared with the accurate solution. Figure 1 presents the deparatures
of the calculated current function curves from tbe accurate solution. In case 1 the values of ~r”*x(c) did not differ by more than 0.00094, except for a few initial integration steps, where a maximal difference of 0.0106 m. As could be expected, in case 2 departures from the accurate solution were smaller and did not exceed 0.0009. This error was also cu 35 times smaller than the maximal estimate 0.0177~ lP predicted from equation (34). Such deviations of the CV curves from the accurate solution are quite minor, of no importance to any expected practical applications of the calculated data. However, the efficiency gain is considerable: in the particular case of A and ZV_ the
0.5
/T
0.0 % X ‘= *n -0.5
cnse 2
s” E (0
/ casi!
-1.0
forward -1.5
scan I
5
‘everse
1
scan
I
I
I
10
15
20
1
25
??O
5
Fig. 1. Departures of the current function x”~x(C), calcu-
lated using the method proposed in the present work, from the accurate solution by the Huhcr method. A = 0.004.
316
LE~~AWK. BIENIASZ’
present algorithm is about 300 times faster than the classical Huber method. Furthermore, the required memory space is negligible, as compared with the product-integration methods, because neither a vector of iV_ current function values nor a vector of N_ constant coefficients (previously needed to recalculate the convolution integral) have to be stored. The high values of the computational times obtained, and the efficiency gain, resulted to some extent from the slowness of the formula translator used in the “ELSIM” program (Bieniasz, 1992). Therefore, the calculations have also been made without using the formula translator and other facilities included in this program. In this case the simulation times were 4.54 and 0.044 min, for the Huber and present methods, respectively, which gives an efficiency gain around 103. In further tests arbitrary perturbations were added during calculations, either to x(NA), Y” or to all these variables simultaneously in order to check the stability of the algorithm. These computations were performed with A = 0.00155, A = 0.004 and A = 0.01. In one of these tests, random noise, having a uniform distribution over the range (-0.005, +O.OOS), was added to all Yr and to x(NA) at each N. The disturbed voltammogram coincided, on average, with the correct one, while the scatter of its values remained at the same level all the time. Similar results were obtained with other perturbations, which seems to demonstrate that the algorithm is practically stable in the sense that any accidental errors, such as machine errors, do not tend to accumulate during the calculations. CONCLUSIONS The results obtained in this work indicate that the proposed method of replacing the kernel (2) of the integrals (1) with expression (11) leads to a highly efficient and sufficiently accurate numerical algorithm of solving the integral equation (3) for the CV current function. Work is underway to extend this approach to other types of integral equations occurring in the theory of the CV method and involving integrals (1). In the simplest case, such an extension requires equation (26) to be rearranged to obtain F(NA), which is formally equal to integral (1) and to substitute this expression into an integral equation. In the next step, the x(NA) values have to be determined from the resulting equations. Preliminary results confirm the validity of this method. A similar treatment can also be used to enable a fast global analysis of cyclic voltammograms, based on semi-integration of experimental data (Bond et al., 1985). The derived
coefficients aj and bi can be potentially useful for a representation of kernels (2) occurring in integral equations in other areas of science and technology (see e.g. Baker, 1978; Miller, 1971; Krasnov et al., 1971).
REFFRENCE!S
Aoki K. & Kato N. (1988) J. Electroanul. 51.
Chem.
245,
Baker C. T. M. (1978) me Nwnericul Treatment of Integral
Equations. Clarendon Press, Oxford. Bard A. J. & Faulkner L. R. (1980) EIecrrochemical Metho&. Fundamentals and Applications. Wiley, New York. Bieniasz L. K. (1991) J. Electroanul. Chem. 504, 101. Bieniasz L. K. (1992) Computers Chem. 16, 11. Bond A. M., H&derson T: L. E. & Oldham K. B. (1985) J. Electroanof, Chem. 191, 75. Britz D. 119881 Digital Simulrrtion iu Electrochemistrv. Springer; Be& Dahlauist G. & BiBrck A. (1974) Numerical Methods. P&ice-Hall, En&wood Cliffs, NJ. Feldberg S. W. (1969) Eiectroanal. Chem. 3, 199. Huber A. (1939) Mona&t. Math. Phys. 47, 240. Kagiwada H. H. 8r Kalaba R. (1974) Znregral Equutions via Zmbedding Methods. Addison-Wesley, Reading, MA, USA. Kernighau B. W. 8r Ritchie D. M. (1978) The C Programming Language. Prentice-Hall, Inc., Englewood Cllls, NJ. Kouchenova N. V. & Maron I. A. (1981) Commrtalional Mathematics. Worked Examples ‘and, Problems with Elements of Theory, p. 323. Mir Publishers, Moscow. Krasnov M., Kiselev A. & Makarenko G. (1971) ProbIerns and Exercises in Integral Equations. p. 166. Mir Publishers, Moscow. Maedonald D. D. (1977) Transient Techniques in Electrochemistry. Plenum, New York. Miller R. K. (1971) Nonlinear Volterra Integral Equations, p. 78. Benjamin, Menlo Park, CA. Mirkin M. V. & Nilov A. P. (1991), Computers Chem. 15, 55. Nicholson R. S. & Shain I. (1964) Anal. Chem. 34, 706. Nicholson R. S. (1965) Anal. Chem. 37, 667. Nicholson R. S. & Olmstead M. L. (1972) In Computers in Chemistv and Instrumentation (Edited by Mattson J. S., Mark H. B. Jr & MacDonald H. C. Jr), Vol. 2, p. 120. Marcel Dekker, New York. Oldham K. B. (1981) J. Electroanal. Chem. 121, 341. Singh T., Dutt J. & Kaur S. (1991) J. Electroanai. Chem. 304. 17. Zoski*C. G., Oldham K. B., Mahon P. J., Henderson T. L. G. & Bond A. M. (1991) J. Electrounal. Chem. . 297, 1.
APPENDIX In order to determine the coefficients ai and bi for expression
(11) various numerical strategies were attempted, including multiparameter minimization procedures, solutions of simultaneous systems of nonlinear equations resulting from collocation conditions and some other procedures. However, all these methods failed, always leading to divergence, floating point errors or other problems. Finally, the following arbitrary recipe (programmed for a computer) was found to be successful. The aim of this procedure was to find an approximation: A,(x)
= 5 oiexp(-LJ;X) i-0
(AlI
for the functionj(x) = x - “* in the largest possible range of X, such that If(x) - A,&)[ =St, where E is a desired accuracy, M and all a, are as small as possible and Iuil c ((II,,, where Ial,, is an assumed value. The procedure starts with M = 0 and further terms are added successively to the sum (Al), leading gradually to better approximations for f(x), For M = 0 the colloealion
Solving the Abel integral equation
317
condition is applied: f(%)
= a, exP(-4%)
W)
and f’(xo)
= [a0 exp(-W&l
(A3)
to the largest value x0 of x to be considerai. The prime in equation (A3) denotes a derivative with tespect to x. Equations (A2) and (A3) enable calculation of 0, and 4. Next, M is set to 1, the diffbrence D,(x)=&) - A,(x) is taken into consideration and the entire range of x: 0 c x c x, is manned, using a very small increment of x. At every point x new coeSkients a, and b, are tentatively calculated, by requiring: D,(x)
=a,exp(-b,r)
(A4)
and D;(x)
= [u, exp(-b,x)l’.
(AS)
Apointx,,atwhichIf(x)-A,(x)l
0 0,
a0
I
0.5
I
1.0
L
1.5
2.0
X Fig. At. Convergence of the suazessive approximations (Al) with gradually increasing M to the function J(x) = x-“~, The curves were calculated using coefficients u, and b, from Table 1. Respective values of M are indicated in the figure, describe the fkmotion behaviour near x = 0, as illustrated in Fig. Al. The best result was obtained with x,, = 30, 1aI_ = l.E + 100 and L = 0.005 and the corresponding (approximate) values of Al,, b, and xi are listed in Table 1.
0097-8485/92 s5.00 + 0.00 1992 Pqanl*n press Ltd
tipyright Q
METHOD OF AN EFFICIENT NUMERICAL SOLVING THE ABEL INTEGRAL EQUATION FOR CYCLIC VOLTAMMETRY LESLAW K. ~IJEMIASZ* titute
of Pbysioal Chemistry of the Polish Academy of Soienas, Molten Salts Laboratory, ul. Zagrody 13. 30-318 Cracow, Poland (Received 23 December i991; tr revised form IO AprU 1992)
Abstract-An e5cient, apparently numerically stable and su5cieatly accurate algorithm of solving the Abel integral equation for cyclic voltamtuetry has been proposed. The method is analogous to that of
replacing the kernel of a Fredholm integral equation by a degenerate one. An approximate expression for the kernel (c -z)-m has been suggested which allows the determination of current function values in the range C E [O.OOlSS, 30.1 of the nondimensioaal time parameter C. Depending on the number of integration steps, the algorithm can be faster by l-3 orders of magnitude than the traditional Huber method. The xne.moryrequiremen ts arc also much lower. The method can he extended to other cyclic vohmmetric integral equations involving ({ - z)-ln kernels. A similar approach ~811be wed for an approximate semi-integration, enabling e.g. a fast global analysis of cyolic voltammograms.
INTRODUCTION Cyclic voltammetty (CV) is one of the most powerful and frequently used transient techniques in electroanalytical chemistry (Bard & Faulkner, 1980). Unfortunately, the discussion of experimental CV data is difficult, because for most electrochemical reaction mechanisms an analytical description of voltammograms is not possible and numerical calculations are necessary. One of the techniques used in this context consists of converting kinetic partial differential equations into corresponding integral equations and solving the latter by computer (Macdonald. 1977). This method of theoretica modelllng is less general than a direct solution of the kinetic equations by e.g. a &rite difference approach (Feldberg, 1969; Britz, 1988). For example it is essentialIy not applicable to problems involving homogeneous reactions of order higher than 1, except in a few cases (Nicholson, 1965). However, it usually enables the expression of the influence of certain parameters (such as potential scan rate, bulk concentrations or diffusion coefficients) in a semi-analytical way. This leads to a reduction of numerical calculations and facilitates applications of the theory. Therefore, the method is still employed (see for example Aoki 8c Kato, 1988; Singh et al., 1991; Mirkin & Nilov, 1991; Bieniasx, 1991), after the period of its widest use in the 1960s and 1970s (Nicholson & Olmstead, 1972 and Refs cited therein).
l
Temporary address: Kcmisk Institut, Aarlun Unlversitat, Laugelandsgade 140, DK-8000 Aarhus C. Denmark.
The CV integral equations often involve convolution integrals: r WC - Z)X(Z) d& s0 where the kernel K({ - z) takes the form: K(I - 2) = (I -z)-1’2,
(I)
(2)
C is a non-dimenaiond time parameter and x(c) denotes a non-dimensional CV current function (see e.g. Macdonald, 1977 for standard definitions). As a reprcacntative example the Abel equation: I
:(~-z)~~~~x(~)~={~+B~~II[-~S(~,C.I~~-’(3)
can be considered, with W
I‘) =
c { 2C. - c
for C Q C, for I > G
(4)
where c, is a non-dimensional switching time. This equation corresponds to a reversible redox reaction: Ox+ne-* Red, which takes place at a planar electrode in semi-infinite diffusion conditions. The parameter 13determines the position of a voltammogram along the F axis. A detailed derivation of equation (3) has been presented by Nicholson and Shain (1964). Numerical methods, traditionally used to solve the CV integral equations have heen outlined by Nicholson and Olmstead (1972). They helong to the class of product-integration methods (Baker, 1978) and consist of the following. The considered range of f: is divided into N_ small subintervals OGlGLX of Iength A each (integration steps). Using this grid, integrals (1) are replaced by finite sums resulting
311
312
LE~~AW K. BIENIASZ
from certain quadrature formulas. One looks for x(f) values at points [ =NA, where N=O,l,..., N,. For any N 2 1 the finite sum for integral (1) contains a term involving x(NA) and term(s) involving previously calculated x(c) values. Since the value of the integral (1) is determined by the integral equation, this allows all the x(NA) values to be calculated one after another. However, for every new N the terms involving previously calculated x(C) values need to be recalculated, because the kernel function depends on c_ As a result, when N increases, the computation of x(NA) requires progressively more and more operations and the entire computational time is proportional to N&/2. The integral (1) can also be considered as the so-called semi-integral of the ~(6) function (Oldham, 1981). For the calculation of semi-integrals similar algorithms are employed (see e.g. Bond et al., 1985). The above methods have been implemented in the recently developed user-friendly PC program “ELSIM”, intended for electrochemical kinetic simulations (Bieniasz, 1992). Although they possess good numerical properties (Baker, 1978; Nicholson k Olmstead, 1972), the above inefficiency with respect to the computational time disadvantageously affects the performance of the “ELSIM” program since personal computers are still relatively slow. Therefore, a search for alternative, faster methods is highly desirable. Faster algorithms for semi-integration can also be of use for the so-called global analysis of cyclic voltammograms (Bond ef al., 1985; Zoski et al., 1991). The method of replacing the kernel function by a degenerate one has been used to solve integral equations in various areas of science and technology (see e.g. Krasnov et al., 1971; Kopchenova & Maron, 1981). In this method one looks for a set of M + 1 functions ui(c), wi(z) enabling the following approximate expression for a kernel function K(c, z) of interest: K(c, z, u 5
Ui(C)Witz)*
should be proportional to N_ (and M). Consequently, it should be much lower than for the traditional algorithms described above. Of course, the success of this method depends on finding a suitable approximation (5) for the kernel. No applications of such an approach to the CV integral equations can be found in the electrochemical literature. Therefore, in the present work the possibility of solving the CV equations in this way is examined. The discussion is focused on equation (3), which allows basic properties of the method to be examined and clearly demonstrated, without any influence of additional factors. The uniform grid of [ values, described above, is assumed. PRINCIPLES (I) Elimination kernel
of the problem of the singularity
The main difficulty, which at fist sight appears to preclude any application of the considered method to equation (3), is that the kernel (2) is singular at z = C_Any attempts to find an expansion (5) which would retain this property and enable a direct use of equation (6) seem rather hopeless. However, within the framework of numerical integration this problem can be overcome with little effort, as shown below. Firstly, similar to the assumptions accepted in the methods described by Nicholson and Olmstead (1972), the following two most convenient approximations (7) and (8) for the current function can be considered over each subinterval z E [(k - l)A, kA], wherek=l,2,...,Nand N>l: thestepfunction approximation x(z) = #A),
Although this method is usually discussed in conjunction with Fredholm integral equations (see e.g. Baker, 1978; Krasnov et al., 1971; Kopchenova & Maron, 198 1; Kagiwada & Kalaba, 1974), one can think of a similar approach in the case of Volterra equations (most CV equations belong to that category). In that case one would have: I i K(5*z)X(x) dx = f50 ui(C) W&)X(Z) dz. (6) s0 s0 The advantage of expression (6) is that with a numerical integration procedure the integrals on the righthand side need not be recalculated for every C, because integrated functions do not depend on C, Therefore, the computational time necessary to calculate integral (1) for all C = NA (corresponding to the previously mentioned uniform grid of C values)
(7)
and the linear approximation x(z) y x(kA) + A-‘h(kA) -
(5)
i-0
of the
xKk - 1)4l(z --k4.
(8)
In the following text they are referred to as cases 1 and 2, respectively. Secondly, for every N b 1 the integral (1) in equation (3) can be separated into two components:
oc
(g - z).+x(z)
f
dz =
s(-* s
(5 -
o
zr”*x(z> dz
c
+
(t; - z)-‘“x(z)
dz.
(9)
i-A
The second component can be easily integrated: (I - z)-Wz)
dz =&(C)
+ WC -
=dxWA)+hxW--1)A)l
41 (10)
where g = 2A’j2, h = 0 for case 1 and g = (4/3)A”‘, h = l/2 for case 2.
Solving the Abel integral equation Table I. Approximate valuesof the ceftkkntsa,, 6, of formula (II), as well as cmrespmding collocation points xl, determined by the proceduredcscriw in the Appendix. For i z 11 the absolute accuraq of a, is poor, but this does not influeactthe BECUI(LC~ of fomlula (11)
i
a,
0
0.3010139437B + 00 0.3987372698E +OO 0.61120347288 + 00 0.100738741OE + 01 0.17047444063 + 01 0.2911638772E -I-01 0.4989617665E + 01 0.8561125026B +Ol 0.15405269088 + 02 0.953950656513 tO2 -0.1766463285E +04 0.99629686768 f 06 -0.14434697498 + I2 0.4391864389B + 22 - 0.4777665570E + 43 0.2962918922E +99
1 2 3 4 : 7 : 10 11 12 13 14 I5
C-A Y? = exp(--bJ)
1 0 0 1 1 2 2 3 3 4 4 4 5 s 5 6
3o.ooooO ll.2mol 3.77001 1.27001 0.43101 0.14601 0.04951 0.01681 0.00453 om305 O.aO230 0.00192 0.00172 0.00162 0.00157 0.00155
for the kernel
The approximation for the kernel needs to have the form (S), which significantly limits the number of possibilities. In addition, the most evident polynomial expansion cannot describe simultaneously rapid variations of function (2) at z ‘Y [ and its slow decay when z Q c. Instead, some combination of power functions and exponential functions seems more appropriate. Therefore, in the present work the following particular expression is proposed: (I; -2)-l”
N f
&exP(-b(l
-Z)),
i-0
exp(--b,NA)
Expression (11) can be now substituted into the first integral in the right-hand side of the equation (9). This gives: = f aiyr, i-0
kA
c
k-1
k-l x
[P,x((k
-
114
(12)
+
ax&AX
(14)
where PI=0
(15)
qi = [l - exp(-b,A)]b;’
(16)
for case 1, or p, = -{exp(-_,A) -[l q,={l-[1
- exp(-b,A)](b,A)-‘}b;’
(17)
-exp(-b,A)](b,A)-‘}b;’
(18)
for case 2. From equation (14) it follows that for any N 2 2 the values of y” and yp-’ are interrelated: Y?=
&[Y?-
’ + P,XW
- 2)A) + q&N
l)A)ls (19)
-
where Rt = exp( - b, A).
(20)
By combining equations (3), (9), (10) and (12) one obtains: x(NA) = g-’
[
1 -hx(W -l)A), (21)
F(NA) - 5 a&“’ 1-O
where F(C) represents the right-hand side of equation (3). Finally, after denoting: YN =
a,yYg
-
1
(22)
Pi= a,p,R,g-’
(23)
Q,= w,R,g-’
(24)
equations (19) and (21) can be rewritten: Y;“= Q((N
- 2)A) +
QdO
- l)A) + R, Y;- ’
(25)
for I’ = 0, 1, . . . , M and x(NA) = g-‘F(NA)
rule
(z ) dz
exp(&z)X(z) do. (13) 1 (k-1)A The last expression is valid for N 2 2. When N = 1, one simply has y F = 0. Using the approximations (7) or (8) one obtains: N--l yp = exp(-b,NA) 1 exp(b,kA) =
(11)
which is relatively simple, involves easily calculated and integrable functions and consequently leads to a convenient, recursive algorithm for the function x(c), as shown below. The procedure of finding constants ‘I~,b, is described in the Appendix. Approximate values of these coefficients are listed in Table 1. With M = 15 and coefficients from Table 1, the absolute accuracy of expression (11) is better than 0.005 for 0.00155 < < -z d 30.0. As will be shown, this range and accuracy are fully sticient for solving equation (3).
(C -z)-“2X(z)dz
cxpV&
I0 N-l
0.16666666678 0.1870952637E + 0.75567844608 + 0.24626869168 + 0.7520127382B + 0.22456284028 + 0.66519076683 + 0.1964557424E + 0.5891772I53E + 0.230959962OE + 0.43779635lOE + 0.8573015812E + 0.16448518558+ 0.3242s129753+ 0.64388113lSE+ O.l486313334E.+
The integration range z E 10, [ -A] of the first integral in the right-hand side of equation (9) does not involve the singularity point. Therefore it is convenient to apply an approximate expression (5) for the kernel (2) to this integral only and to modify appropriately the integration limits in equation (6).
(3) Integration
where
b,
(2) The approximation
313
- hx((N - 1)A) -
;
Y”.
(26)
i-0
ALGORITHM (I) Practical formulation A practical algorithm of solving equation (3) results from the above equations. At the beginning
Lesww K. BIENMSZ
314
(N = 0) an initial value of the current substituted: x(O) = 0.
function
(21)
When N = 1 ([ = A) one uses equation
I VI d Ix I,l~rl~p(-b,A)(b,g)-‘.
(10):
x(A) = g -‘F(A).
Simultaneously, all Yj are set to 0, on the basis of equations (13) and (22). Next, for N 2 2 the two-step reeurrenee procedure can be applied to calculate successively all the x(NA) values. In the fust step the integrals Y;” are updated, according to equation (25) and in the second step the corresponding x(NA) value is calculated from equation (26). It is important that even with a minimal A acceptable with the coefficients a, and b, from Table 1 (A = 0.00155), not all the integrals Yy need to be calculated. Equation (13) implies that
where
(29)
lxlmu is an expected maximal value of lx([)l. /*Determine M =
15;
a necessary
number
of
integrals
Initial
Chi
=
condition
0.0;
Chi_2
for
First
for(i Chi
=
chi);
calculated 0;
i c=
n;
= F(Delta)/g;
Calculation8
for(N
*
2;
i++)Y{i] Chi_l
0
l/
for
N ==
'/
1 f/
= 0.0;
= Chi;
Chi);
for
N c=
<
M--;
/+ output
values
printf(1V%10.5g\nu,
ff
N -=
*/
= Chi;
printf("%10.5g\n'*,
/*
~[i]
whila(-b[M]*Delta-log(b[Wl)+log(fabs(aIlq)) log(O.O005*g/Chimax))
/*
/* output
*/
N z-a 2 */
Nmax
; N++)
{ /*
Update
SUm
=
for(i
integrals
Y[i]
l/
0.0; = 0;
i c-
M;
i++)
{Y[i]-P[i]*chi_z+q[i]*~i_l+R[i]*Y[i]; I*
Set
Chi
=
Chi_2
the
solution
F(N*Delta)/g = Chi_l;
printf(*'%10.5g\nn, 1
Chi_l
sum+=y[i];)
*/ - h*Chi_l
- sum;
= Chi; chi);
(30)
A safe estimation is 1~1, < 1. It can also be assumed that every Yr for which the above upper bound on I YfI is smak than e.g. 0.0005 can be neglected [this value is much lower than the maximal error of the ~(5) values resulting from the kernel approximation (ll), as is shown below]. Taking into account the coefficients a, and b, from Table 1, inequality (30) implies that terms with low i dominate in the sum (12). If so, the necessary number M of terms to be summed up equals 12 for A = 0.00155 [cases 1 and 2), 8forA=O.W(cases1and2)and7(ease1)orS(case 2) for A = 0.01. This result is very convenient, because repeated computations involving coeflicents P,, Q, and R, derived from Usand b, values for i > 9 in Table 1 can lead to underflow and/or uncontrolled roundoff errors. This also reduces the computational effort. A relevant fragment of the code in the “C” language (Kernighan & Ritcble, 1978), designed to solve equation (3) can look as follows:
(28)
IvYI d Ix lrrmx exp(--b,A)b;‘,
(29) can be combined with equation (22)
Inequality to give:
is
/* output
*/
Solving the Abel integral equation An extended procedure of this kind has been included into the “ELSIM” program (Bieniasz, 1992) and served for the calculations described in the Results section. The values of a, and b, from Table 1 were used. (2) Error an&sir Three kinds of errors are expected in the considered algorithm: the error resulting from the approximation (1 l), errors due to the replacement of the integrals by !inite sums and machine errors. In order to estimate the effect of the kernel approximation on x(c), equation (3) can be rewritten with the kernel disturbed by a small constant &:
I
,’ KC - z)-‘fi f
&*~~I) dz = e3,
where x*(C) is a corresponding function. After denoting
(31)
disturbed CV current
6x(1) = x(T) -X’(l)
(32)
equation (31) can be combined with equation (3), rearranged and partially solved analyticalIy to give the approximation: I~~(C)Ic I%l(lX lmsx+I&ax)~C”7~-’
(33)
is a maximal value of Ifx(l)l. This where l&u implies that the error is bounded and that Is(c)1 can increase at most proportionally with Cl”. The maximal possible error can be attained at the maximal considered value {, of [, so that one obtains from formula (33): ILxl_=
“J2C~a-’rlXl_.
1 - laJ25$7c-
(34)
one Ix I, = 1, ~1~~~=0.005and I,=30 gets 14Imu= 0.0177, i.e. the maximal error of the current function, induced by the approximation (1 l), is comparable with the error of this expression (.Q. In practice, predictions (33) and (34) are overestimated and la,(C)1 is much smaller. One reason is that the error involved in the approximation (11) is not constant, as assumed in equation (31), but its value and sign vary with C -z. The errors resulting from the replacement of integrals by finite sums can always be reduced by selecting a sufficiently small A value. The integration step A = 0.01 has been checked and found quite satisfactory for the calculation of cyclic voltammograms (Nicholson & Shain, 1964). With the present algorithm any A 2 0.0015 can be used. This can be important for the extensions of the algorithm to more complicated integral equations. In a good algorithm the machine errors should not accumulate too much during the calculations (see e.g. Dahlquist & Bjiirck, 1974; Baker, 1978). A formal analysis of this property of the considered algorithm on the basis of equations (25) and (26) appears rather complex. The resistance of tbe algorithm to errors of this kind was tested experimentally (see Results).
Assuming
315 RESULTS
The operational testing of the proposed algorithm was performed in the following way. First, equation (3) was solved using Huber’s method (Huber, 1939; Nicholson & Ohm&ad, 1972), implemented in the “ELSIM” program (Bieniasz, 1992), with A = 0.004 and iV_ = 7500. This calculation of the CV curve took 39.21 min, when using an IBM compatible PC with an Intel 80386 processor operating at 33 MHz and an 80387 math coprocessor. The value: ln(0) ~6.5, recommended by Nicholson and Shain (1964) was used and {,- - 15 was assumed. The obtained solution was considered accurate. Next, the present algorithm was applied (both cases 1 and 2), using the same integration step. These calculations took about 0.13 min each. Double precision variables were used in all the calculations. The results were compared with the accurate solution. Figure 1 presents the deparatures
of the calculated current function curves from tbe accurate solution. In case 1 the values of ~r”*x(c) did not differ by more than 0.00094, except for a few initial integration steps, where a maximal difference of 0.0106 m. As could be expected, in case 2 departures from the accurate solution were smaller and did not exceed 0.0009. This error was also cu 35 times smaller than the maximal estimate 0.0177~ lP predicted from equation (34). Such deviations of the CV curves from the accurate solution are quite minor, of no importance to any expected practical applications of the calculated data. However, the efficiency gain is considerable: in the particular case of A and ZV_ the
0.5
/T
0.0 % X ‘= *n -0.5
cnse 2
s” E (0
/ casi!
-1.0
forward -1.5
scan I
5
‘everse
1
scan
I
I
I
10
15
20
1
25
??O
5
Fig. 1. Departures of the current function x”~x(C), calcu-
lated using the method proposed in the present work, from the accurate solution by the Huhcr method. A = 0.004.
316
LE~~AWK. BIENIASZ’
present algorithm is about 300 times faster than the classical Huber method. Furthermore, the required memory space is negligible, as compared with the product-integration methods, because neither a vector of iV_ current function values nor a vector of N_ constant coefficients (previously needed to recalculate the convolution integral) have to be stored. The high values of the computational times obtained, and the efficiency gain, resulted to some extent from the slowness of the formula translator used in the “ELSIM” program (Bieniasz, 1992). Therefore, the calculations have also been made without using the formula translator and other facilities included in this program. In this case the simulation times were 4.54 and 0.044 min, for the Huber and present methods, respectively, which gives an efficiency gain around 103. In further tests arbitrary perturbations were added during calculations, either to x(NA), Y” or to all these variables simultaneously in order to check the stability of the algorithm. These computations were performed with A = 0.00155, A = 0.004 and A = 0.01. In one of these tests, random noise, having a uniform distribution over the range (-0.005, +O.OOS), was added to all Yr and to x(NA) at each N. The disturbed voltammogram coincided, on average, with the correct one, while the scatter of its values remained at the same level all the time. Similar results were obtained with other perturbations, which seems to demonstrate that the algorithm is practically stable in the sense that any accidental errors, such as machine errors, do not tend to accumulate during the calculations. CONCLUSIONS The results obtained in this work indicate that the proposed method of replacing the kernel (2) of the integrals (1) with expression (11) leads to a highly efficient and sufficiently accurate numerical algorithm of solving the integral equation (3) for the CV current function. Work is underway to extend this approach to other types of integral equations occurring in the theory of the CV method and involving integrals (1). In the simplest case, such an extension requires equation (26) to be rearranged to obtain F(NA), which is formally equal to integral (1) and to substitute this expression into an integral equation. In the next step, the x(NA) values have to be determined from the resulting equations. Preliminary results confirm the validity of this method. A similar treatment can also be used to enable a fast global analysis of cyclic voltammograms, based on semi-integration of experimental data (Bond et al., 1985). The derived
coefficients aj and bi can be potentially useful for a representation of kernels (2) occurring in integral equations in other areas of science and technology (see e.g. Baker, 1978; Miller, 1971; Krasnov et al., 1971).
REFFRENCE!S
Aoki K. & Kato N. (1988) J. Electroanul. 51.
Chem.
245,
Baker C. T. M. (1978) me Nwnericul Treatment of Integral
Equations. Clarendon Press, Oxford. Bard A. J. & Faulkner L. R. (1980) EIecrrochemical Metho&. Fundamentals and Applications. Wiley, New York. Bieniasz L. K. (1991) J. Electroanul. Chem. 504, 101. Bieniasz L. K. (1992) Computers Chem. 16, 11. Bond A. M., H&derson T: L. E. & Oldham K. B. (1985) J. Electroanof, Chem. 191, 75. Britz D. 119881 Digital Simulrrtion iu Electrochemistrv. Springer; Be& Dahlauist G. & BiBrck A. (1974) Numerical Methods. P&ice-Hall, En&wood Cliffs, NJ. Feldberg S. W. (1969) Eiectroanal. Chem. 3, 199. Huber A. (1939) Mona&t. Math. Phys. 47, 240. Kagiwada H. H. 8r Kalaba R. (1974) Znregral Equutions via Zmbedding Methods. Addison-Wesley, Reading, MA, USA. Kernighau B. W. 8r Ritchie D. M. (1978) The C Programming Language. Prentice-Hall, Inc., Englewood Cllls, NJ. Kouchenova N. V. & Maron I. A. (1981) Commrtalional Mathematics. Worked Examples ‘and, Problems with Elements of Theory, p. 323. Mir Publishers, Moscow. Krasnov M., Kiselev A. & Makarenko G. (1971) ProbIerns and Exercises in Integral Equations. p. 166. Mir Publishers, Moscow. Maedonald D. D. (1977) Transient Techniques in Electrochemistry. Plenum, New York. Miller R. K. (1971) Nonlinear Volterra Integral Equations, p. 78. Benjamin, Menlo Park, CA. Mirkin M. V. & Nilov A. P. (1991), Computers Chem. 15, 55. Nicholson R. S. & Shain I. (1964) Anal. Chem. 34, 706. Nicholson R. S. (1965) Anal. Chem. 37, 667. Nicholson R. S. & Olmstead M. L. (1972) In Computers in Chemistv and Instrumentation (Edited by Mattson J. S., Mark H. B. Jr & MacDonald H. C. Jr), Vol. 2, p. 120. Marcel Dekker, New York. Oldham K. B. (1981) J. Electroanal. Chem. 121, 341. Singh T., Dutt J. & Kaur S. (1991) J. Electroanai. Chem. 304. 17. Zoski*C. G., Oldham K. B., Mahon P. J., Henderson T. L. G. & Bond A. M. (1991) J. Electrounal. Chem. . 297, 1.
APPENDIX In order to determine the coefficients ai and bi for expression
(11) various numerical strategies were attempted, including multiparameter minimization procedures, solutions of simultaneous systems of nonlinear equations resulting from collocation conditions and some other procedures. However, all these methods failed, always leading to divergence, floating point errors or other problems. Finally, the following arbitrary recipe (programmed for a computer) was found to be successful. The aim of this procedure was to find an approximation: A,(x)
= 5 oiexp(-LJ;X) i-0
(AlI
for the functionj(x) = x - “* in the largest possible range of X, such that If(x) - A,&)[ =St, where E is a desired accuracy, M and all a, are as small as possible and Iuil c ((II,,, where Ial,, is an assumed value. The procedure starts with M = 0 and further terms are added successively to the sum (Al), leading gradually to better approximations for f(x), For M = 0 the colloealion
Solving the Abel integral equation
317
condition is applied: f(%)
= a, exP(-4%)
W)
and f’(xo)
= [a0 exp(-W&l
(A3)
to the largest value x0 of x to be considerai. The prime in equation (A3) denotes a derivative with tespect to x. Equations (A2) and (A3) enable calculation of 0, and 4. Next, M is set to 1, the diffbrence D,(x)=&) - A,(x) is taken into consideration and the entire range of x: 0 c x c x, is manned, using a very small increment of x. At every point x new coeSkients a, and b, are tentatively calculated, by requiring: D,(x)
=a,exp(-b,r)
(A4)
and D;(x)
= [u, exp(-b,x)l’.
(AS)
Apointx,,atwhichIf(x)-A,(x)l
0 0,
a0
I
0.5
I
1.0
L
1.5
2.0
X Fig. At. Convergence of the suazessive approximations (Al) with gradually increasing M to the function J(x) = x-“~, The curves were calculated using coefficients u, and b, from Table 1. Respective values of M are indicated in the figure, describe the fkmotion behaviour near x = 0, as illustrated in Fig. Al. The best result was obtained with x,, = 30, 1aI_ = l.E + 100 and L = 0.005 and the corresponding (approximate) values of Al,, b, and xi are listed in Table 1.