loves integral equation numerical

Progress In Electromagnetics Research, PIER 97, 357{372, 2009 THE CAPACITANCE OF THE CIRCULAR PARALLEL PLATE CAPACITOR OBTAINED BY SOLVING THE LOVE INTEGRAL EQUATION USING AN ANALYTIC EXPANSION OF THE KERNEL M. Norgren and B. L. G. Jonsson Div. of Electromagnetic Engineering, Royal Institute of Technology Teknikringen 33, SE 10044, Stockholm, Sweden Abstract |The capacitance of the circular parallel plate capacitor is calculated by expanding the solution to the Love integral equation into a Fourier cosine series. Previously, this kind of expansion has been carried out numerically, resulting in accuracy problems at small plate separations. We show that this bottleneck can be alleviated, by calculating all expansion integrals analytically in terms of the Sine and Cosine integrals. Hence, we can, in the approximation of the kernel, use considerably larger matrices, resulting in improved numerical accuracy for the capacitance. In order to improve the accuracy at the smallest separations, we develop a heuristic extrapolation scheme that takes into account the convergence properties of the algorithm. Our results are compared with other numerical results from the literature and with the Kirchho® result. Error estimates are presented, from which we conclude that our results is a substantial improvement compared with earlier numerical results. 1. INTRODUCTION The exact capacitance of the circular parallel plate capacitor, with in¯nitely thin plates, remains an unsolved problem in potential theory, in the sense that to this date no explicit analytical solution has been reported. However, the problem can be formulated as a Fredholm integral equation of the second kind, known as Love's integral equation [1], which can be solved numerically. To our knowledge, the up to date most accurate studies of the capacitance at small plate separations are the ones by Wintle and Kurylowicz [2] and by Carlson and Illman [3]. Both of these studies Corresponding author: M. Norgren ([email protected]). 358 Norgren and Jonsson have been used as benchmarks for solutions obtained by other methods; see e.g., [4, 5, 6]. Wintle and Kurylowicz [2] use an El-Gendi method [7] to rewrite the Love equation and apply numerical integration using the Clenshow-Curtis quadrature method [8] to obtain the capacitance. Carlson and Illman [3] solve the Love equation through an expansion of the kernel into a Fourier-cosine series. Later, they have also extended that method to solve the three-plate problem by means of coupled Love type equations [9]. It is known [1] that for small plate separations, a solution obtained via a series expansion of the kernel converges slowly, requiring a large number of expansion terms. To calculate the expansion coe±cients of the kernel, Carlson and Illman [3] use numerical integration. Hence, their method is limited by a combination of the accuracy of the integration and the large number of terms needed. The accumulated errors e®ectively limit the expansion to about 100 terms, which is insu±cient for convergence at very small separations. In this paper, we show that all of the integrals in the series expansion in [3] can be expressed analytically in terms of the well- studied Sine- and Cosine integrals. In this way, we improve the numerical accuracy of the expansion coe±cients up to the accuracy of the evaluations of the Sine- and Cosine integrals, which makes it possible to increase the number of expansion functions considerably. Hence, in our method, the numerical accuracy in the capacitance is mainly limited by the truncation of the number of expansion functions. Thus, we can present improved results for the capacitance at small plate separations; results that surpass the results in [2, 3], both in correctness and in the signi¯cant numbers of digits. Our results are also in excellent agreement with the result by Kirchho®[10], which becomes increasingly accurate when the plate separation tends to zero [11, 12]. The paper is organized as follows: In Section 2, we review the expansion method, originally presented in [3]. In Section 3, we derive the analytical expressions for all the expansion integrals. Numerical results are presented in Section 4 and Section 5 contains some conclusions. 2. PROBLEM FORMULATION AND INITIAL ANALYSIS The circular parallel plate capacitor is depicted in Figure 1. The distance between the circular plates is denoted dand their common radius is denoted a. The model is idealized in the sense that the plates have zero thicknesses. Following the notation used in many of the previous studies of this problem, we let ·=d=adenote the normalized separation between the plates. Progress In Electromagnetics Research, PIER 97, 2009 359 a d Figure 1. The circular parallel plate capacitor (side view slightly from above). The capacitance of the parallel plate capacitor is [3] C= 4"0aZ1 0f(s) ds; (1) where the function f(s) is the solution to the modi¯ed Love integral equation f(s)¡Z1 0K(s; t)f(t) dt= 1; 0·s·1; (2) with kernel K(s; t) =· ¼·1 ·2+ (s¡t)2+1 ·2+ (s+t)2¸ : (3) Note that in the original derivation by Love [1], the kernel and the function are de¯ned in the range ¡1·s; t·1, and the kernel has only one term, but since f(s) can be shown to be even one can instead use the formulations (2) and (3) (an elegant and short derivation of Love's integral equation can be found in [12]). To solve Equation (2) numerically, we follow the approach in [3] and expand the kernel and the unknown function into the Fourier- cosine expansion functions Ãm(s) =p 2¡±m0cos(m¼s); m = 0;1; : : : ; (4) which in our study have been normalized to ful¯l the orthogonality relation Z1 0Ãm(s)Ãm0(s) ds=±mm0; (5) where ±mm0denotes the Kronecker delta function. Note that similarly tof(s) allÃm-functions have vanishing ¯rst derivatives at s= 0. Carrying out the expansions of f(s) and K(s; t), in terms of fÃmg, 360 Norgren and Jonsson we obtain f(s) =1X m=0fmÃm(s); (6) fm=Z1 0f(s)Ãm(s) ds; (7) K(s; t) =1X m=01X n=0KmnÃm(s)Ãn(t); (8) Kmn=Z1 0Z1 0K(s; t)Ãm(s)Ãn(t) dsdt; (9) which yield the following in¯nite linear system of equations for the coe±cients ffng1 n=0: 1X n=0(±mn¡Kmn)fn=±m0; m = 0;1; : : : : (10) From (1), (6) and the orthogonality of (4), the capacitance reduces to C= 4"0af0; (11) where f0is simply the (0,0)-element in the inverse of the matrix with elements ±mn¡Kmn. In the numerical implementation the matrix is truncated into the size 0 ·m; n·N, where we call Nthe truncation number. It is well-known [11, 3] that at small separations, ·, large values of Nare required to obtain convergence of f0in (11). Also, for the accuracy of the result, it is crucial that the matrix elements Kmn, given by the integrals in (9), have been calculated with a high accuracy. For small values of ·the kernel K(s; t) has a pronounced crest at s=t, making it di±cult to evaluate the integrals in (9) numerically. In [3] this problem has been alleviated, by adding and subtracting a suitable term to the kernel, thereby removing the crest in one of the integrals, and making the inner integral in the other integral available for explicit evaluation. A related procedure has also been used earlier in [2]. However, a remaining problem is that for large values of m; n the expansion functions (4) yield rapidly oscillating integrands, which (when integrated numerically) result in slow convergence and poor accuracy. Hence, it would be bene¯cial if all integrals, encountered when expanding the kernel, could be expressed analytically in terms of established functions, readily available for numerical evaluation with a high accuracy. In the next section we will show that this is indeed the case. Progress In Electromagnetics Research, PIER 97, 2009 361 3. EXPANSION OF THE KERNEL INTO SINE AND COSINE INTEGRALS In this section, we derive the analytical expressions for the expansion of the kernel K(s; t). First, we notice from (3) and (9) that the expansion coe±cients have the property Knm=Kmn, resulting in a symmetric matrix. For the special cases when the indices coincide and/or becomes zero, it is advantageous for the numerical evaluation to derive special simpli¯ed expressions for Kmn. In most of the derivations, we utilize the following properties of the Sine and Cosine integrals [13]: Si(z¤) = [Si( z)]¤; (12) Ci(z¤) = [Ci( z)]¤; (13) Si(¡z) =¡Si(z); (14) Ci(¡z) = Ci( z)¡j¼(0<argz < ¼ ); (15) where ¤denotes complex conjugation and j denotes the imaginary unit. 3.1. Case m=n= 0 By elementary integrals, we obtain K0;0=Z1 0Z1 0K(s; t) dsdt (16) =· ¼Z1 0½Z1 0·1 ·2+ (s¡t)2+1 ·2+ (s+t)2¸ ds¾ dt =1 ¼Z1 0· arctanµt+ 1 ·¶ ¡arctanµt¡1 ·¶¸ dt =1 2¼· 4 arctanµ2 ·¶ ¡·lnµ 1 +4 ·2¶¸ : (17) 3.2. Cases m= 0;n > 0 Here the inner integral is the same as in the previous case, resulting in K0n=p 2Z1 0Z1 0K(s; t) cos( n¼t) dsdt =p 2 ¼Z1 0· arctanµt+ 1 ·¶ ¡arctanµt¡1 ·¶¸ cos(n¼t) dt:(18) De¯ning the function I1(¯; ·; ® ) =Z1 0cos(¯t) arctanµt+® ·¶ dt; (19) 362 Norgren and Jonsson it follows that K0n=p 2 ¼[I1(n¼; ·; 1)¡I1(n¼; ·; ¡1)]: (20) By a change of variable, u=t+® ·, and integration by parts, we obtain I1(¯; ·; ® ) =sin(¯) ¯arctanµ1 +® ·¶ ¡1 ¯Z(1+®)=· ®=·sin¡ ¯(·u¡®)¢ 1 +u2du=sin(¯) ¯arctanµ1 +® ·¶ +1 ¯Imn sin¡ ¯(®+ j·)¢³ Ci¡ ¯(®+ j·)¢ ¡Ci¡ ¯(®+ 1 + j ·)¢´ ¡cos¡ ¯(®+ j·)¢³ Si¡ ¯(®+ j·)¢ ¡Si¡ ¯(®+ 1 + j ·)¢´o :(21) The evaluation of the last integral was carried out in the Maple software, and the result was simpli¯ed using the trigonometric- hyberbolic addition formulas and the properties (12) and (13). Insertion of (21) into (20) and further simpli¯cations yield K0n=p 2 n¼2 Imn cos¡ n¼(1 + j ·)¢ Si¡ n¼(2 + j ·)¢ + cos¡ n¼(1¡j·)¢ Si(¡jn¼·) ¡sin¡ n¼(1 + j ·)¢ Ci¡ n¼(2 + j ·)¢ ¡sin¡ n¼(1¡j·)¢ Ci(¡jn¼·)o :(22) 3.3. Cases m6=n;m > 0;n > 0 Here, we obtain Kmn= 2Z1 0Z1 0K(s; t) cos( n¼t) cos( m¼s) dtds =2 ¼Z1 0Z1 0·· ·2+ (s¡t)2+· ·2+ (s+t)2¸ cos(n¼t) cos( m¼s) dtds: (23) De¯ning the function I2(¯; ·; ® ) =Z1 0·cos(¯t) ·2+ (t+®)2dt; (24) it follows that Kmn=2 ¼Z1 0[I2(n¼; ·; ¡s) +I2(n¼; ·; s )] cos( m¼s) ds: (25) Progress In Electromagnetics Research, PIER 97, 2009 363 Again, using Maple and some simpli¯cations it follows that I2(¯; ·; ® )=Imn sin¡ ¯(®+ j·)¢³ Si¡ ¯(®+ j·)¢ ¡Si¡ ¯(®+ 1 + j ·)¢´ + cos¡ ¯(®+ j·)¢³ Ci¡ ¯(®+ j·)¢ ¡Ci¡ ¯(®+ 1 + j ·)¢´o : (26) Using (12){(15), it follows that I2(¯; ·;¡s) +I2(¯; ·; s ) =¡Imn sin¡ ¯(s+ j·)¢ Si¡ ¯(s+ 1 + j ·)¢ +cos¡ ¯(s+ j·)¢ Ci¡ ¯(s+ 1 + j ·)¢ +sin¡ ¯(s¡j·)¢ Si¡ ¯(s¡1¡j·)¢ + cos¡ ¯(s¡j·)¢ Ci¡ ¯(s¡1¡j·)¢o : (27) Inserting (27) into the expression (25) for Kmn, we can write Kmn=2 ¼I3(n¼; ·; m¼ ); (28) where the function I3(¯; ·; ° ) =Z1 0£ I2(¯; ·;¡s) +I2(¯; ·; s )¤ cos(°s) ds =¡1 2Im½Z1 0·³ sin¡ (¯+°)s+ j·¢ + sin¡ (¯¡°)s+ j·¢´ Si¡ ¯(s+ 1 + j ·)¢ +³ cos¡ (¯+°)s+ j·¢ + cos¡ (¯¡°)s+ j·¢´ Ci¡ ¯(s+ 1 + j ·)¢ +³ sin¡ (¯+°)s¡j·¢ + sin¡ (¯¡°)s¡j·¢´ Si¡ ¯(s¡1¡j·)¢ +³ cos¡ (¯+°)s¡j·¢ + cos¡ (¯¡°)s¡j·¢´ Ci¡ ¯(s¡1¡j·)¢¸ ds¾ =¡1 2Im© I4(¯+°; ¯;j·;1 + j·) +I4(¯+°; ¯;¡j·;¡1¡j·) +I4(¯¡°; ¯;j·;1 + j·) +I4(¯¡°; ¯;¡j·;¡1¡j·)ª ; (29) and where I4(q; ¯; z 1; z2) =Z1 0[sin(qs+¯z1) Si(¯s+¯z2)+ cos( qs+¯z1) Ci(¯s+¯z2)] ds =1 q· sin(¯z1+q) Ci¡ ¯(1 +z2)¢ ¡cos(¯z1+q) Si¡ ¯(1 +z2)¢ ¡sin(¯z1) Ci(¯z2) + cos( ¯z1) Si(¯z2) 364 Norgren and Jonsson + cos( ¯z1¡qz2)³ Si¡ (¯¡q) (1 + z2)¢ ¡Si¡ (¯¡q)z2¢´ ¡sin(¯z1¡qz2)³ Ci¡ (¯¡q) (1 + z2)¢ ¡Ci¡ (¯¡q)z2¢´¸ :(30) To derive (30), we utilized the trigonometric addition formulas and the integration formulas 5.31 and 5.32 in [14]. Summarizing, we obtain Kmn=¡1 ¼ Im© I4¡ (n+m)¼; n¼; j·;1 + j·¢ +I4¡ (n+m)¼; n¼; ¡j·;¡1¡j·¢ +I4¡ (n¡m)¼; n¼; j·;1 + j·¢ +I4¡ (n¡m)¼; n¼; ¡j·;¡1¡j·¢ª :(31) 3.4. Cases m=n >0 In cases when m=n >0 the last two terms in (31) are not suitable for numerical evaluation, since in the ¯nal expression in (30) the limit q!0 must be taken. Alternatively, we can start with setting q= 0 in the integrand in (30); the result reduces to Knn=¡1 ¼Im½ I4(2n¼; n¼; j·;1 + j·) +I4(2n¼; n¼; ¡j·;¡1¡j·) + sin(j n¼·)Z1 0³ Si¡ n¼(s+ 1 + j ·)¢ ¡Si¡ n¼(s¡1¡j·)¢´ ds + cos(j n¼·)Z1 0³ Ci¡ n¼(s+ 1 + j ·)¢ ¡Ci¡ n¼(s¡1¡j·)¢´ ds¾ =¡1 ¼Im½ I4(2n¼; n¼; j·;1 + j·) +I4(2n¼; n¼; ¡j·;¡1¡j·) +j sinh( n¼·)³ (2 + j ·) Si¡ n¼(2 + j ·)¢ ¡j·Si(jn¼·)´ + cosh( n¼·)³ (2 + j ·) Ci¡ n¼(2 + j ·)¢ ¡j·Ci(jn¼·)¡j¼´¾ :(32) 4. NUMERICAL CALCULATIONS The here presented improved expansion is most useful at small separations, wherefore we restrict our study to values ··0:01; results obtained at larger separations are readily available in the literature, with some of the most accurate in [2, 3]. The value of the truncation number is dictated mainly by the computational resources at hand. Here, we have used a maximum value of N= 15000, at the smallest separations. Progress In Electromagnetics Research, PIER 97, 2009 365 All our results will be presented in terms of the normalized capacitance C=C=(4"0a), where 4 "0ais the capacitance between two in¯nitely separated disks. Hence, it follows from (11) that C=f0. 4.1. Extrapolation Schemes To improve our numerical results, we have employed extrapolation. 4.1.1. Power Law Model First, we considered a simple power law model for the capacitance: f0(N) =^C ¡¯(N·)¡®; (33) where f0(N) is the result obtained when using the truncation Nin Equation (10), ^Cis the extrapolated estimate in the limit N! 1 ;® and¯(assuming ® > 0; ¯ > 0) are coe±cients that are determined together with ^C. In our initial tests, the extrapolation parameters, ®; ¯; ^C, were determined by ¯tting (33) to f0(N); f0(N=2) and f0(N=3) (using rounded values for the fractions of Nwhen necessary). In Table 1 we present our results, for various separations ·, for two di®erent values of the product N·. An important observation in Table 1 is that for constant values of N·the values of the parameters ®and¯are approximately constant, regardless of the value of ·. This is especially true for N·= 1, but also for the smaller ·-values when N·= 0:2. The results for N·= 0:2 (farther from convergence) and N·= 1 (closer to convergence) together indicate that the extrapolation scheme over-estimates the capacitance. Another observation from the data in Table 1 is that for each value of N·the amount of extrapolation is approximately the same, regardless of the value of ·. 4.1.2. Heuristic Model for Improvement at Low Accuracies Considering decreasing separations ·, the limitation of the truncation Nwill eventually force us to use smaller values of N·, this delays the convergence. However, the properties of the convergence, as demonstrated in Table 1, o®er an opportunity to improve poorly converged results, by using the following heuristic model: f0(N) =~C+h(N·); (34) where the term ¡¯(N·)¡®in (33) has been replaced by a general function h(N·). Using Equation (34), we have developed a heuristic 366 Norgren and Jonsson extrapolation method that is based on the assumption that when going from one ·-value to the next smaller value the function hremains almost unchanged. To verify this, we have compared the h-functions obtained at three di®erent values of ·. We used Equation (34), setting ~C=^C, with ^Ctaken from the seventh column in Table 1. The results are given in Figure 2, as function of ( N·)¡1for the values between the ones in Table 1. The h-functions in Figure 2 essentially illustrate the convergence, with higher degree of convergence at lower ( N·)¡1 values, where the h-values at ( N·)¡1= 1 indicate the amount left to extrapolate. The three curves are essentially spanning over the same distance, especially those at the smaller separations ·= 0:001 and·= 0:0001, for easier comparison. The nearly overlapping curves indicate that the convergence of the matrix inversion is determined mainly by the product N·. Table 1. The extrapolation parameters, ®; ¯and ^C, obtained when ¯tting the power law extrapolation model (33) to f0(N); f0(N=2) and f0(N=3). The estimate of the true capacitance is C= 4"0a^C. N·= 0:2 · f0(N) ^C ® ¯¢102 0.01 80.235 80.539 1.380 3.30 0.005 158.937 159.240 1.431 3.03 0.002 394.778 395.078 1.472 2.81 0.001 787.647 787.946 1.486 2.73 0.0005 1573.217 1573.516 1.492 2.71 0.0002 3929.640 3929.849 1.496 2.69 0.0001 7856.804 7857.102 1.498 2.69 N·= 1 · f0(N) ^C ® ¯¢103 0.01 80.4312 80.4363 2.358 5.15 0.005 159.1384 159.1436 2.377 5.14 0.002 394.9827 394.9878 2.384 5.16 0.001 787.8533 787.8585 2.386 5.16 0.0005 1573.4238 1573.4290 2.388 5.16 0.0002 3929.8466 3929.8518 2.389 5.16 0.0001 7857.0105 7857.0156 2.391 5.15 Progress In Electromagnetics Research, PIER 97, 2009 367 0 1 2 3 4 5×10 -1κ=0.01 κ=0.001 κ=0.0001 h (N κ) (N κ)-0.2-0.15-0.1-0.050 Figure 2. Comparison of the estimated h-functions, obtained from equation (34), for three di®erent orders of the separation ·. The h- function essentially indicates how close to convergence the calculation ofCis, °atness of the curve for high N·values indicates that Chas essentially the distance jh(Nk)jleft to its approximately true value. Now, having con¯rmed the assumption, the heuristic extrapola- tion method works as follows: iOrder the considered ·-values as ·0> · 1> · 2> : : : , with N; f 0(N);~Candhindexed analogously. iiStarting with ·0, for which we already have assumed a fairly good convergence by using a moderate number of expansion functions (see Figure 2), and using the power law (33), with a truncation number N0, we obtain the extrapolation ^C0. This is taken as the starting value, ~C0=^C0, for our improved algorithm, which from now proceeds repetitively: iiiIncrement the index (denoted i). Use equation (34) to ¯nd the h-function from the previous step: hi¡1(Ni¡1·i¡1) =f0;i¡1(Ni¡1)¡~Ci¡1: (35) Letni¡1be a truncation number ful¯lling Ni·i=ni¡1·i¡1: (36) Assuming that hi=hi¡1, it follows from (34) and (36) that f0;i(Ni)¡~Ci=f0;i¡1(ni¡1)¡~Ci¡1=f0;i¡1µ·i ·i¡1Ni¶ ¡~Ci¡1;(37) from which the extrapolated capacitance becomes ~Ci=f0;i(Ni) +~Ci¡1¡f0;i¡1µ·i ·i¡1Ni¶ : (38) 368 Norgren and Jonsson Note that ni¡1, determined from (36), must ful¯l ni¡1·Ni¡1, and if ni¡1does not becomes an integer the last term in (38) must be evaluated by interpolation between the adjacent integer values. ivRepeat step 3 until the ¯nal (smallest) ·-value has been considered. 4.2. Results for the Capacitance At small separations, the ¯rst approximation to the capacitance is the geometric capacitance Cg=¼=(4·). A much better approximation is the result by Kirchho® [10]: Ck¼¼ 4·+1 4lnµ1 ·¶ +1 4[ln(16 ¼)¡1] +o(1): (39) This formula, which has been proved rigourously by Hutson [11], becomes increasingly accurate as ·decreases. Our results are divided into two cases. In the ¯rst case, we used a constant value N·= 3, and considered separations down to ·= 0:0002; the results are given in Table 2. Here, all extrapolations, ^C, were obtained by ¯tting the power law (33) to the f0-values at N; N= 2; N=3. We also present the relative excess E= (~C ¡ C g)=Cg over the geometric capacitance. Our results are compared with the Kirchho® result (39) and with the numerical results in [2, 3]. In the second case, we considered even smaller separations, down to·= 0:00001; the results are presented in Table 3. Since limited computer memory enforced a maximum N= 15000, all extrapolations Table 2. Numerically calculated capacitances, where f0(N) is the non-extrapolated value and ^Cis the extrapolant. The estimated true capacitance is C= 4"0a^C. In columns 2 and 5, only decimals are presented, since the integers are the same as in column 3. For example, at·= 0:01, where ^C= 80:43451, one has f0(N) = 80 :43440. Our method with N·= 3 References · f0(N) ^C E/x87 (39) [2] [3] 0.01 .43440 80.43451 24.1 .42044 80.4342 80.43 0.005 .14169 159.14179 13.1 .13354 159.13 159.1 0.002 .98596 394.98607 5.82 .98206 394.87 395 0.001 .85661 787.85672 3.13 .85443 787.6 787 0.0005 .42707 1573.42718 1.67 .42588 1573 0.0002 .84994 3929.85005 0.73 .84944 3928.9 Progress In Electromagnetics Research, PIER 97, 2009 369 Table 3. Numerically calculated capacitances, where f0(N) is the non-extrapolated value and ~Cis the extrapolant. In columns 2 and 5, only decimals are presented; same integers as in column 3. Our method with N= 15000 References · f0(N) ~C E/x98 (39) [2] 0.0001 .01294 7857.01378 3.86 .01355 7855.9 0.00005 .16055 15711.16855 2.04 .16847 0.00002 .25402 39273.34241 0.87 .34244 0.00001 .05664 78543.42381 0.46 .42390 10 10 10 1010 10 10 -5 -4 -3 -2-4-3-2 κ|C − C |~ k Figure 3. The magnitude of the di®erence between our extrapolated results and the Kirchho® results from (39), as function of the relative separation ·. Circles denote positive and dots negative di®erences. were obtained with our heuristic method. As starting value, we took the^C-value obtained at the smallest separation considered in Table 2, i.e.·0= 0:0002;~C0= 3929 :85005; see Section 4.1.2. At the largest separation, ·= 0:01, our result agrees better with the reference numerical results than with (39), but for the smaller separations our results are more close to (39). In Figure 3, we have plotted the di®erence between our results and the Kirchho® result, which in e®ect is our numerical approximation of the rest term in Equation (39). For 0 :0002···0:01 the di®erence is positive and decreases with ·with the approximate behavior /·0:8, but when reaching the smallest ·-values the di®erence turns negative and increases in magnitude. This cannot by itself be taken as an error, since (39) is not the exact solution, neither a lower bound, but a solution that becomes increasingly accurate as ·decreases. We should also mention the result by Ignatowsky [15, 2], which di®ers from Kirchho®'s result in that the constant term is replaced 370 Norgren and Jonsson by (ln 8 ¡1=2)=4, thereby falling below the Kirchho® result by the amount of approximately 0 :33447. P¶ olya and SzegÄ o [16] have shown that Ignatowsky's result is a sharp lower bound for the capacitance. Thus, we see in Table 3 that at ·= 0:00001 our non-extrapolated result is below this lower bound, which is an indication that convergence has not been reached due to an insu±cient number of expansion functions. However, the extrapolated result is above the sharp lower bound. 4.3. On the Accuracy of the Results Our numerical simulations have shown that for a ¯xed ·the unextrapolated capacitance f0(N) increases with N, and from Tables 1 and 2 it appears that when using the power law extrapolation formula (33) the extrapolated value ^Cdecreases with N. If we in Table 1 use the ^C-values at N·= 1 as references, the extrapolations obtained at N·= 0:2 overestimate the references with about one third of the total amount of extrapolation. Similarly, if we in Table 2 use the^C-values at N·= 3 as references, the extrapolations obtained at N·= 1 (in Table 1) overestimate the references with about one third of the total amount of extrapolation. Applying this rule to the values in Table 2, we conclude that the true values are approximately 4 ¢10¡5 below the extrapolated value. Using the subsequent repeated extrapolation, the error obtained at·= 0:0002 is propagated to the lower ·-values. If the extrapolation scheme was ideal, i.e., if the h-function in (34) was the same regardless of ·, no further contributions to the error would have occured. To get a very rough estimate of the cumulated errors, we can study the similarities between the curves in Figure 2. Given the excellent agreement between the curves at the smaller ·-values, a gross estimate is that the cumulated error is within 1/10 of the amount of extrapolation. Applying this estimate to the values given in Table 3 we obtain the estimate of the maximum error, denoted ¢ C, acquired in each step of the interpolation scheme. The results, rounded upward to one signi¯cant digit, are presented in Table 4. Since we typically have to let N·decrease with ·, the most signi¯cant contribution to the accumulated error at a certain ·-value is acquired in the last step of the extrapolation scheme. Table 4. Estimated order of the cumulated error during each step of the heuristic extrapolation scheme. · 0.0002 0.0001 0.00005 0.00002 0.00001 ¢C4¢10¡58¢10¡58¢10¡49¢10¡34¢10¡2 Progress In Electromagnetics Research, PIER 97, 2009 371 5. CONCLUSIONS By expanding the kernel in the Love equation into a Fourier cosine series, with the coe±cients expressed analytically in terms of Sine and Cosine integrals, we have increased the accuracy of the expansion, making it possible to use considerably larger truncation numbers N than in previous studies. In this way, we have improved the numerical values for the capacitance at small plate separations ·. The present method enable us to consider smaller plate separation distances than has been considered before, while maintaining a high accuracy. Numerical tests indicated that the degree of convergence is determined by the product N·. Hence, for this method of calculating the capacitance it becomes practically impossible to accommodate for a decreasing plate separation by a corresponding increase of the number of expansion functions; cf. Figure 2. To compensate for this problem, we have developed a heuristic extrapolation scheme that uses the information obtained at intermediate separation. At larger separations, ·, than those considered in our study, convergence is obtained for smaller numbers, N, of expansion functions. In such cases, one does not need to use the analytical results in Section 3 for the integrals, since they can instead be calculated accurately by means of numerical methods [3, 2]. In fact, for large ·- values the explicit expressions in Section 3 are unsuitable for numerical evaluation. This can be seen by observing that from (4) and (9) it follows that for any · >0 the matrix elements ful¯l the relation jKmnj ·2K0;0<2 (40) and that when we tested the algorithm for separations 1 ··· 10 we encountered spurious large matrix elements violating the condition (40). The reason is that when N·À1, the individual terms in the expressions given in Section 3 exhibit an exponential growth /eN·and the condition (40) is met by taking the di®erences between such very large terms, but numerically this kind of evaluation leads to cancelation e®ects. ACKNOWLEDGMENT LJ acknowledges the partial ¯nansiation of The Swedish Agency for Innovation System (VINNOVA) project P33516-1. 372 Norgren and Jonsson REFERENCES 1.Love, E. R., \The electrostatic ¯eld of two equal circular co-axial conducting disks," Quart. J. Mech. and Appl. Math. , Vol. 2, 428{ 451, 1949. 2.Wintle, H. J. and S. Kurylowicz, \Edge corrections for strip and disc capacitors," IEEE Trans. Instr. and Meas. , Vol. 34, No. 1, 41{47, March 1985. 3.Carlson, G. T. and B. L. Illman, \The circular disk parallel plate capacitor," Am. J. Phys. , Vol. 62, No. 12, 1099{1105, December 1994. 4.Nishiyama, H. and Nakamura, M., \Capacitance of disk capacitors," IEEE Trans. Comp., Hybrids., and Manuf. , Vol. 16, No. 3, 360{366, 1993. 5.Donolato, C., \Approximate evaluation of capacitances by means of Green's reciprocal theorem," Am. J. Phys. , Vol. 64, No. 8, 1049{ 1054, August 1996. 6.Hwang, C. O. and J. A. Given, \Last-passage monte carlo algorithm for mutual capacitance," Phys. Rev. E , Vol. 74, 027701, 2006. 7.El-Gendi, S. E., \Chebyshev solution of di®erential equations," Comput. J. , Vol. 12, 282{287, 1969. 8.http:==en.wikipedia.org =wiki=Clenshaw-Curtis quadrature. 9.Carlson, G. T. and B. L. Illman, \Series capacitors and the inverse sum rule," Am. J. Phys. , Vol. 70, No. 11, 1122{1128, November 2002. 10.Kirchho®, G., \Zur theorie des kondensators," Monatsb. Acad. Wiss. Berlin , 731{734, 1877. 11.Hutson, V., \The circular plate condenser at small separations," Proc. Camb. Phil. Soc. , Vol. 59, 211{224, 1963. 12.Sneddon, I. N., Mixed Boundary Value Problems in Potential Theory , Wiley, New York, 1966. 13.Abramowitz, M. and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables , Dover, New York, 1964. 14.Gradshteyn, I. S. and I. M. Ryzhik, Table of Integrals, Series and Products , Academic Press, New York, 1980. 15.Ignatowsky, W., \Kreisscheibenkondensator," Acad. Sci. URSS Trav. Inst. Steklov , Series 2, Vol. 3, 1{104, 1932. 16.P¶ olya, G. and G. SzegÄ o, Isoperimetric Inequalities in Mathemati- cal Physics , Princeton, 1951.