Theory of Polarographic Kinetic Currents for Second-Order Regeneration Reactions at Spherical Electrodes 1. Numerical Solution of Finite-Difference Equations D. T. Pence, J. R. Delmastro, and G . L. Booman Idaho Nuclear Corporation, Idaho Falls, Idaho A numerical finite-difference method of solution, for the equations describing homogeneous reactions with kinetic complications as applied to stationary electrode polarography, is extended to include the effects of spherical diffusion at large values of time. The accuracy and speed of computation were improved by the application of the integral equation for the linear combination of concentrations which obeys the pure diffusion equation, and solution by a simultaneous equation, NewtonRaphson iteration method. The best data from both the finite-difference and steady-state integral equation methods of solution are tabulated. Measurements with spherical electrodes in the near steady-state region are shown to give a one electrode equivalent to a thin layer cell, and to permit the determination of rate constants in a nearly time-independent region by comparison with measurements at a planar electrode.
A GENERALIZED NUMERICAL METHOD for solving stationary electrode polarographic boundary value problems was presented earlier (1). The method involves the simultaneous solution of a series of finite-difference equations which describe the diffusion and kinetic processes of a second-order reaction. The time and space coordinates are transformed into dimensionless parameters which automatically compensate for the time-dependent growth of the diffusion layer. The solution of the problem is presented in the form of a plot of generalized working curves of the ratio of the kinetic current to the diffusion current cs. generalized time. Specific application was made to the large potential-step response for the one-half regeneration mechanism, O + e z R k 2R+O+Z
where 0 is the oxidized form, R the reduced form, Z is an electro-inactive species, and k is the second-order homogeneous rate constant. Since the earlier work was published, the method has been improved in several significant respects. Extension of the time region to include the effects of spherical diffusion at large values of time makes kinetic measurements possible at longer, conveniently accessible times. Unique properties of convergent, spherical diffusion for electrode reactions having homogeneous chemical complications are evident from the numerically derived electrode response curves. These properties of large time asymptotic current ratios of less than the two electron value obtained with planar electrodes and the approach to stationary state behavior are presented and supported by results from steady-state and asymptotic solutions. Efficiency of numerical calculations was nearly doubled through elimination of half the difference equations. This was accomplished by solving for the integral equation representation of the dimensionless concentration parameter which (1) G. L. Booman and D. T. Pence, ANAL.CHEM., 37, 1366 (1965).
arose from the linear combination substitution (2,3). Further time saving was obtained through adopting a simultaneous equation, Newton-Raphson iterative method to efficiently obtain a precise fit to the set of non-linear difference equations. In the approach described in detail in this paper, an initial solution is obtained for each step of time advancement using a predicted value of the concentration variable, and the new concentration vector is then improved by the simultaneous equation, Newton-Raphson method. Because of the very good initial values for the concentration vector, one or two iterations suffice. Because higher signal-to-noise ratios can be obtained by measurement of integrated current and additional information obtained in terms of surface processes, kinetic to diffusion charge ratios are presented in addition to current ratios. A new thin-layer technique is suggested based on the highly convergent diffusion properties of spherical electrodes. The practical choice of electrode size and applicable time region are discussed. Although the specific case of a one-half regeneration reaction mechanism is detailed in this paper, extension to the completely general case of a regeneration mechanism of any order producing any fraction regeneration is straightforward. Also, the corresponding preceding reaction mechanisms involve equations closely analogous to those of the parallel mechanism, and in fact are identical for the first order cases, giving the same response curves when normalized. This duality of response for preceding and parallel reaction mechanisms with potential-step measurements in the limiting current region is developed in detail elsewhere ( 4 ) . ELIMINATION OF LINEAR TERMS
The dimensionless variable $ represents the linear combination of concentrations which obey the pure diffusion equation given as Equation A7 in the Appendix. (3) A complete list of nomenclature is given below as Notation Definitions. The analytical solution of Equation 3 with the required boundary conditions (Equations A3 and A4) can be obtained by the Laplace transform method. Dimensionless equations result with time and distance parameter substitutions similar to those previously described ( I ) , except with T being equal to one half the value defined in the earlier paper. Also, the initial condition for flux at the electrode surface at times near zero-ie., no chemical reaction-is (2) M. Mastragostino, L. Nadjo, and J. M. Saveant, Electrochim. Acta, 13, 721 (1968). (3) M. L. Olmstead, R. G. Hamilton, and R. S. Nicholson, ANAL. CHEM., 41, 260 (1969). (4) J. R. Delmastro and G. L. Booman, 1969. unpublished data. VOL. 41, NO. 6, MAY 1969
737
where T is the square matrix of coefficients for the correction terms 6, given by Equations 10 through 13. where R1 is the R coordinate at the electrode surface. An equation for U1 is obtained independent of the J. variable with the aid of boundary Equations A15 and A17.
+
- ~ 6 a2~ =
~~6~
el - v2,i62 - x2a3= e2
- YnSn - I - V,,& -yM8M1 The integral in Equation 5 cannot be evaluated directly because of the singular point where T equals r. Evaluation can be carried out by assuming U(R1, T ) to be a stepwise defined function. URl, rj+
=
rl(7 > 0) + r2(7> 1
T2)
-
Xn&+
vM,i6M
1
=
=
(10) (11)
8,
8 M
(12) (13)
The coefficients Vn,tinclude the derivative of the kinetic terms K, in Equation A47. The index i is 1-for the first iteration and the ( U,,, + I),, are the values for U obtained by the simultaneous solution of Equation A48. The V,,iterms are defined by
(6)
+ r3
(7
> T 3 ) + . . . . + r, (7 > T,)
With x representing the constant ratio of T,+ l / T j ,and with the initial condition that U(R1, )0, = 0.5, Equation 5 becomes stepwise defined as
where n varies from 2 through M , and k = 2 for n = 2, k = 0 for n = 3 through ( M - l), and k = M for n = M . This application of the Newton-Raphson method involves calculation of the residual vector 8 from Equation 8, the new V,,ffrom Equation 14 to give the new, almost tridiagonal coefficient matrix T, and then the simultaneous solution $ Equation 9 by Gaussian elimination. Corrected values of U are obtained from
(Un,, +di
= (un,j+
d i
-1
+ 8,
(15)
++
where UlV1 is the initial value for U1 at T 5 T,,,,,, U ~ ,isZ the value for U1 at T,,,,, plus the first increment in T,UI,,+ 1 is the value of U 1 at the next serial T ( T j + = xT,). The weighting factors X, are defined by Equation A33. Equation 7 is very suitable for the U 1boundary equation in the set of finite-difference equations. Each time advancement requires the calculation of an additional weighting factor for the sum representing the integral. After many time advancements, the time required to perform the summations represents a significant fraction of the total calculation time. Calculations are still much faster than including the set of simultaneous finite-difference equations for the J, variable, because the time steps are advanced logarithmically (typically 80 steps per decade) and the summation is truncated at large values of the time variable when additional terms become insignificant. NEWTON-RAPHSON ITERATIVE METHOD
The Newton-Raphson Method (5) is applied to the matrix relationship given in the Appendix as Equation A48 in the form e
+
+
w - s = e +
where the vector 8 contcns the residuals resulting from the inexact fit of the vector U to the set of di@ence equations. For the first NewLon-Raphson iteration, U and the kinetic terms K, in the S vector equations are defined using the values of ( U,,,+1 ) 1 obtained from Equation A46. The NewtonRaphson correction terms are obtained from the solution of the matrix equation +
+
m=e
(9)
( 5 ) J. B. Scarborough, “Numerical Mathematical Analysis,” 4th ed., John Hopkins Press, Baltimore, 1958.
738
ANALYTICAL CHEMISTRY
This cycle of calculating 8,V, and c i s repeated until the 6, are sufficiently small. The criteria chosen for completion of the iteration cycles were that the percentage changes in the I ) ~and the (U2,jA1)i were diffusion current ratio, the (UL,,+ less than 0.1%. When these tests were met, the U for the next time point was calculated going back to Equation A48 to obtain the first approximation for the next set of NewtonRaphson iterations with Equation 9. As before ( I ) the ratio of the kinetic to diffusion current is calculated from c
, CHARGE RATIO EQUATIONS
One of the major changes involves the incorporation of kinetic to diffusion charge ratios in addition to current ratios for comparison with experimental data. The reasons for this change are due primarily to instrumentation considerations. In the measurement of charge, the random noise tends to cancel itself out over the integration time period, but in the measurement of current, little analog smoothing is possible. This results from the lower effective bandwidth of the integrator system. Also, for step-functional controlled-potential polarographic measurements, the introduction of a large load resistor and a current follower, which are necessary for good sensitivity, will reduce the frequency response of the potential control circuit (6). With the integrator approach, control circuit bandwidth is retained as the input terminal of the integrator is a wide-bandwidth virtual ground. In addition to obtaining better data, information is also obtained on the initial charge transfer at the electrode, permitting experimental determination of adsorption and related surface charge effects. (6) D. T. Pence and G . L. Booman, ANAL.CHEM., 38, 1112 (1966).
~
Charge ratios were obtained by integration of the electrode current function before selection of the more coarsely spaced values given in the Tables. ASYMPTOTIC LIMITS
Since publishing our first work on this method, additional investigation has been made concerning the upper and lower limits of the generalized working curves. Because the type of reaction being considered involves the change from one to two electron behavior between short and long times at planar electrodes, one would expect the lower and upper kinetic to diffusion current ratio limits to be one and two, respectively. However, for spherical diffusion the upper limit is less than two and is highly dependent on both the degree of sphericity involved and the magnitude of the rate constant. Electrode Response at Very Short Times. In the generalized working curves of kinetic to diffusion current ratio cs. normalized time presented earlier, the asymptote of the current ratio does not appear to approach one as its limit as T approaches zero. This bias, of about one per cent, was noted by Feldberg (7). The approach of the short-time asymptote to unity is dependent on the thickness of the solution segment being considered as equivalent to a semi-infinite layer and can be shown in the following way. At sufficiently short times, CR and Co are completely diffusion controlled and there is negligible kinetic reaction. Under these conditions, Equation A12 reduces to
d2U - + - - =RodU dR2 2 dR Equation 17 can be solved analytically with the following boundary conditions U U
= =
0.5,R 0.0,R
= =
0 10 < T
Figure 1. Kinetic-to-diffusion current ratios for one-half regeneration mechanism Dashed lines indicate stationary-state limiting values. Kineticsphericity factor Q of 0.00 in curve A , 0.06 in B, 0.20 in C, and 1.00 in D size until the probability of escape balances the rate of supply from the electrode surface reaction. The electrode current, during the time this stationary state is being approached, is a measure of the fast kinetic process. This technique would be closely related to the thin film methods developed for planar electrodes, but has the very significant advantage of requiring only one small electrode instead of two electrodes placed very close together, and the material produced at the second electrode is well removed from the small spherical working electrode. Range of Q Values. Very real restrictions exist on the degree of sphericity which can be attained in a practical experiment. The region of sphericity readily available can be seen by expressing the sphericity factor Q, in terms of the ratio of real time to the dimensionless time parameter from Equations A10 and A13.
The spherical geometry constant denoted by a, has values between about 0.03 for a normal hanging drop electrode of 0.05-cm2 surface area, to a value of about 1.0 for an electrode of 50-micron diameter, which is the aperture size of an average dropping mercury electrode capillary. Assuming an observation time of about 10 seconds, the largest Q value would occur in the rising response region. Thus, with a T value of 1, Q values of 0.1 to 3.2 result for the above two electrode sizes. In the rounding-off region, approaching a flat response with a T value of 100, Q values of 0.01 to 0.32 result. In the very flat response region with a T value of 10,000, Q values are in the 0.001 to 0.032 range, giving essentially linear diffusion behavior with electrodes of normal size. Although very small electrodes would be expected to give some difficulty in routine handling, an advantage would accrue in the lower controller bandwidth and current requirements for fast pulse work. Steady-State Solution. A numerical method of solution has been developed and reported ( 9 ) for the equations describing the reaction being considered but under steady-state conditions. This steady-state method provides current and charge ratios which accurately reflect the nature of the electrode response in the range of T of about 100 to infinity. The
1 0 109 0
Figure 2. Kinetic-to-diffusion charge ratios for one-half regeneration mechanism Dashed lines indicate stationary-state limiting value. Kinetic-sphericity factor Q of 0.00 in curve A , 0.06 in B, 0.20 in C, and 1.00 in D numerical solution reported above and the steady-state solution overlap sufficiently to allow the reporting of accurate current and charge ratios throughout the entire range of T. Although the values for the current and charge ratios near T of 100 do not overlap exactly, the agreement is reasonably good. Because the purpose of this paper is to present the most accurate current and charge ratios available, considerable effort was made to reconcile the differences. The asymptotic limits and the values of the current ratios for various values of Q at the dividing point T shown in Table I, provided the basis for reconciliation. Both methods have limitations with respect to the calculation of asymptotic ratio limits. The finite-difference method of solution has an inherent error which is directly related to the approximations introduced by the finite-difference equations. These approximations to a continuous function by finite cell size, cell spacing, and diffusion layer thickness were minimized by devising a “moving boundary” method, accomplished by choosing a diffusion layer thickness of 700 to 7000 drop radii and redividing this solution segment repeatedly during the numerical calculation until the concentration difference in the first two cells was much less than 0.1%. The calculation of the ratio limit at infinite T with the use of the integral equation steady-state method of solution contains an error which is less than 0.25% at Q less than 0.03, but increases with increasing Q. The error is a result of a limitation in an empirically introduced spherical diffusion correction term and is discussed elsewhere (9). In order to take advantage of the merits of both methods of solution, it was necessary to devise the normalizing procedure described above to cover the region of T desired for the indicated values of Q. In general, the data presented in Table I1 and I11 for T of 100 and below were obtained with the use of the finite-difference method of solution, and the data above T of 100 were obtained with the use of the steadystate method of solution. The actual dividing point value of T is given in Table I. Spacing of Diffusion Cells. Although cell spacing in the method reported here is much less of a problem than in methods where the space coordinate is independent of the time parameter, difficulty is still experienced with choice of parameters to adequately cover the concentration profile. A VOL. 41, NO. 6, MAY 1969
743
primary difficulty is the behavior of spherical diffusion in the stationary state which of itself requires 5 X lo5 diffusion cells to meet the requirement of 0.1 % difference in concentration from the electrode surface to the middle of the first cell and a concentration of 0.1 % of the surface concentration at the far boundary. This behavior can be seen from the stationary state equation
Hence a very large solution segment must be considered at long times which is much less linear than the error function relationship existing when linear diffusion is operative. Another difficulty is in matching the diffusion cells to the actual concentration profile of the reacting species. An algorithm giving variable cell spacing would be desirable which could, for example, divide an appropriate segment of solution into 1000 cells, each differing in concentration from the adjacent cells by 0.1 %. Some better “minimum error” criterion may be conceivable. Under the conditions chosen for obtaining the solution to the time-dependent finite-difference problem, namely 800 cells with the far boundary parameter, F, equal to 2.4, the concentration from the electrode surface to the middle of the first cell was found to vary from 0.34% at T = 0.001, to 1.26% and 1.74% at T = 100 with a Q value of 0.01 and 1.00, respectively. Further effort would seem desirable to give a more nearly optimum mesh responsive to changes in the calculated concentration profiles and at the same time to simplify programming through a compiler approach (15). With the time-dependent finite-difference equations, the “movingboundary” method was tried, but failed to significantly improve the results because with each time advancement, the errors due to the poorly defined gradient at the electrode surface are propagated into bad concentration gradients in the region where the new boundary values are taken. Thus, the “moving-boundary” method fails when the concentration gradients are not well defined at a significant distance toward the electrode from the far boundary. For the stationary-state case discussed above, the “moving boundary” technique was very successful, quite probably because of the simplicity of the time-independent solution. NOTATION DEFINITIONS
Constant defined by Equation A26. Spherical geometry constant, D l / z / r o . Defined by Equation A29. Constant defined by Equation A30. Concentration of oxidized form. Bulk concentration of oxidized form. Concentration of reduced form. Stepwise change in value of U at electrode surface. Distance between finite-difference diffusion cells. Change in concentration of U at electrode surface between two successive time advancements. Diffusion coefficient. Constant defined by Equation A31, n = 2,M. Correction terms in Newton-Raphson iterations. n = l,M. Vector of correction terms. ?I
The error function, Erf(x)
nx
;*Joe-ffzdcu.
= -
(15) P. Lindblad and H.Degn, Acta. Chem. Scand., 21, 791 (1967).
744
ANALYTICAL CHEMISTRY
The complement of the error function, 1 - Erf(x). Constant defined by Equation A32. One half the distance to the far boundary for the timedependent case. Defined in Equations A27 and A28. Constant defined by Equation A34. Defined by Equation A35. Increases with each time advancement. Residuals from fitting result of iteration to set of difference equations. n = 1,M. Vector of residuals. Electrode current in presence of kinetic complication. Electrode current controlled completely by diffusion. Index for time advancement. The kinetic contribution in the finite-difference equation. n = 2,M. The homogeneous rate constant for the second-order reaction. Weighting factors for the stepwise representation of the convolution integral occurring in Equation 5 . k = 1,j. Value of distance index, n, at finite-difference cell next to far boundary. Index for distance parameter, R , from electrode surface. n = 1,M. A constant, reflecting the effective electrode sphericity, defined by Equation A13. Dimensionless distance parameter. n = 1,M. Contribution of the equations for the previous time advancement. n = 1,M. Vector of S,. Dimensionless time parameter for each time advancement, starting withj = 1. Integration variable for time coordinate. Value of time for each stepwise definition of U at the electrode surface. Dimensionless concentration variable, giving normalized concentration of the reduced form. U concentration at electrode surface. U concentration at the n’th cell from the electrode surface. U concentration in cell next to far boundary. U concentration in cell chosen at the new boundary in the “moving boundary” method. General notation for U concentration in cell n at timej. Vector of U concentration variable at a single time. Square matrix of coefficients of 6, for solution of Newton-Raphson correction terms. The matrix of coefficients of Un,j+ 1. Coefficients of 6,, calculated for each Newton-Raphson iteration, i. n = 2,M. Vector of coefficients V,,,. Constants defined by Equations A38-40. k = 0,2, and M. Constants defined by Equation A41 for n = 2, and by Equation A42 for n = 3, M - 1. Constants defined by Equation A43 for N = 2, by Equation A44 for n = 3, M - 1 . Time advancement ratio between two adjacent time steps. Dimensionless concentration variable giving the linear combination equation. Concentration of $ at the electrode surface. Electro-inactive species. Sum of weighting factors multiplied by the change in surface concentration of U. Equation A36. APPENDIX
The finite-difference equations for the particular case of the one-half regeneration mechanism (cf. Equations 1 and 2), under large potential-step conditions are derived from the following defining equations in spherical coordinates.
-EUi,j+
I
+ Uz,j+ + S i i
=
+ Y Z U I , ~1+- W Z ~ ? , 1J + XzUa,j+ 1 -YnUn-
-y,?-fuM-
(A31
+--aCR =o ar
ro
To linearize one of the equations and to simplify the form, the following substitutions are made.
J,
=
+
r(1r0
2c0 2 CO" c R )
1,j+
1
- w.+fud4,3+1
Si = GUl,, ro 2C0*
S2 =
These substitutions reduce the defining equations to
Sn
=
0
(A19
+ Sn = 0
To place the equation in dimensionless form, and to permit the space coordinate to automatically compensate for the time dependent growth of the diffusion layer, new time and distance variables, T and R are introduced.
=
0
+H
(A21)
(A22)
+
sLtf=
~ l
- U3,jl - K?
+
C [ u n- 1 , j - 2un,j un+ Dn[un+ l , j - un- I,,] - Kn
Un,,
r
ar2
- 3U2,j
3
- 2r,kCo* U2
+ s.'d
4 + 3C[2u1,, - 3U2,J + u
U2,3
202 - -[4U1,j
L D -az+ at ar2 at
1
=
The notation U,,,means the n'th concentration cell going from n = 2 next to the electrode (n = 1 at the electrode surface) to M for the last concentration cell next to the far boundary, and j is the index of the time point, starting with j = 1 at the chosen initial time at which the finite-difference time advancement is begun. The Sn terms contain the values for U at the previous time point, except for the kinetic contribution Kn. The Sn terms are
u = I(&)
aU = D-a2U -
s 2
('420)
r = ro =
+
(AW
- WoUn,j+ 1 - xnUn+ I,]+
1,j+ 1
The boundary conditions are
r
0
1 3
- 4c[U.w,1 - -UM
u ~ , j
+
1 ~D.M[UM,] -UM3
I,,]
(A23)
+
1 ~ 1
- IJ]
(A24) -
- KM
6425)
In Equations A1 8 through A25, the following definitions apply. ('49)
T
=
kC*t AR
Equation AI and A8 become
a$ = _1 aT
a2$ _
T aRz
R a$ +-2T aR
=
2F/(M - 1)
X = Tj+ i/Tj
(A27)
The constant F is chosen large enough to make the boundary condition Equation A16 hold to the required accuracy. The choice of an F value of 2.4 used to obtain the tabulated data, is discussed above under electrode response at very short times. Because of the way the finite-difference concentration cells are defined, 2F is the value for R at the far boundary
RM+1
=
2F
('428)
Note the removal of the factor of 2 from Equation A10 and from Equation A13 as compared to the notaa factor of l / tion previously used (1). The boundary conditions at the electrode surface in the dimensionless coordinate system become T=0\$=0 R>, O(U = 0
> O)fil R=O
('431)
+ u1 = 1
E
At the far boundary, R.M + 1
=
1
Xk
')I(*)aR
R~ =
('432)
The X k are weighting factors for the stepwise representation of the convolution integral occurring in Equation 5.
The equality of fluxes at the electrode surface is expressed by R > = 0
+ Q(AR)T,!& + (AR)Xi
(g)Ri + ($1
UJQT112
=
(2Xk - x - 1Y2 xk
(A17)
Set of Difference Equations. The elimination of the linear combination variable $ is discussed above. The resultant set of equations in U are expressed in finite-difference equation form as
Q is derived from the representation of the convolution integral at the (j 1)th time point.
+
VOL. 41, NO. 6,MAY 1969
745
j +
fi
=
2
k=
1
Xk 2
(AU)j-
k+ 2
('436)
In Equation A36, AU is the change in value of the surface concentration of U between the ( j 1)th and the jth time intervals. ( A U ) j = U l , j + 1 - U1,j (A37) W, = 1 2A (A39
+
by solution of Equation A48. Gaussian elimination efficiently solves Equation A48 because is nearly tridiagonal. RECEIVED for review October 31, 1968. Accepted February 7, 1969. Work supported by the United States Department of the Interior, Office of Saline Water.
+
Xz
=
2 --(2A 3
+ B2)
('441)
Correct ion Generalized Numerical Method for Stationary Electrode Polarography: Application to Reactions Involving Second-Order Homogeneous Chemical Complications In this article by Glenn L. Booman and Dallas T. Pence [ANAL.CHEM.,37, 1366 (1965)], the following errors appear: On page 1367, the last boundary condition equation after Equation 2 should read r = ro vice r = 0; and the third term on the right hand side of the first equation preceding Equation 15 1
The kinetic contributions are obtained by using a predicted value of U at the ( j 1)th time interval, denoted by (U,,, + and are calculated by linear extrapolation
+
Weighting completely to the future value of U in the kinetic term was found to require less iterations and introduced negligible bias compared with other choices. This completes the equation set required to solve simultaneously for the ( U n , j +1)1 vector. The values of (U,,,+ 1)1 are used as the first trial vector in the Newton-Raphson iteration method. The set of simultaneous equations can be represented in matrix form by +
+
s (A481 where is the square matrix of coefficients of the Un,jk,l terms appearingin Equations A18 through A21, the vector S is defined+by Equations A22 through A25, and the unknown vector U contains the values of Un,,+ 1 which are obtained @U =
746
ANALYTICAL CHEMISTRY
should read __ On page 1368, the right hand side of 2T6R2' Equation 17 should have a AT included as a product in the numerator ; U2 should be added to the right hand side of Equation 21 ; the first equation preceding Equation 22 should read
the first term on the right hand side of the first equation preceding Equation 24 should be positive. The first term on the 2 right hand side of Equation 24 should read - - (-2C 3 Dx)+Av--l. On page 1369, the sign of in the S1expression in the Equation 25 set should be negative. On page 1370, the k ninth line should read kl = k2 = -. On page 1371, Equa2 tions A2, A3, and A4 should all have T as the upper limit of integration; the second and third sentences in the first paragraph following Equation A4 should read, "In Equation 17, is first assumed to be equal to U,. A solution is obtained, and the value of is substituted back into the K term"; the U , term in the third sentence of the following paragraph should be U n , On page 1372, the superscript for D in the the secfirst equation is incomplete and should read D1'2; ond equation on this page should read ZCa vice
+
n2
n,,
u,,
