- Email: [email protected]
A general purpose algorithm for the integration of functions with power-law singularities of arbitrary exponent
Computer Physics Communications ELSEVIER
ComputerPhysics Communications92 (1995) 225-228
A general purpose algorithm for the integration of functions with power-law singularities of arbitrary exponent G.J. Mata, E. Pestana Departamento de Fisica, Universidad Sim6n Bolfvar, Apartado 89000, Caracas 1080A, Venezuela
Received 10 February 1995
Abstract We describe a method for the integration of functions with power-law singularities of arbitrary exponents. The method does not require a priori knowledge of the exponent of the singularity. This fact makes the method useful in situations where such knowledge is not readily available, as might be the case in numerical computations in which the exponent is a nontrivial function of the input parameters. We have found this method particularly useful when applied to Green function calculations of spectral distributions at surfaces and interfaces.
I. Introduction The standard methods used for the numerical evaluation of singular integrals require a detailed knowledge of the limiting form of the integrand as the singularity is approached. This knowledge is used to devise algorithms that extract the singularity or are optimized - as in gaussian integration - for functions with the known limiting form. In some cases one may lack detailed knowledge of the of the limiting form of the integrand. This could happen, for example, when the function to be integrated is the output of a numerical procedure. We have encountered this situation when using Green functions methods in surface and interface problems to obtain properties such as a total energy or a magnetization [1,2]. The calculation of these quantities involves an integral of the appropriate Green function over the energy bands. This calculation is complicated by the presence of band-edge
singularities that generally follow a power-law with an exponent that is a function of the parameters of the problem and of the point in the Brillouin zone. The value of the exponent as a function of the parameters can be calculated only for simple systems. This fact precludes the use, of the standard methods that require a priori knowledge of the exponent. We have developed an algorithm that works whenever the limiting form of the integrand is an integrable power. The specific value of the exponent is not needed and can be found, in fact, as a by-product of the calculation. The method rests upon the fact that near the singularity a sequence of integrals over intervals that are progressively smaller and closer to the singularity has simple scaling properties. Each integral in the sequence can be evaluated by standard methods. The method here described should be useful for other problems involving the repeated evaluation of integrals of functions with power-law singularities of
0010-4655f95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0010-4655(95)00069-0
G.I. Mata, E. Pestana/ ComputerPhysicsCommunications92 (1995) 225-228
226
arbitrary exponent. It might also be generalized for other classes of singularities for which simple scaling sequences are found, In the following section we describe the method in detail. Finally we present the results of a FORT R A N implementation,
where 1 - 2`'-1 a = ~b(0)
b k = ~b'( ~k)
, 1- a 1 -2"-2
(2.9) (2.10)
2- a
and q = 2 '~-1.
2. The m e t h o d
Our method rests on the observation that I k becomes a geometric sequence for large k (Eq. (2.8)). With this in mind we define the sequence
Let f he a function of the form f = ~b(x) x - " ,
(2.1)
where a < 1 and f ( x ) has a continuous derivative in [0, 1]. We want to compute the integral 1 = folf(x) dx,
(2.2)
under conditions in which the function ~b and the exponent a are not known a priori. We can write oo
I=
(2.11)
_ k _ l f ( x ) dx,
(2.3)
or
=
qk
Ik
(2.12)
Ik- 1
Clearly limk_,= qk = q" At this point we already have an algorithm for the computation of the integral I, namely: choose some k such that qk -= q, compute the first k terms in (2.4), and approximate the remainder of the series by a geometric series. We can do much better, however, as there is a way to accelerate the convergence of qk- To this end we define still another sequence: 1 - ot 1 - 2 a-2 ~ ' ( ~ k ) ck = 2 - ot 1 - 2 `'-1 ~b(O)
(2.13)
and write
co
I = ]~ I , , k=O
(2.4)
qk
=
( l+2-k
with
Ck--2Ck-1
)
-7- ----XzT-k+ 1 q1 ± Ck_lZ
(2.14)
Let c = lim k _, ~ c k. Then Ik =
-k-, f ( x ) f2 .k
dx.
(2.5)
We use Taylor's theorem to find Ik=
f2
-k-, [~b(0) + x c h ' ( ~ ( x ) ) l x - ` "
(2.15)
For large enough k:
dx
(2.6)
qt ~ (1 - 2 - k c ) q ,
(2.16)
or
and we use the weighted mean-value theorem to find I k = ~b(0) 2 - ( 1 - ` ' ) k f 1 y-`" d r "1/2 + 2-(2-`')k~b'(ek) f11/2 y - ` ' + ' d y ,
c k - 2 c k_ 1 lim 2 -k = 2-%. k---,~ 1 + C k _ l 2-~+1
qk q ~ 1 - 2-% "
(2.17)
In addition (2.7)
c --- 2 k qk - qk-1 . 2 qk - qk- 1
(2.18)
where ~k is some point in the interval [0, 2-k]. Finally
In the manner of Aitken's accelerating scheme [3] we introduce (2.18) into (2.16) to obtain the sequence:
I, = (a + b k 2 - k ) q k,
qk = 2qk - qk-1,
(2.8)
(2.19)
GJ. Mata, E. Pestana / Computer Physics Communications 92 (1995) 225-228 w h i c h should converge to q faster than qk" At this point w e can write the nth approximation to the integral I as: i(n)=
A s a stopping criterion w e require that I I (n) I(n-1) I < ~, where e is s o m e desired accuracy. This is equivalent to
In
n-i 1 ]~ ik+i n _ k=O 1 -4,
(2.20)
227
Io-1'~.-,
1 -qn
< e.
(2.21)
1 - c~._,
s u b r o u t i n e SINGINT(f, eps, int, ier) c c This routine a p p r o x i m a t e s the integral from x=0 to x=l of c a function f(x) with a p o w e r law s i n g u l a r i t y at x=0. c f(x) c c c c c c c
is of the
derivative eps is int is ier is adquad
over
form f=~(x)x "~ where
~(x)
the t o l e r a n c e to which the integral is evaluated. the c a l c u l a t e d a p p r o x i m a t i o n to the integral. an e r r o r c o n d i t i o n variable. is an a d a p t i v e i n t e g r a t i o n procedure.
real i0, iold, inew, int external f p a r a m e t e r (kmax=50) c I n i t i a l i z e s e q u e n c e of integrals epsi=0.5*eps call a d q u a d ( 0 . 5 , 1.0 , f, epsi, epsi=0.5*epsi call a d q u a d ( 0 . 2 5 , 0.5 , f, epsi, epsi=0.5*epsi call adquad(0.125, 0.25, f, epsi, qold=iold/i0 qnew=inew/iold qt n e w = 2 * q n e w - q o l d k=3 int=i0 + iold + inew/(l-qtnew) s = i0 + iold + inew I0
has a c o n t i n u o u s
[0, i] a n d ~ < i.
i0) iold) inew)
continue iold=inew qold=qnew qtold=qtnew a=2.**(-k-l) b=2.**(-k) epsi=0.5*epsi call adquad(a, b, f, epsi, inew) qnew=inew/iold qtnew:2*qnew-qold int=s+inew/(l-qtnew) s=s+inew del=inew/(l-qtnew)-iold*qtold/(l-qtold) k=k+l if( ( a b s ( d e l ) .ge. e p s ) . a n d . ( k .it. kmax) ) go to if(k .it. ier=0 else ier=l end if return end
kmax)
I0
then
Fig. 1. A FORTRAN subroutine that implements the algorithm to compute an approximation to an integral of the form f = 4p(x)x -~. The parameter f is the name of the function procedure, eps is the required tolerance for the integral, int is the calculated approximation to the integral, and ier is an error condition variable (on output ier = 0 indicates successful completion of the procedure, ier = 1 indicates that the procedure failed to converge).
G.J. Mata, E. Pestana / Computer Physics Communications 92 (1995) 225-228
228
In practice one uses standard integration routines to evaluate the nonsingular integrals I k. Formulae (2.12), (2.19) and (2.20) are then used to evaluate the integral. In the next section we discuss a FORTRAN implementation of this algorithm as well as a representative application.
In Table 1 we present the results for the integral of f(x) = eX(x + 1 / 2 ) / v ~ - evaluated to a tolerance of ,~-----10 - 4 . It is readily seen that the 4~ sequence converges to the exact value (2.11), faster than the qk sequence. Evaluation of f(x) over ten subintervals was required to achieve the result I = 2.718348, which is within the prescribed tolerance from the exact value e = 2.718282.
3. An implementation In Fig. 1 we show a FORTRAN routine that uses the algorithm described above to calculate an approximation to an integral of the form given by (2.2). The integral (2.5), I k, of f(x) over each subinterval is calculated using an adaptive integration procedure adquad [3], to a tolerance ~k = 2k+ 1 ,
(3.1)
so that the error will be bounded by 6. Table 1 Results from the evaluation of the integral of f ( x ) = eX(x + 1 / 2 ) / V~- to a tolerance of e = 10 -4. In successive columns we display the number of the current subinterval, the value of the integral, lk, over the corresponding subinterval, the sequence qk, and the accelerated sequence '~k. The approximation to the total integral I is displayed in the last column. k
Ik
qi
t~k
I
3 4 5 6 7 8 9 10
0.1345089 8.373562E-02 5.541987E-02 3.788698E-02 2.633704E-02 1.846409E-02 1.300008E-02 9.172659E-03
0.5572270 0.6225286 0.6618434 0.6836353 0.6951475 0.7010694 0.7040735 0.7055849
0.6536254 0.6878301 0.7011582 0.7054272 0.7066598 0.7069914 0.7070775 0.7070964
2.706013 2.720425 2.721372 2.719960 2.719014 2.718583 2.718412 2.718348
4. Summary We have presented a method for integrating functions with power-law singularities of arbitrary exponents. The procedure is based upon the fact that the integration interval can be partitioned in such a way that the integral is the sum of a sequence I k which becomes geometric for large k and can be accelerated. The resulting algorithm is automatic in the sense that it does not require a priori knowledge of the exponent of the singularity.
References [1] E. Pestana, G.J. Mata and N. Pantoja, Physica B 165, 166 (1990) 451. [2] G.J. Mata and E. Pestana, Phys. Rev. B 42 (1990) 885. [3] R.L. Burden, J.D. Faires, A.C. Reynolds, Numerical Analysis, (Prindle, Weber and Schmidt, Boston, 1981).
ComputerPhysics Communications92 (1995) 225-228
A general purpose algorithm for the integration of functions with power-law singularities of arbitrary exponent G.J. Mata, E. Pestana Departamento de Fisica, Universidad Sim6n Bolfvar, Apartado 89000, Caracas 1080A, Venezuela
Received 10 February 1995
Abstract We describe a method for the integration of functions with power-law singularities of arbitrary exponents. The method does not require a priori knowledge of the exponent of the singularity. This fact makes the method useful in situations where such knowledge is not readily available, as might be the case in numerical computations in which the exponent is a nontrivial function of the input parameters. We have found this method particularly useful when applied to Green function calculations of spectral distributions at surfaces and interfaces.
I. Introduction The standard methods used for the numerical evaluation of singular integrals require a detailed knowledge of the limiting form of the integrand as the singularity is approached. This knowledge is used to devise algorithms that extract the singularity or are optimized - as in gaussian integration - for functions with the known limiting form. In some cases one may lack detailed knowledge of the of the limiting form of the integrand. This could happen, for example, when the function to be integrated is the output of a numerical procedure. We have encountered this situation when using Green functions methods in surface and interface problems to obtain properties such as a total energy or a magnetization [1,2]. The calculation of these quantities involves an integral of the appropriate Green function over the energy bands. This calculation is complicated by the presence of band-edge
singularities that generally follow a power-law with an exponent that is a function of the parameters of the problem and of the point in the Brillouin zone. The value of the exponent as a function of the parameters can be calculated only for simple systems. This fact precludes the use, of the standard methods that require a priori knowledge of the exponent. We have developed an algorithm that works whenever the limiting form of the integrand is an integrable power. The specific value of the exponent is not needed and can be found, in fact, as a by-product of the calculation. The method rests upon the fact that near the singularity a sequence of integrals over intervals that are progressively smaller and closer to the singularity has simple scaling properties. Each integral in the sequence can be evaluated by standard methods. The method here described should be useful for other problems involving the repeated evaluation of integrals of functions with power-law singularities of
0010-4655f95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0010-4655(95)00069-0
G.I. Mata, E. Pestana/ ComputerPhysicsCommunications92 (1995) 225-228
226
arbitrary exponent. It might also be generalized for other classes of singularities for which simple scaling sequences are found, In the following section we describe the method in detail. Finally we present the results of a FORT R A N implementation,
where 1 - 2`'-1 a = ~b(0)
b k = ~b'( ~k)
, 1- a 1 -2"-2
(2.9) (2.10)
2- a
and q = 2 '~-1.
2. The m e t h o d
Our method rests on the observation that I k becomes a geometric sequence for large k (Eq. (2.8)). With this in mind we define the sequence
Let f he a function of the form f = ~b(x) x - " ,
(2.1)
where a < 1 and f ( x ) has a continuous derivative in [0, 1]. We want to compute the integral 1 = folf(x) dx,
(2.2)
under conditions in which the function ~b and the exponent a are not known a priori. We can write oo
I=
(2.11)
_ k _ l f ( x ) dx,
(2.3)
or
=
qk
Ik
(2.12)
Ik- 1
Clearly limk_,= qk = q" At this point we already have an algorithm for the computation of the integral I, namely: choose some k such that qk -= q, compute the first k terms in (2.4), and approximate the remainder of the series by a geometric series. We can do much better, however, as there is a way to accelerate the convergence of qk- To this end we define still another sequence: 1 - ot 1 - 2 a-2 ~ ' ( ~ k ) ck = 2 - ot 1 - 2 `'-1 ~b(O)
(2.13)
and write
co
I = ]~ I , , k=O
(2.4)
qk
=
( l+2-k
with
Ck--2Ck-1
)
-7- ----XzT-k+ 1 q1 ± Ck_lZ
(2.14)
Let c = lim k _, ~ c k. Then Ik =
-k-, f ( x ) f2 .k
dx.
(2.5)
We use Taylor's theorem to find Ik=
f2
-k-, [~b(0) + x c h ' ( ~ ( x ) ) l x - ` "
(2.15)
For large enough k:
dx
(2.6)
qt ~ (1 - 2 - k c ) q ,
(2.16)
or
and we use the weighted mean-value theorem to find I k = ~b(0) 2 - ( 1 - ` ' ) k f 1 y-`" d r "1/2 + 2-(2-`')k~b'(ek) f11/2 y - ` ' + ' d y ,
c k - 2 c k_ 1 lim 2 -k = 2-%. k---,~ 1 + C k _ l 2-~+1
qk q ~ 1 - 2-% "
(2.17)
In addition (2.7)
c --- 2 k qk - qk-1 . 2 qk - qk- 1
(2.18)
where ~k is some point in the interval [0, 2-k]. Finally
In the manner of Aitken's accelerating scheme [3] we introduce (2.18) into (2.16) to obtain the sequence:
I, = (a + b k 2 - k ) q k,
qk = 2qk - qk-1,
(2.8)
(2.19)
GJ. Mata, E. Pestana / Computer Physics Communications 92 (1995) 225-228 w h i c h should converge to q faster than qk" At this point w e can write the nth approximation to the integral I as: i(n)=
A s a stopping criterion w e require that I I (n) I(n-1) I < ~, where e is s o m e desired accuracy. This is equivalent to
In
n-i 1 ]~ ik+i n _ k=O 1 -4,
(2.20)
227
Io-1'~.-,
1 -qn
< e.
(2.21)
1 - c~._,
s u b r o u t i n e SINGINT(f, eps, int, ier) c c This routine a p p r o x i m a t e s the integral from x=0 to x=l of c a function f(x) with a p o w e r law s i n g u l a r i t y at x=0. c f(x) c c c c c c c
is of the
derivative eps is int is ier is adquad
over
form f=~(x)x "~ where
~(x)
the t o l e r a n c e to which the integral is evaluated. the c a l c u l a t e d a p p r o x i m a t i o n to the integral. an e r r o r c o n d i t i o n variable. is an a d a p t i v e i n t e g r a t i o n procedure.
real i0, iold, inew, int external f p a r a m e t e r (kmax=50) c I n i t i a l i z e s e q u e n c e of integrals epsi=0.5*eps call a d q u a d ( 0 . 5 , 1.0 , f, epsi, epsi=0.5*epsi call a d q u a d ( 0 . 2 5 , 0.5 , f, epsi, epsi=0.5*epsi call adquad(0.125, 0.25, f, epsi, qold=iold/i0 qnew=inew/iold qt n e w = 2 * q n e w - q o l d k=3 int=i0 + iold + inew/(l-qtnew) s = i0 + iold + inew I0
has a c o n t i n u o u s
[0, i] a n d ~ < i.
i0) iold) inew)
continue iold=inew qold=qnew qtold=qtnew a=2.**(-k-l) b=2.**(-k) epsi=0.5*epsi call adquad(a, b, f, epsi, inew) qnew=inew/iold qtnew:2*qnew-qold int=s+inew/(l-qtnew) s=s+inew del=inew/(l-qtnew)-iold*qtold/(l-qtold) k=k+l if( ( a b s ( d e l ) .ge. e p s ) . a n d . ( k .it. kmax) ) go to if(k .it. ier=0 else ier=l end if return end
kmax)
I0
then
Fig. 1. A FORTRAN subroutine that implements the algorithm to compute an approximation to an integral of the form f = 4p(x)x -~. The parameter f is the name of the function procedure, eps is the required tolerance for the integral, int is the calculated approximation to the integral, and ier is an error condition variable (on output ier = 0 indicates successful completion of the procedure, ier = 1 indicates that the procedure failed to converge).
G.J. Mata, E. Pestana / Computer Physics Communications 92 (1995) 225-228
228
In practice one uses standard integration routines to evaluate the nonsingular integrals I k. Formulae (2.12), (2.19) and (2.20) are then used to evaluate the integral. In the next section we discuss a FORTRAN implementation of this algorithm as well as a representative application.
In Table 1 we present the results for the integral of f(x) = eX(x + 1 / 2 ) / v ~ - evaluated to a tolerance of ,~-----10 - 4 . It is readily seen that the 4~ sequence converges to the exact value (2.11), faster than the qk sequence. Evaluation of f(x) over ten subintervals was required to achieve the result I = 2.718348, which is within the prescribed tolerance from the exact value e = 2.718282.
3. An implementation In Fig. 1 we show a FORTRAN routine that uses the algorithm described above to calculate an approximation to an integral of the form given by (2.2). The integral (2.5), I k, of f(x) over each subinterval is calculated using an adaptive integration procedure adquad [3], to a tolerance ~k = 2k+ 1 ,
(3.1)
so that the error will be bounded by 6. Table 1 Results from the evaluation of the integral of f ( x ) = eX(x + 1 / 2 ) / V~- to a tolerance of e = 10 -4. In successive columns we display the number of the current subinterval, the value of the integral, lk, over the corresponding subinterval, the sequence qk, and the accelerated sequence '~k. The approximation to the total integral I is displayed in the last column. k
Ik
qi
t~k
I
3 4 5 6 7 8 9 10
0.1345089 8.373562E-02 5.541987E-02 3.788698E-02 2.633704E-02 1.846409E-02 1.300008E-02 9.172659E-03
0.5572270 0.6225286 0.6618434 0.6836353 0.6951475 0.7010694 0.7040735 0.7055849
0.6536254 0.6878301 0.7011582 0.7054272 0.7066598 0.7069914 0.7070775 0.7070964
2.706013 2.720425 2.721372 2.719960 2.719014 2.718583 2.718412 2.718348
4. Summary We have presented a method for integrating functions with power-law singularities of arbitrary exponents. The procedure is based upon the fact that the integration interval can be partitioned in such a way that the integral is the sum of a sequence I k which becomes geometric for large k and can be accelerated. The resulting algorithm is automatic in the sense that it does not require a priori knowledge of the exponent of the singularity.
References [1] E. Pestana, G.J. Mata and N. Pantoja, Physica B 165, 166 (1990) 451. [2] G.J. Mata and E. Pestana, Phys. Rev. B 42 (1990) 885. [3] R.L. Burden, J.D. Faires, A.C. Reynolds, Numerical Analysis, (Prindle, Weber and Schmidt, Boston, 1981).