numerical integration

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

numerical integration

Michael Meyer
Greetings,

Sorry, the last message was sent by mistake! Here it is again:
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(i,s)) }
                       g <- function(s){ return(M(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()


Many thanks,

Michael Meyer

        [[alternative HTML version deleted]]


______________________________________________
[hidden email] 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.
Reply | Threaded
Open this post in threaded view
|

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

______________________________________________
[hidden email] 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.