An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20131217/53889974/attachment.pl>
Problem to solve an integral equation
4 messages · Aurélien Philippot, David Winsemius
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
An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20131217/a39ae36b/attachment.pl>
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