Greetings,
 
I encounter a strange problem computing some numerical integrals on [0,oo).
Define
$$
M_j(x)=exp(-jax)
$$
where $a=0.08$. We want to compute the $L^2([0,\infty))$-inner products
$$
A_{ij}:=(M_i,M_j)=\int_0^\infty M_i(x)M_j(x)dx
$$
Analytically we have
$$
A_{ij}=1/(a(i+j)).
$$
In the code below we compute the matrix $A_{i,j}$, $1\leq i,j\leq 5$, numerically
and check against the known analytic values.
 
When I run this code most components of A are correct, but some are zero.
I get the following output:
 
[1] "Dot products, analytic:"
         [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 6.250000 4.166667 3.125000 2.500000 2.083333
[2,] 4.166667 3.125000 2.500000 2.083333 1.785714
[3,] 3.125000 2.500000 2.083333 1.785714 1.562500
[4,] 2.500000 2.083333 1.785714 1.562500 1.388889
[5,] 2.083333 1.785714 1.562500 1.388889 1.250000

[1] "Dot products, numeric:"
         [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 6.250000 4.166667 3.125000 2.500000 2.083333
[2,] 4.166667 3.125000 2.500000 2.083333 0.000000
[3,] 3.125000 2.500000 2.083333 0.000000 0.000000
[4,] 2.500000 2.083333 0.000000 0.000000 0.000000
[5,] 2.083333 1.785714 1.562500 1.388889 1.250000
 
 
Undoubtedly the code is suboptimal.
What is wrong with it?
 
a  = 0.08  # alpha
M <- function(j,s){ return(exp(-j*a*s)) }
# Inner product in $L^2([0,+\oo))$
#
innerProduct <- function(f,g){
 
     integrand <- function(s){ return(f(s)*g(s)) }
     return(integrate(integrand,lower=0,upper=Inf)$value)
}

# The 5x5 matrix of the L^2-inner products of $M_1,\dots,M_5$.
#
numericalDotProductMatrix_M <- function(){
     A <- matrix(0,5,5)
         for(i in seq(1:5))
              for(j in seq(i:5)){
  
                       f  <- function(s){ return(M_j(i,s)) }
                       g <- function(s){ return(M_j(j,s)) }
   
                       A[i,j] <- innerProduct(f,g)
                       if(i<j) A[j,i] <- A[i,j]
              }
     return(A)
}

# The 5x5 matrix of the L^2-inner products of $M_1,\dots,M_5$.
#
dotProductMatrix_M <- function(){
     
        A <- matrix(0,5,5)
         for(i in seq(1:5))
              for(j in seq(1:5)) A[i,j] <- 1/(a*(i+j))
         return(A)
}

testNumericalIntegration <- function(){
 
     A <- dotProductMatrix_M()
     print("Dot products, analytic:")
     print(A)
 
    # why doesn't the numerical integration work
     B <- numericalDotProductMatrix_M()
     print("Dot products, numeric:")
     print(B)
}

testNumericalIntegration()
	[[alternative HTML version deleted]]

