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() > ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.