The object of this paper is to revisit a family of improper integrals first investigated by Jaeger (Jaeger 1942 Proc. R. Soc. Edinb., A 61, 223–228). Of particular interest is a sub-class related to axisymmetric diffusion as it occurs in contemporary electrochemistry, specifically chronopotentiometry and chronoamperometry; applications that demand precision well beyond what can be achieved by simple interpolation of the tabulated solutions given by Jaeger and coworkers. To that end, the paper first outlines the less known numerical techniques necessary to solve such integrals and then employs them to obtain numerical solutions over a broad range of the temporal parameter τ, which includes the asymptotic approximations for small and large τ. Computed also are time integrals not previously calculated. More useful to practitioners, however, are approximations to the integrals that are easy to evaluate and sufficiently simple to be manipulated analytically, and the remainder of the work is devoted to such approximants. Each is constructed from a base function which captures the precise asymptotic behaviour at small and large τ plus a correction function, which is fitted either by a half-range Fourier sine series or an inverse power series in τ. Inclusion of sufficiently many terms in each series allows the approximants to realize an accuracy concordant with that of the numerical solutions. The approximants may also be integrated with respect to τ to obtain appropriate time integrals for use in convolution algorithms.
In a series of papers on the flow of heat from circular cylinders, Jaeger and co-workers derive a class of integrals of the form (Jaeger 1942; Jaeger & Clarke 1942; Carslaw & Jaeger 1959) 1.1a 1.1b and 1.1c where 1.1d and Jν and Yν are Bessel functions of order ν of the first and second kind.
Jaeger & Clarke (1942) further derive asymptotic expressions for for various values of the coefficients p and q for small and large values of the temporal parameter τ and comment that values for intermediate τ must be deduced numerically. Proceeding numerically, however, is not straightforward, first because of the singularity introduced in equation (1.1a) by u−1 in the limit , and second because the integrals span the semi-infinite domain. Of course, Jaeger and co-workers knew how to cope with such situations and with a mechanical calculator managed to evaluate for various p, q, although the numerical techniques by which they obtained them are not specified. In testimony to their prowess, their tabulated results remained the benchmark results until the work of Peng et al. (2002) and Britz et al. (2010) decades later.
Meanwhile, integrals in the class 1.1(a)–(d) have arisen in a variety of physical situations dominated by axisymmetric diffusive transport on scales ranging from microchemical (Longinotti & Corti 2007) to geophysical (Becker et al. 2004; Fitzmaurice et al. 2004). These include heat transfer (Carslaw & Jaeger 1959), groundwater flow (Peng et al. 2002; Yeh & Wang 2007) and electrochemistry (Dornfeld & Evans 1969; Aoki et al. 1985a,b,c; Kovach et al. 1985; Szabo et al. 1987; Singleton et al. 1989). In turn, this plethora of applications has demanded solutions to 1.1(a)–(d) to a precision well beyond what can be achieved by simple interpolation of the tabulated solutions given by Jaeger and coworkers. This same demand has led to further attempts to solve 1.1(a)–(d) with p=0 and q=1 numerically (Peng et al. 2002; Yang & Yeh 2007; Britz et al. 2010).
Our intent here is to detail the powerful techniques Jaeger and coworkers likely used to solve 1.1(a)–(d) and our focus is on two integrals of particular interest in electro-analytical chemistry: the first relevant to chronoamperometry (Aoki et al. 1985b, Szabo et al. 1987) and the second relevant to chronopotentiometry (Dornfeld & Evans 1969; Aoki et al. 1985c). That said, the techniques and the computing requirements necessary to exploit such techniques will likely not be at the fingertips of all interested researchers and for that reason our second objective is to provide simple form approximate solutions to 1.1(a)–(d). Of course, approximations too have been attempted before, most notably by Aoki et al. (1985a,b,c), but contemporary needs demand a level of precision beyond what their approximations can achieve. Our intent is to provide approximations of precision comparable with that of our numerical solutions of , which in turn we take care to ensure are accurate.
The integral relevant to chronoamperometry, I(0,1;τ), describes the axisymmetric diffusive transport from a columnar cylindrical conductor in circumstances where a potential (in the conductor) is applied and the resulting temporal current monitored. The current may also be expressed as the convolution of the time derivative of I with the surface concentration of the electroactive component of the fluid in which it is immersed (Mahon & Oldham 2001). Tabulated values of I were first generated by Jaeger & Clarke (1942) and more recently by Peng et al. (2002) and Britz et al. (2010). Aoki et al. (1985b) also evaluated I and proposed an approximate solution to it that has, relative to their numerical solution of I, a maximum error of 1 per cent; Fang et al. (2009) also provided a simple approximation with an accuracy of about 1 per cent. Here, we calculate both I and .
The time integral , on the other hand, is important to chronopotentiometry in which the current is applied and the resulting potential is monitored. In this case, the surface concentration can be obtained from a convolution of the measured Faradaic current with G (Mahon & Oldham 2001). Aoki et al. (1985c) have likewise evaluated and developed an approximate solution to it, albeit accurate to 3 per cent (relative to their numerical solution of ). Convolution algorithms can also be constructed which employ the higher time integral (Oldham 1986) and because of that we calculate in addition to G and . Indeed, convolutive modelling (Mahon & Oldham 1998, 1999, 2000, 2004, Mahon et al. 2002) is a powerful tool for simulating electrochemical experiments. Alternatively, global analysis (Bond et al. 1985, 1994a,b; Mahon 2010) can be used to extract key thermodynamic, kinetic and transport parameters from a single voltammetric experiment where the measured Faradaic current is processed to obtain the surface concentration of electroactive species.
We proceed as follows: we begin (in §2) by revisiting the more difficult of the cases we consider, viz. the case p=0, q≠0 to yield I(0,1;τ), and study in detail the numerical difficulties associated with it and how they are overcome. We then go on to evaluate , followed, in §3 by the p≠0, q=0 cases, G(1,0;τ), and . In §4, we construct simple approximants to I and G that closely recover our numerical results and can be integrated with respect to τ to yield , and . We discuss our results in §5.
2 The numerical solution of the Jaeger integral I(0,1;τ)
We begin with the more difficult of the cases under consideration, viz. the case p=0, q=1 to yield I(0,1;τ), which is defined in equation (2.1). As mentioned above, Aoki et al. (1985b) and others have made attempts to solve equation (2.1), but to date only Peng et al. (2002) and Britz et al. (2010) treat the singular kernel in a rational way. In essence, both split the integration of equation (2.1) into a series of proper integrals that together span the half space and which are each carefully evaluated numerically, the exception being the integral containing the singularity, which is evaluated analytically over a tiny range in the vicinity of the singularity (see §2a(i)). Numerical and procedural differences aside, both studies concur to five decimal places over τ∈[10−2,103], while Britz et al. (2010) provide results to six decimal places over a far greater range of τ, viz. τ∈[10−4,104]. Accordingly, both studies concur with Jaeger & Clarke (1942), typographical errors notwithstanding.
Different here is that we desingularize equation (2.1) by subtracting off the singularity and then adding it back to a complimentary principal-valued integral, which is evaluated analytically over the half-space. Then, as the integrand in the desingularized integral is continuous over the half-space, we can evaluate it directly using Laguerre–Gauss quadrature, which is specific to the semi-infinite domain. In §5, we provide results to 10 significant figures which concur to the appropriate number of decimal places with the previous results of Peng et al. (2002) and Britz et al. (2010).
(a) Subclass 1: p=0, q=1
(i) Case 1: the integral I(0,1;τ)
In evaluating the integral 2.1 we first note that because and as , then the integrand goes as as , which while singular at u=0, is in fact integrable. But unlike Peng et al. (2002) and Britz et al. (2010), who evaluate this integrand over range [0,ϵ], where ϵ≪1, we desingularize the integrand by subtracting and then adding a complimentary principal-valued integral evaluated over the half-space (see Pullin & Phillips 1981, Tjan & Phillips 2007, 2008). Specifically, we write 2.2 Of course, J0(u) and Y0(u) are also zero at multiple points, but because the zeros do not coincide and because 2.3 the integrand in the first term on the right-hand side of equation (2.2) (with α=0) is continuous for for all τ≥0 and thus readily integrable.
More subtle, however, is the behaviour of the integrand of the complimentary integral in the double limit . But the details need not concern us. Rather, because the complimentary integral is identical to equation (2.1) evaluated in the limit , we can employ the asymptotic solution for (2.1) given by Jaeger & Clarke (1942), viz. for sufficiently small τ 2.4 and write in general that 2.5
Our remaining task is to reduce the now well-behaved integral in equation (2.5) to quadrature and to this end, we employ generalized Laguerre–Gauss quadrature which is specific to the semi-infinite interval (Hildebrand 1979, §8.6). This is one of a class of Gaussian quadratures in which the weights wk and abscissa uk are chosen to render the integral exact if the functional f(u) is continuous and polynomial. To wit, we herein employ the weight function uαe−u with α>−1 and write 2.6 where uk is the kth zero of the generalized Laguerre polynomial of degree m and the aforementioned α in equations (2.3) and (2.6) is the order of the polynomial, as (Abramowitz & Stegun 1965, §22.5.17; Hildebrand 1979, p. 394); finally, wk is the weight which can be expressed in a variety of equivalent forms in terms of , for example, 2.7 where Γ is the Gamma function. Finally, the error E at any location ξ along the interval is 2.8
In the context of equation (2.6), therefore, and in view of equation (2.5), we write 2.9 and use double-precision arithmetic to evaluate the sum in equation (2.6) over m terms. Observe that we write wkeuk as . This is done for numerical reasons. In short, because uk can vastly exceed unity for k≫1, we can encounter a situation in which euk exceeds the largest number the computer can handle. In contrast, uk is large as wk is small, so that stays within machine limits. Unknown initially, however, is an appropriate value for m and so we increase it progressively until there is no change in the numerical solution to a specific number of significant figures, a typical value being m=500 for convergence to five significant figures over all considered values of τ. Accordingly, for the results to converge to 10 significant figures over all considered τ, we set m=50 000 with α=0. These were calculated using double-precision arithmetic on a 64-bit Linux PC using a gfortran compiler. Several thousand values of τ are evaluated in less than a second.
Of course, our results must also recover Jaeger’s (1942) asymptotic solution for , viz. 2.10 and do, as we see in figure 1. Here , β=π2/6−γ2 and γ=0.57722, etc., is Euler’s constant. In fact to highlight the asymptotic behaviour, we plot our numerical results as τ1/2I(0,1;τ) along with the asymptotic solutions (2.4) and (2.10) noting from equation (2.4) that as , shown as a solid horizontal line.
(ii) Case 2: , the first integral with respect to τ
Our task now is to evaluate 2.11 using the methods outlined above. Beginning with the asymptotic behaviour, we find on integrating equation (2.4) with respect to τ* from ϵ>0 to τ and taking the limit as , that 2.12 while on integrating equation (2.10) with respect to τ, we deduce as that 2.13 where Ei(y) is the exponential integral and the constant Ω0=3γ/2 (compare §4b(i)).
Our numerical results for along with the asymptotic solutions (2.12) and (2.13) are plotted in figure 2 noting, from equation (2.12), that as . As is clear from the figures, our numerical solution credibly spans between the asymptotic solutions and depicts the appropriate asymptotic behaviour in each limit.
3 The numerical solution of the Jaeger integral G(1,0;τ)
The second class of integrals of interest are those for which p=1 and q=0 to realize G(1,0;τ). Beyond the initial work by Jaeger (1942), little if any work would appear to have been done on G(1,0;τ), although Aoki et al. (1985c) investigated the time integral of G, viz. and provide an approximation for it for τ<4 in terms of powers of τ. Here, we consider in detail the broad range of integrals, viz. G, and .
(a) Subclass 2: p=1, q=0
(i) Case 1: the integral G(1,0;τ)
We first consider G(1,0;τ), viz. 3.1 Here, we note that because J1∼0 and Y1∼u−1 as and that because zeros in J1(u) do not coincide with zeros in Y1(u), the integrand is continuous in allowing us to employ Laguerre–Gauss quadrature directly. Then, in the context of equation (2.6) 3.2
Because equation (3.4) further indicates that τG(1,0;τ)∼π2/8 as , shown as a dashed straight line, we plot our numerical results in figure 3 as τG(1,0;τ) along with the asymptotic solutions (3.3) and (3.4). The numerics closely recover the asymptotics as they should.
(ii) Case 2:
Of interest now is the integral with respect to τ of equation (3.1), viz. 3.5 Here too the integrand is well behaved over the half space and is thus readily integrable, a feature that allowed Aoki et al. (1985c) to proceed using Simpson’s rule, albeit over the finite rather than semi-infinite domain.
From equation (3.3), the asymptotic behaviour as follows as 3.6 and from equation (3.4) as 3.7 Numerical results for are plotted in figure 4 and are again seen to concur closely with the asymptotic solution (3.6).
(iii) Case 3:
Finally, we integrate equation (3.1) a second time with respect to τ of equation (3.1) to find: 3.8 From a numerical viewpoint, we would rather not have the divergent stand-alone τu2 term in the numerator and to obviate that we move the divergence out of the integral by writing v2=τu2 to find: 3.9 In this form, we readily find that the kernel is continuous for all τ>0 and we may proceed as before. Specifically, in the context of equation (2.6), we write, 3.10
Here, the asymptotic behaviour for small τ readily follows from equation (3.3) as 3.11 and from equation (3.7) as 3.12 where the constants are We plot our numerical results for in figure 5 along with the asymptotic solutions (3.11) and (3.12). Again, it is clear from the figures that our numerical solutions credibly span between the asymptotic solutions.
4 Inner and outer approximations
While the evaluation of is relatively straightforward using the techniques outlined in §2, some researchers would instead prefer easy-to-use straightforward formula. Our intent now is to provide just that. Indeed, our mandate is to provide approximations to better than 0.1 per cent of the numerical values of that are readily integrable to , . In order to achieve this, we find it helpful first to split the range of τ into two parts and then introduce approximations for in each, with an appropriate overlap. To that end, we introduce an inner range that spans 0<τ≤2 and an outer range over .
Our simple approximations must, of course, recover the asymptotic behaviour of in each range and for that reason, we tailor our fit around a base function B(τ), which is antecedent to the asymptotic form in each region. We then curve fit the deviation D(τ) of the base function from , as . By default, D is then necessarily zero at τ=0 and as and is typically not zero elsewhere. The question then is by what basis we should fit the deviation, bearing in mind we require something readily integrable, which is preferably rapidly convergent and naturally exhibits the aforementioned zeros. Our first thought is a power series in τ and in fact that worked well in all cases, albeit expressed as circular functions in the inner region and an inverse series in the outer region.
Specifically, a natural choice in the outer region is the form 4.1 which is asymptotic to zero as , the coefficients An being chosen to ensure for all τ≥τ0. Of course, we must specify B prior to evaluating An but our mindset is that B will be a function not dissimilar to the asymptotic form of (see §4b(i) below).
In the inner region, on the other hand, we are dealing with a finite rather than a semi-infinite domain and here it is preferable to craft B to ensure that D has a further zero in the overlap region, as this allows us to employ circular functions as a basis to fit D. To this end, we again craft our base function B from the first few terms of the asymptotic form. We then manipulate the coefficient of one term to ensure that D has a further zero which we place in the overlap region at ca τ=τ0. We then fit D with a Fourier series.
Fourier series are simple to calculate and converge rapidly with an increasing number of terms N, provided the function to be fitted is continuous (Hildebrand 1979, §9.1). Fourier series are also readily integrable and, because the coefficients then go as bn/nω and thus diminish with increasing n (see equation (4.1)), the series retains its convergence properties. Thus, we can approximate and then integrate to obtain an approximation for , , etc. Finally, because design D has zeroes at each end of its domain, we treat it as an odd function in a half-range sine expansion of period T=2τ0 and frequency ω=2π/T truncated to N terms, as 4.2 where the coefficients 4.3
Finally, a curious reader may well ask why we do not also employ Fourier series to fit D in the outer region, albeit with a change of variable to say τ−1 to render the domain finite. The answer, of course, is that we can. But we are then required to integrate terms of the form , which give rise to a class of special functions that belie our quest for simplicity. Alternatively, of course, we could work from say and differentiate with respect to τ, which does leave us with elementary functions. But unfortunately, the coefficients then go as nbn/ω resulting in oscillatory behaviour.
We proceed now to seek approximations for I, namely inner approximations in §4a and outer approximations in §4b.
(a) Inner region approximations 0≤τ≤2
(i) The functions I and
As outlined above, we construct our basis function B(τ) from, say, the first M+1 terms of the asymptotic series (2.4) with the coefficient of the M+1th term set initially to an unknown value, ϱ say. M too is unknown initially but a little experimentation indicated that a good term to manipulate in equation (2.4) is the linear term, so we set M=3 to yield 4.4 Then, since I(0,1;τ) is known accurately from §2a(i), the deviation D(τ) follows as 4.5 Of course, by design D(0)=0, but we should also like D(τ0)=0 say at τ0=2 and solve for the value of ϱ that ensures this through equation (4.5); we find ϱ=0.050555≈1/(2π2). The ensuing deviation is plotted in figure 6. Observe that it has the features we desire: continuous with zeros at each end of the range; a form ideal for fitting with a half-range sine expansion in which T=4. Substitution of D (from equation (4.5)) into equation (4.3) then yields the coefficients for bn, the first 15 of which are listed in table 1.
Our approximation for I(0,1;τ) over 0≤τ≤2 is then 4.6 and on integration with respect to τ from τ=0, we have over 0≤τ≤2, 4.7 To ensure that our approximations depict the correct asymptotic behaviour as , we observe from equation (4.6) that and from equation (4.7) that as they must.
One question, however, remains: how many terms must we include in equations (4.6) and (4.7) to satisfy our tolerance goal; in short, what value must N be to recover the numerical solution to within 0.1 per cent? To answer this, we refer to figure 7a where we plot the absolute percentage error between our simple approximation (4.6) and the numerical solution. First, with N=1, we see that the error is about 1 per cent. This rapid convergence is due largely to our effort in choosing a basis function, which ensured that D have a form not dissimilar to a simple sine function. Setting N=2 lowers the error to around 0.2 per cent and we get close to, and then well surpass, our 0.1 per cent goal with N=3 and N=4, respectively. Finally, the error falls to ca 0.01 per cent with N=10. The findings carry over to equation (4.7).
(ii) The functions G, and
We approximate G by entirely similar means. Here, we choose as our base function 4.8 so, over 0≤τ≤2, 4.9 4.10 and 4.11 Here, we find that D has a zero at τ0=2 when ϱ*=−0.16321 and the coefficients b*n to fit D are as given in table 1. Furthermore, our 0.1 per cent error requirement is realized with N=8 for G, and ; the error for G is shown in figure 7b.
(b) Outer region approximations
Our approach here is largely the same as in §4a, in that we again choose a base function built from the asymptotic behaviour and approximate the deviation of the integral from it with a power series, here in inverse form. Different here, however, is that we do not require the deviation to have two zeros; instead, we require that it relaxes monotonically to zero only as .
(i) The functions I and
For reasons of accuracy, we would like the base function B to be the leading order term in our approximation and the deviation D a smallish correction to it. To that end, our ideal base function B(τ) should mimic I to within a few percent over the range of interest and recover the asymptotic behaviour depicted by equation (2.10) as . Of course, our first trial function for B(τ) is by default equation (2.10), but we have no reason to expect that it will remain within a few percent of I over the full range of interest. Indeed, realizing this criterion can be deceptively difficult and, of all the cases considered, the base function for I turned out to be the most obstinate to construct. The reason in this case is because the asymptotic series (2.10) is an inverse power series not in τ but , which falls below unity once τ<e1+2γ/4, rendering the y−2 and y−3 terms ever more dominant as . We could, of course, drop these terms in constructing B, but convergence towards equation (2.10) is then slow because y is logarithmic in τ. After many variants, we chose to retain equation (2.10) as is and to counter the sensitive behaviour as with a simply integrable additive function (determined from intuition and trial) that vanishes promptly as . That said, we were never able to quell the sensitivity to τ=1 and instead restrict our domain to , where . Our basis function then reads 4.12 with .
Our remaining task is to curve fit D(τ) and, in view of equation (2.10), an appropriate choice would be an inverse power series in y, but alas integrals of inverse y give rise to special functions. Thus, we instead choose a simpler inverse series in τ, albeit one beginning with τ−3/2 to obviate a divergent term on integration. Then, finally, we have 4.13
Unfortunately, the integration of equation (4.12) still introduces the exponential integral, but that is unavoidable. We then find 4.14 where we find numerically that Ω0=0.868300. The coefficients An are determined via singular value decomposition (Press et al. 1986, §2.9, §14.3) and are listed in table 2. Note that in contrast to the Fourier series fits in which we are at liberty to use as many terms in the series as we wish, here we must employ all terms in the series. A plot of the absolute percentage error relative to the numerical solution is given in figure 8, where again we see that the percentage error is well within our upper bound but only for τ<100. For larger τ, the percentage error peaks at around 0.2 near τ=500 and then falls as the approximation relaxes to the asymptotic solution near τ=1000.
(ii) The functions G, and
In contrast to the previous case, our choices for the base function B(τ) and fit to D(τ) for G are straightforward: to wit, we employ the asymptotic behaviour (4.6) as is for our base function and a simple inverse power series to approximate the deviation. As in §4a(i), the coefficients are determined via singular value decomposition and are listed in table 2. Our approximations then read over (): 4.15 4.16 and 4.17 A plot of the absolute percentage error relative to the numerical solution is given in figure 8. Here, we again comfortably satisfy our 0.1 per cent requirement.
The methods employed herein to obtain numerical solutions to specific cases of 1.1(a)–(d) and then approximate them are generic to the broader class of integrals 1.1(a)–(d) and will no doubt be useful for other integrals which arise in the physical sciences. In all examples, our numerical solutions for closely recover appropriate asymptotic solutions and, for the case I(0,1;τ), previous numerical solutions (Peng et al. 2002; Britz et al. 2010) and can thus be considered credible. Accordingly our approximations recover the numerical solutions over the full range of τ with, for the most part, a precision well within our 0.1 per cent goal. Further, as our approximations obviate the need for a broad table of I(p,q;τ) as given by Jaeger & Clarke (1942), we will not provide one, but we do provide small tables of values for use as a check by those wishing to employ our techniques. Specifically, check values for the numerical algorithm to solve directly are given to 10 significant figures in table 3, while those for our approximations are provided in tables 4 and 5.
The most difficult case to approximate was the outer region for I and this in turn is carried through to the error that at τ≈500 peaks at about 0.2 per cent. The reason is partly because we chose a readily integrable simple inverse series, but as the often large coefficients to fit it imply, this is not an optimal choice. Because of this, readers comparing our approximate results for I(0,1;τ) in table 4 with our results found from numerical integration table 3, will observe a slight disparity at τ=102 and τ=103, our approximate values being lower, viz. the approximate being (0.851958, 0.618064) at respective values of τ, compared with our numerical (0.852635, 0.619230) values, which concur with that of Jaeger’s (0.853, 0.620) and of course Peng et al. (2002) and Britz et al. (2010). Nevertheless, our approximations at larger τ precisely recover the asymptotic solution (2.10); for example, at τ=104, we find 0.483014 compared with the asymptotic value of 0.482994. Finally, we address how to apply our approximations to examples specific to convolution analysis with cylindrical electrodes in a forthcoming companion paper (Mahon & Phillips in preparation).
WRCP would like to thank David Lucy and Sergey Suslov for their helpful comments and the Australian Research Council for their support through grants LP088388 and DP1093517.
- Received May 12, 2011.
- Accepted July 5, 2011.
- This journal is © 2011 The Royal Society