[R] Solution to differential equation

Fri Dec 17 12:01:37 CET 2010



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@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.
>