[R] Problem to solve an integral equation

[David Winsemius]

On Dec 17, 2013, at 11:25 AM, Aurélien Philippot wrote:

> thanks David for your answer.
> 
> Sorry for the integrand; as you noticed, a parenthesis was missing. 
> Here is another integrand which is properly defined:
> 
> integrand<- function(C,x){-((5000000+100000*x+10000*max(0,x-50))^(-1)-(5000000+100000*x+C)^(-1))*
>                             1/(x*0.20*sqrt(10)*sqrt(2*pi))*exp(-0.5*((log(x/50)-0.10*10)/(0.20*sqrt(10)))^2)}
> 
> 
> In the general line of the solution that you propose, I guess you start by generating a sequence of values for C for which you evaluate the objective function. Do you then generate a new sequence between any two numbers of the first sequence for which the output function changes sign? 
> 

As I said, the usual route is to use sapply or Vectorize and it is up to you to provide an approportiate domain,

Vout <- Vectorize(output)

# Your function does not appear to have a solution in the domain you are investigating:

which ( Vout(seq(0, 1000000, by=10000)  ) < 0  )
integer(0)

# So the output of uniroot.all seems appropriate:

uniroot.all(Vout, c(0,1000000))
numeric(0)

-- 
David.
> 
> 
> 
> 
> 2013/12/17 David Winsemius <dwinsemius at comcast.net>
> 
> On Dec 17, 2013, at 8:53 AM, Aurélien Philippot wrote:
> 
> > Dear R experts,
> >
> > I am trying to find numerical solutions for an integral equation.
> >
> >
> > Here is an example:
> >
> > I started by defining the integrand, as a function of x and C, where x is
> > the variable of integration and C is the parameter I am interested in:
> >
> >
> > integrand<- function(C,x){-((100000*x)^(-1)-(1000*x+C)^(-1))*
> > 1/(0.20*sqrt(2*pi))*exp(-0.5*(x-0.10)/(0.20))^2)}
> 
> Error: unexpected ')' in:
> "integrand<- function(C,x){-((100000*x)^(-1)-(1000*x+C)^(-1))*
> 1/(0.20*sqrt(2*pi))*exp(-0.5*(x-0.10)/(0.20))^2)"
> 
> Taking out the last parenthesis resulted in a variety of error messages including:
> 
> > output(0)
> Error in integrate(integrand, C = C, lower = 0, upper = Inf) :
>   extremely bad integrand behaviour
> > output(-10)
> Error in integrate(integrand, C = C, lower = 0, upper = Inf) :
>   the integral is probably divergent
> 
> 
> 
> > Then, I integrate with respect to x, and define an objective function
> > (output) with respect to C:
> >
> >
> > output<- function(C) integrate(integrand, C=C, lower=0, upper=Inf)$value
> >
> >
> 
> Generally these problems are solved with the use of something along the lines of
> 
>     sapply( seq(...), output)  # across an appropriate domain
> 
> Or using the `Vectorize` function.
> 
> But since you did not supply tested code for the integrand, I am not proceeding further.
> 
> --
> david.
> 
> 
> > I would like all the values of C (if any) for which the integral is equal
> > to 0.
> >
> > I did the following:
> >
> >
> > library(rootSolve)
> >
> > uniroot.all(output, c(0,1000000))
> >
> >
> > but that does not work ("evaluation of function gave a result of wrong
> > length"). I cannot either plot the function using plot or curve. I don't
> > understand what is going on?
> >
> > Thanks in advance for any advice!
> >
> > Aurelien
> --
> 
> David Winsemius
> Alameda, CA, USA
> 
> 

David Winsemius
Alameda, CA, USA