sorry, wanted to CC list hit wrong button no caffeine
> > From: rvarad...@jhmi.edu
> > To: rvarad...@jhmi.edu
> > Date: Thu, 16 Dec 2010 22:37:17 -0500
> > CC: r-help@r-project.org; msamt...@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@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: rvarad...@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@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: rvarad...@jhmi.edu
> > >
> > >
> > > -----Original Message-----
> > > From: r-help-boun...@r-project.org [ On
> > > Behalf Of Ravi Varadhan
> > > Sent: Wednesday, December 15, 2010 5:08 PM
> > > To: 'Scionforbai'; 'mahesh samtani'
> > > Cc: r-help@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: rvarad...@jhmi.edu
> > >
> > >
> > > -----Original Message-----
> > > From: r-help-boun...@r-project.org [ On
> > > Behalf Of Scionforbai
> > > Sent: Wednesday, December 15, 2010 12:28 PM
> > > To: mahesh samtani
> > > Cc: r-help@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@r-project.org mailing list
> > >
> > > PLEASE do read the posting guide
> > > and provide commented, minimal, self-contained, reproducible code.
> > >
> > > ______________________________________________
> > > R-help@r-project.org mailing list
> > >
> > > PLEASE do read the posting guide
> > > and provide commented, minimal, self-contained, reproducible code.
> > >
> >
> > ______________________________________________
> > R-help@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.
>
______________________________________________
R-help@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.
