Solving complex integrals in R

Question

I would like to solve a double integral in R of the following form:

where:

b0= function(time){0.05*sin(0.1*time)+0.14}
bi=function(time){0.05*sin(0.1*time)+0.12}

It is a quite complicated equation and I'm not sure if I am solving it in the proper way.

Here is what I`ve tried:

I`ve divided the equation in small parts:

wi<-Vectorize(function(n,t,bi,di){di-bi(n+t)},"n")

InnerIntegral = function(tau,t,bi){
  0.1*exp(integrate(wi,lower=0,upper=(tau),t=t,bi=bi,di=0.1)$value)
}

Pext<-function(t,bi,Text){
  integrate(Vectorize(InnerIntegral,c('tau')),lower=0,upper=Text-t,t=t, bi=bi)$value
}

 PrIntegral<-Vectorize(function(t,b0,bi,Text){
     b0(t)*(1-Pext(t,bi=bi,Text=Text))
    },"t")


 #For T=100
    T=100
    integrate(Vectorize(PrIntegral,'t'),lower=0,upper=T,b0=b0,bi=bi,Text=100)$value

This gives me -27.77913, but I'm not sure if it is ok. Someone with experience using integrate in R could help me please? Maybe I'm interpreting in the wrong way the equation...

Thanks!

In your last line of code you Vectorize(PrIntegral,'t') that is already vectorized by t by definition. With that change, the code returns -27.77913 for me — Roman Zenka, Feb 12 '18 at 16:26
I count three different variables of integration: t, tau, and n — IRTFM, Feb 12 '18 at 16:59
I would try giving these to Wolfram Alpha to see what it could give you in closed form. It's worth a try. There are three integrals here. Start with the inner one and work your way out. — duffymo, Feb 12 '18 at 18:49
Actually my code give the same results as @RomanZenka, -27.77913, but I was not sure if I was doing it ok or I had to do something more complicated! Thank you very much for answering! — user3436882, Feb 12 '18 at 22:06

score 2 · Answer 1 · answered Feb 12 '18 at 16:28

I cleaned this up and it worked for me. No need to pass bi and b0 functions.

b0 <- function(time){0.05*sin(0.1*time)+0.14}
bi <- function(time){0.05*sin(0.1*time)+0.12}

wi <- Vectorize(function(n,t,di){di-bi(n+t)},"n")

InnerIntegral <- function(tau,t){
   0.1*exp(integrate(wi,lower=0,upper=(tau),t=t,di=0.1)$value)
}

InnerIntegralVec <- Vectorize(InnerIntegral,c('tau'))

Pext <- function(t,Text){
   integrate(InnerIntegralVec,lower=0,upper=Text-t,t=t)$value
}

PrIntegral <- function(t,Text) {
   b0(t)*(1-Pext(t,Text=Text))
}

PrIntegralVec <- Vectorize(PrIntegral, "t")


#For T=100
T=100
integrate(PrIntegralVec,lower=0,upper=T,Text=100)$value

Solving complex integrals in R

1 Answers1