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]]
