Advanced Search

Question: How to evaluate a one-dimensional improper integral of an irregularly oscillatory function numerically?

Question:How to evaluate a one-dimensional improper integral of an irregularly oscillatory function numerically?

sursumCorda 637 Maple 2023
integration numeric shill
May 11 2023
1

The highly oscillatory integrand is: 

(sin(x + sqrt(x)) + x*BesselJ(0, x^2))/(1 + x): # int(%, x = 0 .. infinity, numeric); # Note that the upper limit of `x` is NOT finite.

Unfortunately, Maple still cannot return a result in a few minutes. For instance, 

> 

restart;
infolevel[`evalf/int`] := 1:
> 

smartplot((sin(x+sqrt(x))+x*BesselJ(0, x^2))/(1+x))

 
> 

timelimit(0.1e3, evalf(Int((sin(x+sqrt(x))+x*BesselJ(0, x^2))/(1+x), x = 0 .. infinity)))

Control: Entering NAGInt

trying d01amc (nag_1d_quad_inf)
Control: d01amc failed
evalf/int/control: NAG failed result = result
evalf/int/improper: integrating on interval 0 .. infinity
evalf/int/improper: applying transformation x = 1/x
evalf/int/improper: interval is 0 .. 1 for the integrand:

                        /     (1/2)\            /    2\
                     sin\x + x     / + x BesselJ\0, x /
                     ----------------------------------
                                   1 + x               

and interval is 0 .. 1 for the integrand:

                                                /   1 \
                                         BesselJ|0, --|
                        /       (1/2)\          |    2|
                        |1   /1\     |          \   x /
                     sin|- + |-|     | + --------------
                        \x   \x/     /         x       
                     ----------------------------------
                                 /    1\  2            
                                 |1 + -| x             
                                 \    x/               

Control: Entering NAGInt
Control: trying d01ajc (nag_1d_quad_gen)
Control: d01ajc failed
evalf/int/control: NAG failed result = result

evalf/int/control: singularity at left end-point
evalf/int/transform: series contains {x^(1/2), x^(3/2), x^(5/2), x^(7/2), x^(9/2), x^(11/2), x^(13/2), x^(15/2), x^(17/2), x^(19/2), x^(21/2), x^(23/2), x^(25/2)}

evalf/int/singleft: applying transformation x = x^2
evalf/int/singleft: interval is 0 .. 1. for the integrand:

                  /   / 2            \    2        /    4\\  
                2 \sin\x  + csgn(x) x/ + x  BesselJ\0, x // x
                ---------------------------------------------
                                        2                    
                                   1 + x                     

evalf/int/control: Applying simplify/ln, integrand is 2*(sin(x^2+x)+x^2*BesselJ(0,x^4))/(1+x^2)*x
evalf/int/CreateProc: Trying easyproc
evalf/int/CreateProc: Trying makeproc
evalf/int/ccquad: n = 2 integral estimate = .7758263476953
                    n = 6 integral estimate = .8097413335347
evalf/int/ccquad: n = 18 integral estimate = .8097470386638
                                  error = .8097470386638e-11
From ccquad, result = .8097470386638 integrand evals = 19 error = .8097470386638e-11
Control: Entering NAGInt
Control: trying d01ajc (nag_1d_quad_gen)
Control: d01ajc failed
evalf/int/control: NAG failed result = result

evalf/int/control: singularity at left end-point

evalf/int/transform: series contains {1/x^(3/2), x^(1/2), cos(1/x^2+1/4*Pi), sin((x+x^(1/2))/x^(3/2)), sin(1/4*(4+Pi*x^2)/x^2)}
evalf/int/singleft: applying transformation x = x^2
evalf/int/singleft: interval is 0 .. 1. for the integrand:

                   /   /        /1\  \                    \
                   |   |1 + csgn|-| x|                    |
                   |   |        \x/  |  2          /   1 \|
                 2 |sin|-------------| x  + BesselJ|0, --||
                   |   |      2      |             |    4||
                   \   \     x       /             \   x //
                 ------------------------------------------
                                 3 /     2\                
                                x  \1 + x /                

evalf/int/control: Applying simplify/ln, integrand is 2*(sin((1+x)/x^2)*x^2+BesselJ(0,1/x^4))/x^3/(1+x^2)
evalf/int/control: Applying simplify/trig, integrand is (2*sin((1+x)/x^2)*x^2+2*BesselJ(0,1/x^4))/(x^3+x^5)
evalf/int/control: singularity at left end-point

evalf/int/transform: series contains {cos(1/x^4+1/4*Pi), sin((1+x)/x^2), sin(1/4*(4+Pi*x^4)/x^4)}
evalf/int/transform: no transform found
evalf/int/series: integrating on 0 .. .2493468135214 the series:

       2 (%2)     2 (%2) x   /2 (%2)      \  3   /  2 (%2)      \  5   
      --------- - -------- + |------- - %3| x  + |- ------- + %3| x    
        (1/2)       (1/2)    |  (1/2)     |      |    (1/2)     |      
      Pi      x   Pi         \Pi          /      \  Pi          /      
                                                                       
           /       2 (%2)      \  7   /     2 (%2)      \  9           
         + |- %3 + ------- - %4| x  + |%3 - ------- + %4| x            
           |         (1/2)     |      |       (1/2)     |              
           \       Pi          /      \     Pi          /              
                                                                       
                                                                       
           /2 (%2)           \  11   /          2 (%2) \  13           
         + |------- - %4 - %5| x   + |%4 + %5 - -------| x             
           |  (1/2)          |       |            (1/2)|               
           \Pi               /       \          Pi     /               
                                                                       
           /       2 (%2)      \  15   /     2 (%2)      \  17    / 19\
         + |- %5 + ------- + %6| x   + |%5 - ------- - %6| x   + O\x  /
           |         (1/2)     |       |       (1/2)     |             
           \       Pi          /       \     Pi          /             
                                                                       
                                                                       
                      /        4\                                      
             (1/2)    |4 + Pi x |                                      
      %1 := 2      sin|---------|                                      
                      |     4   |                                      
                      \  4 x    /                                      
                    /1 + x\   (1/2)                                    
      %2 := %1 + sin|-----| Pi                                         
                    |  2  |                                            
                    \ x   /                                            
             (1/2)    /1    1   \                                      
            2      cos|-- + - Pi|                                      
                      | 4   4   |                                      
                      \x        /                                      
      %3 := ---------------------                                      
                      (1/2)                                            
                  4 Pi                                                 
                        /        4\                                    
               (1/2)    |4 + Pi x |                                    
            9 2      sin|---------|                                    
                        |     4   |                                    
                        \  4 x    /                                    
      %4 := -----------------------                                    
                       (1/2)                                           
                  64 Pi                                                
                (1/2)    /1    1   \                                   
            53 2      cos|-- + - Pi|                                   
                         | 4   4   |                                   
                         \x        /                                   
      %5 := ------------------------                                   
                        (1/2)                                          
                  512 Pi                                               
                           /        4\                                 
                  (1/2)    |4 + Pi x |                                 
            1371 2      sin|---------|                                 
                           |     4   |                                 
                           \  4 x    /                                 
      %6 := --------------------------                                 
                          (1/2)                                        
                  16384 Pi                                             

evalf/int/CreateProc: Trying easyproc
evalf/int/CreateProc: Trying makeproc

evalf/int/CreateProc: Trying procmake
evalf/int/control: applying double-exponential method
evalhf mode unsuccessful -- retry in software floats
evalf/int/quadexp: applying double-exponential method

Error, (in tools/sign) time expired

  
> 

time() =

3999.500

 (1)

 =

%?

Download oscillatoryInt.mws

I tried to specify another method (as described in evalf/int), but it seems that this doesn't work. Any ideas?
Please Wait...