Michael Meyer <mjhmeyer <at> yahoo.com>
writes:>
Check your logic. The following lines show that integrate *does* return the
correct values:
a = 0.08 # alpha
M <- function(j,s){ return(exp(-j*a*s)) }
A <- matrix(NA, 5, 5)
for (i in 1:5) {
for (j in i:5) {
f <- function(s) { return(M(i,s)) }
g <- function(s) { return(M(j,s)) }
integrand <- function(s){ return(f(s)*g(s)) }
A[i, j] <- integrate(integrand,lower=0,upper=Inf)$value
if (i != j) A[j, i] <- A[i, j]
}
}
A
# [,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
Try setting "A <- matrix(NA, 5, 5)". You'll see that the wrong
entries in
matrix A are still NA, that is not yet computed.
Regards, Hans Werner
> 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()
>