[R] numerical integration

[Hans W Borchers]

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()
>