[R] Solution to differential equation

[Ravi Varadhan]

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.

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 <rvaradhan at jhmi.edu>
Date: Thursday, December 16, 2010 4:11 pm
Subject: RE: [R] Solution to differential equation
To: 'Ravi Varadhan' <rvaradhan at jhmi.edu>, 'Scionforbai' <scionforbai at gmail.com>, 'mahesh samtani' <msamtani at gmail.com>
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.
>