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 [mailto: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 [mailto: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 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. ______________________________________________ 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.