sorry, wanted to CC list hit wrong button no caffeine
> > From: rvaradhan at jhmi.edu
> > To: rvaradhan at jhmi.edu
> > Date: Thu, 16 Dec 2010 22:37:17 -0500
> > CC: r-help at r-project.org; msamtani at gmail.com
> > Subject: Re: [R] Solution to differential equation
> >
> > One small correction to my previous email:
> >
> > k1 and k2 can be any positive, real numbers, except that `k2'
cannot be an integer. This is not a problem, since the integral can be easily
obtained when `k2' is an integer.
> >
> > Before developing the power series approach, I was thinking that it
might be possible to use incomplete beta function. However, the exponents of the
two factors in the integrand of the beta function have to be (strictly) greater
than -1. In your integrand, one of the exponents, -k2, can be less than -1, and
the other exponent is equal to -1. Hence, the incomplete beta function would not
work, which is why I had to develop the power series approach.
> >
>
> did you see my earlier post with link to wolfram integrator? Where i also
> requested anyone wanting to get rid of a copy of G&R Integral Tables to
> contact me off list since a dog really did eat mine? I think it came up
> with "F" or hypergeometrics. The denominator reduces to something
like
> y*(a-y)^n depending on one choice of variables. Analytic of course is
> in the eye of the beholder but if you can get it in terms of known
> functions presumably you can find existing code to evaluate or
> even better known identities for manipulation.
>
> http://www.mail-archive.com/r-help at r-project.org/msg120242.html
>
>
>
> > Let me know how this works for you.
> >
> > Ravi.
> >
> > ____________________________________________________________________
> >
> > Ravi Varadhan, Ph.D.
> > Assistant Professor,
> > Division of Geriatric Medicine and Gerontology
> > School of Medicine
> > Johns Hopkins University
> >
> > Ph. (410) 502-2619
> > email: rvaradhan at jhmi.edu
> >
> >
> > ----- Original Message -----
> > From: Ravi Varadhan
> > Date: Thursday, December 16, 2010 4:11 pm
> > Subject: RE: [R] Solution to differential equation
> > To: 'Ravi Varadhan' , 'Scionforbai' , 'mahesh
samtani'
> > Cc: r-help at r-project.org
> >
> >
> > > Hi Mahesh,
> > >
> > > Here is a solution to your problem. I have used the power-series
approach.
> > > You can solve for any positive value of k1 and k2 (except that k2
> > > cannot be
> > > unity). I think this is correct, but you can compare this with a
numerical
> > > ODE solver, if you wish.
> > >
> > > pseries <- function(u, k2, eps=1.e-08) {
> > >
> > > fn <- function(x) x * log(u) - log(eps) - log(x)
> > >
> > > N <- ceiling(uniroot(fn, c(1, 1/eps))$root)
> > >
> > > n <- seq(1, N+1)
> > >
> > > u^(-k2) * sum(u^n / (n - k2))
> > >
> > > }
> > >
> > >
> > > logistic.soln <- function(t, k1, k2, Rm, R0) {
> > >
> > > C <- k1 * Rm^k2
> > >
> > > y0 <- pseries(R0/Rm, k2)
> > >
> > > rhs <- C*t + y0
> > >
> > > ff <- function(x, k2) pseries(x, k2) - rhs
> > >
> > > ans <- try(uniroot(ff, c(1.e-03, 1-1.e-03), tol=1.e-06,
k2=k2)$root,
> > > silent=TRUE)
> > >
> > > y <- if (class(ans) !="try-error") ans * Rm else Rm
> > >
> > > y
> > > }
> > >
> > >
> > >
> > > t <- seq(0, 10, length=200)
> > >
> > > # R0 > 0, Rm > R0, k1 > 0, k2 > 0, k2 != 1
> > >
> > > k1 <- 0.5
> > >
> > > k2 <- 1.5
> > >
> > > Rm <- 2
> > >
> > > R0 <- 0.2
> > >
> > > system.time(y <- sapply(t, function(t) logistic.soln(t, k1,
k2, Rm, R0)))
> > >
> > > plot(t, y, type="l")
> > >
> > >
> > > Hoep this helps,
> > > Ravi.
> > >
> > > -------------------------------------------------------
> > > Ravi Varadhan, Ph.D.
> > > Assistant Professor,
> > > Division of Geriatric Medicine and Gerontology School of Medicine
Johns
> > > Hopkins University
> > >
> > > Ph. (410) 502-2619
> > > email: rvaradhan at jhmi.edu
> > >
> > >
> > > -----Original Message-----
> > > From: r-help-bounces at r-project.org [ On
> > > Behalf Of Ravi Varadhan
> > > Sent: Wednesday, December 15, 2010 5:08 PM
> > > To: 'Scionforbai'; 'mahesh samtani'
> > > Cc: r-help at r-project.org
> > > Subject: Re: [R] Solution to differential equation
> > >
> > > This integral is NOT easy. Your solution is wrong. You cannot
integrate
> > > term-by-term when the polynomial is in the *denominator* as it is
here.
> > >
> > > I am not sure if there is a closed-form solution to this. I
remember
> > > you
> > > had posed a problem before that is only slightly different from
this.
> > >
> > >
> > > Previous formulation:
> > >
> > > dR/dt = k1*R*(1-(R/Rmax)^k2); R(0) = Ro
> > >
> > > Current formulation:
> > >
> > > dR/dt = k1*(R^k2)*(1-(R/Rmax)); R(0) = Ro
> > >
> > >
> > > Note that the difference is that the exponent `k2' is in
different places.
> > > Initially I thought that this should not make much difference.
But
> > > appearances of similarity can be quite deceiving. Your previous
formulation
> > > admitted a closed-form solution, which I had shown you before.
But this
> > > current formulation does not have a closed-form solution, at
least
> > > not to my
> > > abilities.
> > >
> > > For the current formulation, you have u^k2 * (1-u) in the
denominator
> > > of the
> > > integrand. This cannot be simplified in terms of partial
fractions.
> > > In the
> > > case of previous formulation, the denominator was u * (1 - u^k2),
which
> > > could be simplified as (1/u) + u^(k2-1)/(1-u^k2). So, we are
stuck,
> > > it
> > > seems.
> > >
> > > We can expand 1/(1 - u) in a power-series expansion, provided |x|
<
> > > 1, and
> > > then integrate each term of the expansion. However, this is not
very
> > > helpful
> > > as I don't know what this series converges to.
> > >
> > > May be I am missing something simple here? Any ideas?
> > >
> > > Ravi.
> > >
> > > -------------------------------------------------------
> > > Ravi Varadhan, Ph.D.
> > > Assistant Professor,
> > > Division of Geriatric Medicine and Gerontology School of Medicine
Johns
> > > Hopkins University
> > >
> > > Ph. (410) 502-2619
> > > email: rvaradhan at jhmi.edu
> > >
> > >
> > > -----Original Message-----
> > > From: r-help-bounces at r-project.org [ On
> > > Behalf Of Scionforbai
> > > Sent: Wednesday, December 15, 2010 12:28 PM
> > > To: mahesh samtani
> > > Cc: r-help at r-project.org
> > > Subject: Re: [R] Solution to differential equation
> > >
> > > > I am trying to find the analytical solution to this
differential equation
> > > >
> > > > dR/dt = k1*(R^k2)*(1-(R/Rmax)); R(0) = Ro
> > >
> > > > If there is an analytial solution to this differential
equation
> > > then it
> > >
> > > It is a polynomial function of R, so just develop the expression
and
> > > when you get the two terms in R (hint: one has exponent k2, the
other
> > > k2+1) you can simply apply the basic integration rule for
polynomials.
> > > It literally doesn't get easier than that.
> > >
> > > the resulting function is:
> > >
> > > k1/(k2+1) * R^(k2+1) - K1/[Rmax*(k2+2)] * R^(k2+2) + Ro
> > >
> > > scion
> > >
> > > ______________________________________________
> > > R-help at r-project.org mailing list
> > >
> > > PLEASE do read the posting guide
> > > and provide commented, minimal, self-contained, reproducible
code.
> > >
> > > ______________________________________________
> > > R-help at r-project.org mailing list
> > >
> > > PLEASE do read the posting guide
> > > and provide commented, minimal, self-contained, reproducible
code.
> > >
> >
> > ______________________________________________
> > R-help at r-project.org 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.
>