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.
