Convergence acceleration by the Ford-Sidi W(m) algorithm

From Wikiversity
Jump to navigation Jump to search


Computation of infinite integrals involving non-oscillating or oscillating functions (trigonometric or Bessel) multiplied by Bessel function of an arbitrary order[edit]

Although we deal with infinite integrals we will use transformation for infinite series [2] because these integrals are evaluated by accelerating the convergence of sequence of partial sums. Still using the transformation for infinite integrals [2] produces excellent results.The integrals are of the form : where c(x) is non-oscillating function or oscillating function (trigonometric or Bessel) and w(x) is Bessel  function of an arbitrary order.

These partial sums (as proposed by A. Sidi) are computed by integration between consecutive zeros of the Bessel function (see[6]) by a low order Gaussian quadrature (see[7]), although other numerical quadrature methods can be employed.

The Ford - Sidi algorithm routine is applied to these partial sums for convergence acceleration.

The transformation is implemented very efficiently by the Ford-Sidi algorithm [1] . Special case with m=1 (algorithm) is the Sidi W algorithm [4].

This Ford-Sidi algorithm [1] is cheap mainly because set of sophisticated recursions that enable to compute the g function and the rest of the algorithm only up to the value m and not to the full length of LMAX+1 as needed by direct application of the Ford-Sidi algorithm.

We propose here  Ford-Sidi algorithm computed by direct application of the Ford-Sidi algorithm, even if it costs more

than the Ford-Sidi algorithm [1] because :

1. It is simple and short as it does not include those sophisticated recursions.

2. The values of the g function above m up to LMAX+1 are possible easily using the g(k)=t*g(k-m) (see 3.10 p. 1219 in [1]).

3. Most of the arrays dimensions are reduced to 1-D.

In computing the g function the simplest sampling of r(l)=l is chosen. It enables easily the efficiency of the code by saving previous values of the partial sums. This sampling is also very useful for alternating series.

The proposed Ford-Sidi algorithm produces accurate results even with simplest sampling of  r(l)=l and small value of m=2.

This algorithm copes successfully with:

  1. Integrals with non oscillating function multiplied by Bessel function of arbitrary order.
  2. Weber Schafheitin discontinuous integrals with oscillating  functions as, for example, cos(x) and  cos(x)/x multiplied by Bessel function of order one  with short range of r=2 while most extrapolation algorithms (except Wynn's Epsilon) fail.
  3. Integrals with oscillating function as Bessel multiplied by Bessel function of arbitrary order.
  4. Integrals that involve Bessel function of very high order (see also [6] and the transformation of Sidi [3]).
  5. Integrals as (27)(28) that there integrands are not integrable at infinity and need to be defined in the sense of "Abel sums". (see[5]).

30 integrals with 6 different ranges of r=0.05, 1.0 , 2.0, 4.0, 10.0 and 100.0 are tested successfully.

The proposed Ford-Sidi algorithm recursive scheme is :

Rewriting with mm=n-(k+1) :

Formulating the arrays in  1-D with mm=n-(k+1) :

with S=partial sum Nmax=number of partial sums

the simple and short semi formal code for the Ford - Sidi algorithm is :


 1 // MAIN PROGRAM
 2 
 3 S = 0
 4 
 5 //// calling FUNCTION WMALG that computes the Ford-Sidi W(m) algorithm for convergence acceleration.
 6 // and also calling FUNCTION MLTAG that computes the g funcion and the partial sums 
 7  
 8 
 9 
10 
11 END PROGRAM
12 
13 //Routine for computing Ford-Sidi W(m) algorithm
14 FUNCTION WMALG
15 // prevent zero denominator - not exists in [1]. 
16 IF (ABS(G1) >= 1D-77) THEN
17  FSA(N) = S / G(1)
18  FSI(N) = 1 / G(1)
19  FOR I=2 TO Nmax + 1 DO
20   FSG(I,N) = G(I) / G(1)
21  ENDFOR 
22  ELSE
23  FSA(N) = S 
24  FSI(N) = 1 
25  FOR I=2 TO Nmax + 1 DO
26   FSG(I,N) = G(I) 
27  ENDFOR 
28  ENDIF 
29  FOR K=0 TO N-1 DO
30   MM=N-(K+1)
31   D = FSG(K+2,MM+1) - FSG(K+2,MM)
32   FOR I=K+3 TO Nmax+1 DO
33    FSG(I,MM) = (FSG(I,MM+1) - FSG(I,MM)) / D
34   ENDFOR 
35   FSA(MM) = (FSA(MM+1) - FSA(MM)) / D
36   FSI(MM) = (FSI(MM+1) - FSI(MM)) / D 
37  ENDFOR 
38 ENDFOR
39 // prevent zero denominator
40 IF (ABS(FSA(0)) >= 1D-77) THEN
41 APPROX = FSA(0) / FSI(0)
42 ELSE
43 APPROX = 1D-77
44 ENDIF
45 
46 END FUNCTION

with AN = sequence element S = partial sum   the simple and short semi formal code for computing the partial sums and the g function is :

//routine for computing partial sums and g function
FUNCTION MLTAG
 // making partial sums S  
 //the simplicity of this part of the code is enabled because the choice of the simplest sampling
 //R(L)=L (see B.5 p. 1228 in [1]). It also enables efficiency of the code by saving previous values
 //of the partial sums (S).
 AN=
 S =S+AN
 
 // making g function up to M
 FOR K=1 TO M DO
//computing the sequence elements AN(needed for the g function)

//efficiency of this part of the code can be also achived by saving previous results of the sequence elements
//AN that were computed already for the prtial sums S.
  AN=  G(K)
 ENDFOR

 // forward difference
 FOR I=2 TO M DO
  FOR J=M TO I STEP=-1 DO
   G(J)=G(J) - G(J-1)
  ENDFOR 
 ENDFOR

 T=1/(N +1)

 // making g function from M to Nmax+1 (see 3.10 p. 1219 in [1]).
 FOR K=1 TO Nmax+1 DO
  IF K <= M THEN
   G(K)=G(K) * ((N+1) ^ k)
  ELSE
   G(K)=T * G(K-M)
  ENDIF
 ENDFOR 
ENDFOR

END FUNCTION

where for integrals with complex functions (1), (4) and (6)

The 30 tested integrals, with exact values, are:[edit]

(1)   rapidly convergent - complex function

(2) rapidly convergent

(3) slowly convergent

(4) slowly convergent - complex function

(5) algebraically divergent

(6) algebraically divergent - complex function

(7) oscillatory kernel

(8) oscillatory kernel

(9) Bessel function order greater than 1

(10) products of Bessel functions

(11) products of Bessel functions

(12) products of Bessel functions

(13) products of Bessel functions

(14) rapidly convergent

(15) rapidly convergent

Integrals with r = 1.0

(16) divergent oscillatory (see[5][7])

(17) (see[7])

(18)

(19)

(20) big order of Bessel function

(21) very big order of Bessel function

(22) oscillatory kernel - product of Bessel functions

(23) oscillatory kernel - product of Bessel functions (see [3])

(24) see [6]

(25) big order of Bessel function (see [6])

(26) very big order of Bessel function (see [6])

(27) (see[3])

(28) divergent oscillatory (see[5][7])

(29) very oscillatory (see[8])

(30) (see[3])


Computed results of the integrals compared with exact values - Excellent accuracy[edit]

IN=Integral No.  R=Radius(Range)  ACCUR=Requested Accuracy COMPR=Real Computed Result COMPI=Imaginary Computed Result EXACTR=Real Exact Result  EXACTI=Imaginary Exact Resu lt

IN    R      ACCUR         COMPR            COMPI         EXACTR                EXACTI     
 1   0.05   0.1E-09     0.3535533216    -0.3532409596    0.3535533216       -0.3532409596       
 2   0.05   0.1E-09    0.2495322244E-01  0.000000000     0.2495322244E-01   0.000000000   
 3   0.05   0.1E-09     20.00000000      0.000000000     20.00000000            0.000000000   
 4   0.05   0.1E-09     19.29318268     -0.6824013726    19.29318268            -0.6824013726     
 5   0.05   0.1E-09   -0.1558940106E-10   0.000000000    0.000000000           0.000000000      
 6   0.05   0.1E-09    -7999.770489       9.764355802    -7999.770489          9.764355802      
 7   0.05   0.1E-09   -0.2504697287E-01   0.000000000    -0.2504697287E-01   0.000000000      
 8   0.05   0.1E-09    0.1681838313E-14   0.000000000     0.000000000            0.000000000      
 9   0.05   0.1E-09    0.6234411536E-03   0.000000000     0.6234411536E-03    0.000000000      
10   0.05   0.1E-09   -0.2336820089E-13   0.000000000     0.000000000           0.000000000      
11   0.05   0.1E-09   -0.1966156371E-14   0.000000000     0.000000000           0.000000000      
12   0.05   0.1E-09   -0.1925418594E-26   0.000000000     0.9196825005E-10   0.000000000      
13   0.05   0.1E-09    0.3861415026E-04   0.000000000     0.3861415026E-04   0.000000000
14   0.05  0.1E-09  0.1851080515E-02  0.000000000   0.1851080515E-02  0.000000000  
15   0.05  0.1E-09  0.7401237782E-01 000000000      0.7401237782E-01  0.000000000          


 1   2.00   0.1E-09     0.2457791604     -0.1928180249E-01   0.2457791604     -0.1928180249E-01  
 2   2.00   0.1E-09     0.2763932023       0.000000000       0.2763932023      0.000000000      
 3   2.00   0.1E-09     0.5000000000       0.000000000       0.5000000000      0.000000000  
 4   2.00   0.1E-09    0.1895626091E-01   -0.1200712156      0.1895626091E-01  -0.1200712156      
 5   2.00   0.1E-09     0.9174218254E-13   0.000000000       0.000000000       0.000000000    
 6   2.00   0.1E-09   -0.5389270093E-01   0.6576733896E-01  -0.5389270093E-01  0.6576733896E-01  
 7   2.00   0.1E-09     0.5000000000       0.000000000       0.5000000000      0.000000000     
 8   2.00   0.1E-09     0.8660254038       0.000000000       0.8660254038      0.000000000     
 9   2.00   0.1E-09     0.1708203932       0.000000000       0.1708203932      0.000000000     
10   2.00   0.1E-09     0.2500000000       0.000000000       0.2500000000      0.000000000     
11   2.00   0.1E-09     0.6250000000E-01   0.000000000       0.6250000000E-01  0.000000000     
12   2.00   0.1E-09     0.5731056602E-05   0.000000000       0.5731052806E-05  0.000000000     
13   2.00   0.1E-09     0.4696496560E-04   0.000000000       0.4696496450E-04  0.000000000     
14   2.00  0.1E-09  0.4266924586E-01  0.000000000   0.4266924586E-01  0.000000000                         
15   2.00  0.1E-09  0.2297574777E-01  0.000000000   0.2297574777E-01  0.000000000     


 1   4.00   0.1E-09    -0.1344288395E-01  0.2631845721E-01  -0.1344288395E-01  0.2631845721E-01  
 2   4.00   0.1E-09     0.1893660937       0.000000000       0.1893660937       0.000000000      
 3   4.00   0.1E-09     0.2500000000       0.000000000       0.2500000000       0.000000000      
 4   4.00   0.1E-09    -0.1405775698E-01  -0.4552202582E-02 -0.1405775698E-01  -0.4552202582E-02  
 5   4.00   0.1E-09     0.2293554563E-13   0.000000000       0.000000000        0.000000000     
 6   4.00   0.1E-09     0.2558970306E-02   0.3574319813E-02  0.2558970306E-02   0.3574319813E-02  
 7   4.00   0.1E-09     0.2500000000       0.000000000       0.2500000000       0.000000000      
 8   4.00   0.1E-09     0.9682458366       0.000000000       0.9682458366       0.000000000      
 9   4.00   0.1E-09     0.1478525782       0.000000000       0.1478525782       0.000000000      
10   4.00   0.1E-09     0.2500000000       0.000000000       0.2500000000       0.000000000      
11   4.00   0.1E-09     0.3906250000E-02   0.000000000       0.3906250000E-02   0.000000000      
12   4.00   0.1E-09     0.4230955773E-04   0.000000000       0.4230954559E-04   0.000000000      
13   4.00   0.1E-09     0.7499947354E-04   0.000000000       0.7499947370E-04   0.000000000    
14  4.00  0.1E-09  0.3200000000E-01  0.000000000   0.3200000000E-01  0.000000000                    
15  4.00  0.1E-09  0.6400000000E-03  0.000000000   0.6400000000E-03  0.000000000                       


 1   10.00  0.1E-09   -0.3962101630E-08  -0.9735953642E-08  -0.3962101582E-08  -0.9735953664E-08  
 2   10.00  0.1E-09    0.9004962810E-01   0.000000000        0.9004962810E-01   0.000000000      
 3   10.00  0.1E-09    0.1000000000       0.000000000        0.1000000000       0.000000000 
 4   10.00  0.1E-09    0.5990701076E-04  -0.6020541162E-04   0.5990701076E-04  -0.6020541162E-04  
 5   10.00  0.1E-09    0.3664870892E-14   0.000000000        0.000000000        0.000000000    
 6   10.00  0.1E-09   -0.9092300944E-05   0.6231542441E-06  -0.9092300945E-05   0.6231542437E-06  
 7   10.00  0.1E-09    0.1000000000       0.000000000        0.1000000000       0.000000000      
 8   10.00  0.1E-09    0.9949874371       0.000000000        0.9949874371       0.000000000      
 9   10.00  0.1E-09    0.8149379340E-01   0.000000000        0.8149379340E-01   0.000000000      
10   10.00  0.1E-09    0.1000000000       0.000000000        0.1000000000       0.000000000      
11   10.00  0.1E-09    0.9999999998E-04   0.000000000        0.1000000000E-03   0.000000000      
12   10.00  0.1E-09    0.3729561332E-03   0.000000000        0.3729561345E-03   0.000000000      
13   10.00  0.1E-09    0.3869011865E-03   0.000000000        0.3869011875E-03   0.000000000    
14  10.00  0.1E-09   0.8787397112E-02  0.000000000  0.8787397112E-02   0.000000000  
15  10.00  0.1E-09  -0.6610702415E-03  0.000000000  -0.6610702415E-03  0.000000000 


 1   100.00  0.1E-09   -0.2222804096E-13  0.8523120307E-13   0.000000000        0.000000000      
 2   100.00  0.1E-09   0.9900005000E-02   0.000000000        0.9900005000E-02   0.000000000     
 3   100.00  0.1E-09   0.1000000000E-01   0.000000000        0.1000000000E-01   0.000000000 
 4   100.00  0.1E-09   0.2980559246E-12  -0.3626946695E-12  -0.4851871203E-34  -0.1952579141E-32  
 5   100.00  0.1E-09   0.1216602506E-15   0.000000000        0.000000000        0.000000000     
 6   100.00  0.1E-09   -0.7847883994E-17 -0.3434899908E-17  -0.1345888854E-34   0.1434515653E-34 
 7   100.00  0.1E-09   0.1000000000E-01   0.000000000        0.1000000000E-01   0.000000000     
 8   100.00  0.1E-09   0.9999499987       0.000000000        0.9999499987       0.000000000     
 9   100.00  0.1E-09   0.9801499937E-02   0.000000000        0.9801499938E-02   0.000000000     
10   100.00  0.1E-09   0.1000000000E-01   0.000000000        0.1000000000E-01   0.000000000     
11   100.00  0.1E-09   0.9999978411E-08   0.000000000        0.1000000000E-07   0.000000000     
12   100.00  0.1E-09   0.6601490111E-13   0.000000000        0.2152044132E-37   0.000000000      
13   100.00  0.1E-09  -0.4478087354E-14   0.000000000        0.4211988247E-27   0.000000000 
14  100.00  0.1E-09  0.9986515138E-04  0.000000000   0.9986515172E-04  0.000000000    
15  100.00  0.1E-09  -0.9959575539E-06  0.000000000 -0.9959575826E-06  0.000000000    

     
     
16    1.00   0.1E-09    -1.000000000      0.000000000        -1.000000000        0.000000000      
17    1.00   0.1E-09    0.4210244382      0.000000000         0.4210244382       0.000000000      
18    1.00   0.1E-09    0.4210244382      0.000000000         0.4210244382       0.000000000      
19    1.00   0.1E-09     1.000000000      0.000000000         1.000000000        0.000000000      
20    1.00   0.1E-09    0.9897054531E-01  0.000000000         0.9897054531E-01   0.000000000      
21    1.00   0.1E-09    0.9998998075E-02  0.000000000         0.9998999700E-02   0.000000000      
22    1.00   0.1E-09    0.6366197724      0.000000000         0.6366197724       0.000000000      
23    1.00   0.1E-09    0.2450350646      0.000000000         0.2450350646       0.000000000     
24    1.00    0.1E-09    0.2596307983        0.000000000      0.2596308154        0.000000000     
25    1.00    0.1E-09    0.9266646414E-01    0.000000000      0.9266646571E-01    0.000000000    
26    1.00    0.1E-09    0.9992004781E-02    0.000000000      0.9607746870E-02    0.000000000      
27    1.00    0.1E-09    1.0000000000    0.000000000      1.0000000000    0.000000000      
28    1.00    0.1E-09    9.0000000000    0.000000000      9.0000000000    0.000000000     
29    1.00    0.1E-09    2.637633938    0.000000000      2.627160401844    0.000000000     
30    1.00    0.1E-09  0.2630108981E-01 0.000000000      0.2660899796E-01 0.000000000      

References[edit]

[1] W.F. Ford and A. Sidi. An algorithm for a generalization of the Richardson extrapolation process. SIAM J. Numer. Anal., 24:1212–1232, 1987.

[2] D. Levin and A. Sidi. Two new classes of nonlinear transformations for accelerating the convergence of infinite integrals and series. Appl. Math. Comp., 9:175–215, 1981. Originally appeared as a Tel Aviv University preprint in 1975.

[3] A. Sidi. Extrapolation methods for oscillatory infinite integrals. J. Inst. Maths. Applics., 26:1–20, 1980.

[4] A. Sidi. An algorithm for a special case of a generalization of the Richardson extrapolation process. Numer. Math., 38:299–307, 1982.

[5] A. Sidi. Extrapolation methods for divergent oscillatory infinite integrals that are defined in the sense of summability. J. Comp. Appl. Math., 17:105–114, 1987.

[6] A. Sidi. Computation of infinite integrals involving Bessel functions of arbitrary order by the D¯-transformation. J. Comp. Appl. Math., 78:125–130, 1997.

[7] A. Sidi. A User-Friendly extrapolation method for oscillatory infinite integrals. Math. Comp., 51,249–266, 1988.

[8] A. Sidi. The numerical evaluation of very oscillatory infinite integrals by extapolation. Math. Comp., 38, 517-529, 1982.