[R] Solution to differential equation

Fri Dec 17 18:45:39 CET 2010

I computed the solution using 3 different methods: (1) power-series, (2)
hypergeometric function, and (3) quadrature using `integrate'.

All three of them give same results for 0 < k2 < 1.  However, for k2 > 1,
the hypergeometric method does not work, but both quadrature and
power-series methods do well.  For integer values of `k2' only the
quadrature method works.  So, the bottom line is that the quadrature method
might be the best option since it works for all positive values of `k2' and
is also quite fast.

My code is attached in the text file.

Hope this is helpful,
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 [mailto:r-help-bounces at r-project.org] On
Behalf Of Mike Marchywka
Sent: Friday, December 17, 2010 6:02 AM
To: r-help at r-project.org; msamtani at gmail.com
Subject: Re: [R] Solution to differential equation

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

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