The accurate evaluation of a particular Fermi-Dirac integral

ELSEVIER

Computer Physics Communications I01 (1997) 47-53

The accurate evaluation of a particular Fermi-Dirac integral N. Mohankumar, A. Natarajan Safely Research and Health Physics Group. lndira Gandhi Centre for Atomic Research. Kalpakkam, Tamilnadu, India 603102 Received 12 July 1996; revised 22 November 1996

Abstract In this paper, we provide a simple and very accurate numerical method to evaluate a particular two-dingnsional Fermi-Dirac integral that arises in the plasma transport theory. The proposed method involves suitable coccdina~ and variable ,.ransformations which eliminate variable limits of integration. The resulting inner imegral is evaluated by a

trapezoidal integration with pole correction. The outer integral is evaluated by the IMT scheme. We provide reasonable error bounds. The resui:s agree very well with the results of Fullerton and Ringer, who employ a combination of quadrmm:, asymptotic methods and Chebyshev fitting. PACS: 02.70. + d: 52.25.Fi; 97.10.Cv Keywords: Degenerate el~tron gases; Stellar structure; Numerical q~_L=~__ratute;IMT scheme; Trapezoidal scheme with pole cm'rectioo

1. Introduction Considerations of statistical mechanics are needed to describe the transport and thermal properties of degenerate electron gases such as those occurring in stellar atmospheres or in condensed matter like solids. The solutions to such problems are often conveniently expressed in terms of the standard Fermi-Dirac integrals or the generalised Fermi-Dirac integrals [ 1,2]. The latter are defined as

Iv( a, O) = fo t" ( 1 _.___e,_. + Ot/2)l + t/2

dr.

(])

The parameters a and 0 are related to the degree of degeneracy and the temperature respectively of the medium. ~, takes integer and half integer values. A 0010-4655/~/$17.00 Copynght© t ~ PII S0010-4655(96)00166-X

P.b~shed by E~ev~ ~

similar but two-dimensional integral arises in wobiems dealing with partially degenerate plasmas [3,4]. This integral is defined as ~:

H.-fo dt

tv

:g

51t

l + e r _ ~ f t ds l + e ~ _ . ,

(2)

where /~ takes integer values. The set of pL and J, values of common interest are indicated in Table !. The solution procedure for the statistical mechanics problems referred above is usually iterative, requiring repeated evaluation of the Fermi-Dirac integrals for different values of the parameters. Since these integrals cannot be expressed in closed form, one resorts to functional flUs or to direct numerical evaluation. Recently, the authors have reported an efficient numerical algorithm in the case of half order (i.e., half integer u values) Fermi-Dirac inte8.V.A. ngh¢ n m r ~ .

N. Molumk~mar. A. Natarajan / Computer Physics Communications !01 (:0o7) "7 5a

48

2. I. The IMT scheme

Tab4e I Set o f ~ and r va/ues o f c o m m o n interest 9z

v

o.o

0.o.o.5. 1.5. 2.5 -0.5.0..5. i.5 - 0.5.0.5 -05

I0

20 3.0

g~!s [5,6]. The algorithm is based on the following two interesting observations about the trapezoidal rule. (i) If the integrand is analytic but for some poles in ar infinite strip around the real axis ano ,,ties off like e x p ( - x - ' ) as x tends to ( + x ) , then the trapezoidal rule over the full infinite range can be very. accurate even with large step sizes. (ii) Though the presence of poles within the strip of analyticity, and close to the line of integration, can severely affect the accuracy, it can be compensated exactly by including certain correction terms involving the residues of the poles. The reported algorithm thus offers, p drops, the best attainable combination of accuracy and computational economy. The aim of this paper is to extend these ideas for the evaluation of the two-diniensional Fermi-Dirac integral in Eq.

(2). The proposed numerical scheme involve~, alr,o the application of an additional quadrature scheme, known as the IMT scheme. This again is a trapezoidal scheme but after a specific variable transformation. It is more robust than the widely used Gaussinn methods in handling certain singularities [7]. In order to apply these qurdratt,re rules for the integral in F_,q. (2), the integral has to be cast into an appropilate form by means of a suitable transformation of variables. The paper identifies these transformations. The accuracy of the trapezoidal quadrature schemes with any kind of variable transformation needs to be well supported by suitable error analysis. This point is discussed in the appendix to the paper.

This numerical quadrature scheme was introduced by Iri, Moriguti and Takasawa [7]. This is essentially a trapezoidal scheme after a particular variable transformation. The integral, the variable transformation, and the approximation to the integral are defined below #=

(~)

-]
fou0 ' ( t )

x--~b(u)=

d',

q=

/o Ie -

I/tll

-t)

dr,

(4) dx I du = qS'(u) = --e q

-II't'-";,

~O(u) = f [ x = ok(t,)] ~'(t,),

(5)

i N-i

#

:So 4,(,,)

,, = ~

--/

J N

d. =

,, ,.z, +<,,,>.

r =,6{u ~

~- x l J .

,v-i

1

!= E w,f( x,),

(6)

wi= -~b'(ui).

I-1

Here N is the number of quadrature steps. The transformation derivative 4/(u) is such that the transformed integrand ¢,(u) and all its derivatives vanish at both the ends of the interval, namely at 0 and i. The classical Euler-Maclaurin formula expresses the discretisation error, in terms of the odd order derivatives at both the ends of the interval. Thus the derivatives vanishing idcrt:ica!!y facilitates the rapid asymptotic convergence of the trapezoidal sum to the integral. Suppose that at the 0 and 1 ends the integrand behaves like x " and (1 - x ) ~, respectively, where a and 13 are both greater than - 1 . Then, following Iri et al. [7], the discretisation error b.a hlr this method can be given as

E~ ~ , , , ( N , ,~) + b,( iv, t~), 2. An outline of the IMT scheme and the trapezoidal scheme with pole correction

with (tl+

~ ( u , ,,) = Before we describe our method, we briefly outline the IMT scheme and t ~ trapezoidal scheme with pole correction which are used in the outer and inner integrations, respectively.

(7)

I) I/4+a e - ¢41r(~1" I)N

(eq)"+t(2~rN) 3/4+" + ,)u

+

+

(s)

N. Mohankumar A. Natarafi, n / Cnml,uter Physics Cnmrnu:Jivati¢,rt~ tOl (1997) 47-53

Here a and b are some constants and the expression for ¢(N, ~ ) is exactly same as that of ¢ ( a . N). except that a is replaced by ft.

2.2. The trapezoidal scheme with pole correction We folio,.,, Bialeckl [8] to give an estimate of the error of the trapezoidal quadrature with pole correction. More details on this method and its relation to the classical Euler-Maclaurin fommla and the Abel-Plana formula cart be found in McNamee [9]. Let the required integral over ( - x , •) be

l= L=/( x) dx.

(9)

Define C to be the contour which encloses an infinite rectangular region of width 2d. symmetric about the x axis. This contour has two vertical sections, one at - ~ c and the other at + ~c. Assume that f(Z) is analytic but for some poles within C. The region enclosed by C may be called as the strip of analyticity for the function f(z). Define g ( z ) as below,

=

I

flz)

(2rri)

( z -- x) s i n ( , r z / h )

.

(lO)

Application of the Residue theorem to the integral of g ( z ) around C results in f( x)

-

,

sin[(

~ f( ktt)

lh)( X - - kh)] x

. . . .

-

kh)]

49

Here "h' is the step-size :nd E~ is the discretisation error which is the difference between the actual integral and the sum of trapezoidal and pole correction terms. If the nature of f ( = ) is assumed to be such that in Eq. (! 2) the contributions of the vertical ~ctions of C vanish, then E d can be written as [5.6] E,, -" e - -'''a/h

M,

(14)

where M is given by M=

(e -'~''/~' - - e - "~.,/h)

f(r-i a) +

)

(e~2,,,/,,_e_;,,t/h)

dr.

(15)

This expression shows that as long as the expression M. involving the integral, grows algebraically, the exponential multiplier e -2=a/h will dominate over M and attenuate the error Ea. Note that the parameter of importance is the ratio (d/h). For instance, a value of 6 for this ratio is usually adequate for convergence upto 14 digits. While a large d permits large h. it should be borne in mind that the number of pole correction terms go up thereby. therefore an optimum balance has to be struck. For an appropriately chosen d, we fix h as d/6, ensuring the discretisation error to be below the machine limits. Armed with these preliminaries, we proceed to evaluate the required integral as follows.

+ sin(~rx/h) ~ R e , ; g( : ) ] .

! sin(Trx/h)

3. The transformation and evauation of H,~

a~

f ( z ) dz

. ( z - x) sin(~'z/n)

2rri

(11)

Here the residues are evaluated at z~. which are the poles of f ( z ) . The above expression can be integrated with respect to x leading to

t=rcos(O),

s=rsin(O),

r"~-$2+l 2. (16)

With these transformations, tile integral becomes

-tp

~

The first step is the transformation from the Cartesian (t. s) plane to the polar (r, 0) plane dofined as below.

,=-~

~.a,

- E f ~ Res

_=

--3g

g(z)l:, sin(.,rx/h)dx

(12)

f~ri,~

t.+l

x dr ft.12

cos" ( 0 ) sin *"( 0 ) d 0

. ,ft.* [i +

t L~sin(Trx/k) d.ff~,

Ea = 7.'a'i

=

f(-)d( r,, - x ) s i n ( , r z / h )

(13)

(iT) This transformation ha,~ the adv :ntage of changing the variable limit of the inner integral into a fixed

50

N Mohank~umar, A. Natarajan / Computer Physics Communications lO! f ! 997) 47-53

Table 2 Vai-aes o f the integral H,~, for sample values o f v, tt and a . The entries in fl'~ last column indicate lhe number o f digits upto which Ihe values o f the present method match with the values reported in [3] u

~

o~

Integral

- 0.5 0.5 -- 0.5 - "~.5 - 0.5 - 0.5 - 0.5 - 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

!.0 1.0 i.0 IO ! .0 ! .0 I0 ! .9 1.0 !.0 1.0 1.0 i .0 1.0 1.0 1.0

- 200.0 - | ~.0 0.0 I00,0 200.9 300.0 400.0 500.0 - 200.0 - 100.0 0.0 100.0 200.0 300.0 400.0 500.0 200.0 - 100.0 0.0 100.0 200.0 300.0 400.0 500.0 - 200.0 - 100.0 0.0 t00.0 200.0 300.0 400.0 500.0 - 200.0 - ! 00.0 9.0 I f~ }.0 2G/3.0 3( 0.0 4O0.0 500.0 - 200.0 - 100.0 0.0 100.0 200.0 300.0 400.0 500.0

0.30003864136393 - 173 0.216807 ! 3519199D - 86 0.76771541408831D + 00 0.800327 !6 2 2 4 2 7 6 D + 05 0.45259473747301D + ',~ 0.1247133459 ! 292D + 07 0.25600657069448D + 07 0A4722094378602D + 07 0.10501352447684 - 173 0 . 7 5 8 8 2 4 9 7 3 1 6 8 1 0 D - 87 0.34088960831M98D + 00 0.19058526330546D + 07 0.21553014249678D + 08 0.89082586548580D + 08 0.24381828482601D + 09 0.5324093~472960D + 09 0.10126304!47080 - 173 0 . 7 3 1 7 2 4 0 8 1 3 4 8 6 7 D - 87 0.38311560058603D + 00 0.88955293319367D + 08 0.20116998700088D + 10 0.12471794626953D + ! i 0.45513221446220D + ! i 0.12422968366913D + 12 0.64508307892987 - 173 0 . 4 6 6 t 3 5 3 4 0 6 6 0 9 0 D - 86 0.18162026039910D + 0! 0.57208354971578D + 07 0.64(,,68330727302D + 08 0.26726483334480D + 09 0.731481 [ 4926623D + 09 0.1597264t),t)90346D + I 0 0.26628429~,~2505 - ! 73 0.19241633249 ~;41D - 86 0.941114208730S9D + 90 0.14836572571854D + u'~ 0.33534455920793D + l0 0.20788019349731D + I I 0.75858855923919D + 1 ! 0.20705557183250D + 12 0.1991506a820640 - 172 0.14390573598447D - 85 0.59 i 49187902653D + 01 0.44542750221086D + 09 0.100622014(X)501D + I I 0.62369192931543D + 11 0.227587 ! 0524399D + 12 0.62118512016168D + 12

1.5

1.5 1.5 !.5 1.5 1.5 !.5 ! .5 -0 5 -0.5 -0.5 0.5 -0.5 - ( .5 0.5 - 0.5 0.5 0.5 0.5 0.5 0.5 0.~ 0.5 0.5 -0.5 -0.5 0.5 0.5 -0.5 0.5 0.5 -0.5 -

i .0

1.0 !.0 1.0 1.0 1.0 1.0 ! .0 2.0 2.0 2.0 2.G 2.0 2.0 2.0 2.0 2,0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0 3.0

10.99 10.99 10.98 10.34 !0.3 10.47 11.31 10.76 9.7 I 9.71 ! 1.01 10.71 10.44 9.90 9.74 9.72 9.69 9.69 10.61 9.64 9.94 9.43 9.42 9.5 I ! 1.85 11.86 I 1.48 10.25 I0.01 10.02 10.22 I 0.53 i 0.20 I 0.20 11.95 10.06 10.70 9.88 9.82 9.88 I I. 17 ! I. 17 I ! .09 10.15 10.22 10,43 11.42 10.52

N. Mohankumar, A. Natarajan / Computer Physics Communicatmns !0I t 1997) 47-53

one. We next make the substitutiot: r = u 2 and then interchange the order of integrations to get

H`'~,= f~/2 dO cos"(0) sin~'(0) ~./4 u '-(~+~')+ 3

du

×/~_ [1 +¢=:~o:,0,-o11, +e.~.,o,o,-.] " (is) With g and t, being integer and half integer, respectively, the inner integrand is an even function of u which allows us to use the trapezoidal scheme with pole correction as described earlier. The integrand has simple poles at the set of points z~k and Zzk, given by

i

z,k = +

a :1: i(2k + I)¢r cos(0)

]1/2

'

(19)

a + i ( 2 k + 1)1; ,/2

Z2k = :t:

(20)

sin(0)

Here, k takes nonnegative mteger values. For a given value of 'J, the inner integral can be approximated as u '('+*')÷ 3 du

]\ l,

o][,

e~(o)

51

=

., 2~v+ ~-*-!) AZt *-2t -- E, 2 ~ i n ( O ) [ l + C|a+=12t+S)~'i~ot(0)- " ] " (24)

The infinite trapezoidal sum in Eq. (22) is truncated to N I terms. Denoting this truncated sum by SNI(O), the outer intcgr~iion appro~,imatcd by an N point IMT quadrature sum is as follows: N

H`'~, = E wi[cos( 0j)] `"[sin( 0j)]~' [ S,vOj)( /-1

+ R,( Oi) + R2( 01)] .

(25)

The above sum gives the approximation to the required integral. Table 2 gives the sample values of the integral over a range of parameter values. There is a minimum 9 digit agreement between our results and the sample results of Fullerton and Rinker [3] who claim about I0 digit accuracy for the Chebyshev fit they employ. In the appendix we prove that the results can be easily trusted up to 10 digits.

ol

-:s~(o) + e, + e,.,

(21)

4. C v n c l u s i o n

where

:c

ld~ v + ~)+ 3

u,=ih.

(22)

and R1(0) and R,(O) correspond to the correction terms for the poles at z=, and z2,, respectively, and are given by ci~z|tStt/h

A,t = s;.n[.n.Ztk/h] , St, = S g n [ I m ( z i , ) ] , R,(O) = ~"lk -

E

k

2 ~o~(O)[! + =('+'~"+"']='~'-°]'

(23) ei~ :z=$z=Ih

~2k --"

sin[ lrz2t/h ] "

S2t = S g n [ l m ( z , , ) l ,

We have indicated a numerical method which can be used to evaluate a particular two-dimensional Fermi-Dirac integral used in plasma transport theory, with economy and accuracy. The success of the present method is partly based on the exploitation of the Euler-Maclaturin formula. The economy of the inner integration comes ~o,n the pole correction scheme, appik;able to the full infinite range wag~zoidal rule. The nature of the current integrand is such that after the transformation r == u 2 it is an even function of the variable u, which enables us to choose the full infinite range. For general integrands this may not be the case. For example, the present method is not applicable to the one-dimensional Fermi-Dirac integrals of integer order and the Bose-Einstein integrals. However, in all such cases, the integrals can be estimated accurately, by splitting the range of integration and applying IMT like schemes as indicated in [!0|.

N. Mtn~ankumar.A. Natarajan/ ComputerPhysics Communications !01 (1997)47-53

52

Appendix A Some discussion concerning the possible accuracy of the proposed algorithm for the evaluation of the two-dimensional Fcrmi-Dirac inte~al is pertinent here. Recall that the algorithm employs the IMT rule for the outer integration and the trapezoidal rule with pole correction for the inner integration. Error bounds for the two quadrature schemes are known individually. We onl~, need to assess the overall error for the combined application of both the schemes. There are two major sources of error that we need to consider, (0 the error due to the di.~cretisation process of the outer integration assuming the inner integral ~o be known exactly, an,d (ii) the error resulting from the integration of the error terms of d~e ianer integratio~ The discussion emphasi~s how by an appropriate choice of step size and truncation limits, these errors car, be made ncgbgibly small comp?xcd to machine precision limits. In such a situation, the overall error is essentially decided by the round-off errors. We define gw(u, 0), g2(u, 0). #/l(O) and ¢t2(0) below,

As already discussed in the main paper, choosing h < d / 6 = ¢r/6 would make the discretisation error term Ea(O) smaller than the machine precision limit, say y. A bound on the truncation error may be obtained as follows. Note thin

t e l ( u , O)

g2(t',

g , ( u , O ) = [I + e " ' " ~ ° ' - " ] ,

(A.2)

q,,( 0 ) = (cos( 0 ) )" (sin( 0 ) ) ", =

¢,:(0} =

u ~ ' * ~)+-~

f _.g,(,,.

(A.3)

du

o) g.(,,, o )

( A.4)

With these definitions, the integral becomes

--/,,.j2q,,(o) ¢'2(o)

(A.5)

dO.

"w/4

Replacing the inner integral by the truncated trapezoidal sum of N I terms plus the pole correction terms as in Eq. (25) results in two error terms, viz. disc~etisation error Ed(0) (which is similar to Ed(0) as given by Eqs. (14) and (15)) and the series truncation error E~(0) given by E,(o) = 2h z-Nl+l

g,(,,,,

o)

(A.7)

The truncated trapezoidal sum can be replaced approximately as

[

E,(O) --2

~c~2( t'+

u,+ 3 C • ~: +2,. du,

u o = hN i. (A.8)

Define 12j ~ a as follows: x

l,.j+| -- fu 1121+1 e-u: dl~.

(A.9)

o

Integration by parts leads to the following recursion: I _-u z

(A.I)

~e -":[c°~(O)+'in(e)l+2a _< c-(,.,')+ 2.

!,1+1 = -ic

g t ( u , O)-" [I + e"" c ° " ° ' - ° ] ,

0)]-'

"j

u 5 + j l 2 j _ ,.

(A.10)

For j equal to - ½, the inequality be~o~ is true,

io= {.% e-" a

du<~

(A.II) /40

We can combine the recursion and the last inequality to get the required bound in closed form. As an $ example, for the case (v + / z ) equal to ~ or alter5 nately for j equal to ,7, the explicit bound on Et(0) can be obtained as E t ( O ) < 2e2a-'d[½uos +

Tg

• (A.12)

A choice of u o ~ 8 in the case of negative a and ( 2 a + 64) i/2 for positive a is adequate to make E, much less than y. The error E;(O) introduced due to the approximation of me inner integral is given by

E,(O)= ['/2 ,(O)[E,(O) +Ed(O)] dO. "Ir/4

g.(,,, o) (A.6)

(A.13)

IV. Mohankumar. A. Natarajan / Computer Physics Communications I01 (1997) 47-53

By taking the maximum of the sum of the error terms out of the integral and replacing the lower limit of integration by zero, we can get the bound on

the E~( O) as !

E l ( 0 ) < Max[ E,(0) + Ed( 0)] ~ B ( ( t,+ ! ) / 2 . (/~ + I ) / 2 ) ,

(A.14)

where B denotes r|~ beta func:,ion. Evidently, Et is less than % Next, we consider the first source of error, i.e the discretisation error due te the approximation of the outer integral by IMT quadratme as N

n..= E

+ R,(o,) +

)=J

One can easily show that 1;"

.h(02) = - - - c o t ( 0 j ) . gl

(^.!7)

With O~ as ¢r/4 and a large a which is of t ~ order of 103 , 02 is of the order of 10 -3 , which is several orders larger than the distance of the saddle point. For 0 t close to the ¢r/2 end, the pole can become closer to the saddle point. However, a simple evaluation of the residue of this pole shows that it rapidly approaches zero, s i ~ e the pole term con,'ains e m~#,) in the denom)nator. Hence we conclu,k that the poles of g:(u, 0) pose no problem, ,% similar conclusion can be arrived at for the poles of the function g2(u, 0).

(A.15) wi~re % is the weight functie- Using Eqs. (7) and (8) with a and /3 taken as zero and N equ,~! to 128, one, obtains easily the asymptotic error estimate to be less than 3' hcrc again, Some remarks concerning the a~.ymptotic error estimate are necessary here with respect to the present problem. Formulae (7) and (8) are based on the assumption that the dominant contribution to the integral expression for the error comes from the saddle points located near the ends of the interval in the complex 0 plane [7].This holds good if the poles of the intcgrand (other than the ehd points of the interval) lie sufficiently above the contour. For a 128 point I M T quadrature, the saddle point lies at a distance which is of t~ order of I0-s above the real axis in the 0 plane. Since the terms ~t(8), RI(0) and R,(:9) have no poles except at the end~ of interval, we need to cow,sider only the poles of the individual terms of the trapezoidal sum S~,~(0), or equivalently the zeros of the functions gi(u, 0) and g2(u, 0). In particular, we need to examine critically the nearest zeros. Consider g~(u. 0) first. Denoting the pole O by (01 + iO2). Lhe .-elation fixing the nearest pole is given by 2 COS( 0 t + i 0 2 ) "-

53

a

+ ifr,

(A.16)

Since the discretisation and truncation errors are made smaller than the machine precision~ the accuracy of the calculation is limited only by the round-off errors. The total number of summations is typically of the order of 10 5. Assuming an error of 10-16 for each step, the worst possible round-off error is of the order of 10-it. Hence we can safely trust 10 digits. This error estimate is v~,~y conservative, the actual error being much less.

References [ I ] B. Pichon, Comput. Phys. Conmmn. 55 (I 989) 127. [2] R.P. Sagar, Comput. Phys. Commun. 66 (1991) 271. [3] L.W. Fullerton and R.A. Rinker, Comput. Phys. Commm, 39 (1986) 181. [4] M. ~ , Phys. ~(ev. 170 (1968) 3u6. [51 N. Molumkumar a ~ A. N.atarajan. Phys. Status Solidi b 188 (1995) 635. [6] N. Mohankumar and A. Nalatajan, Aslrophys. J. 458 (1996) 233. [7] M. Iri. S. Motiguti and Y. Takasawa, J. Compm. A l ~ . Math. 17 (t98~) 3. [8] B. Bialecki, BIT 29 (1989) 464. [9] J. McNam~, Math. Cornpm. 18 (1964) 368. [10] A. Natam)an and N. Mohankumat. ComptN. Phys, Commmt. 76 (1993) 48.