Yazhini February 2016

How to solve this ode equations in R?

I like to solve the following ordinary differential equations in R by Runge kutta method.

dy1 <- (-51.33) * ((1-y[2]) / y[1])

dy2 <- 1.54 * y[1] * (1-y[2]) - 2.14 * y[2]

When y1 becomes zero, then dy1 will become infinity. To avoid this, I need to write R code stating when y[1] becomes less than 0.001, stop y[1] derivative and keeps zero. I pasted the R code below:

yini <- c(1,0)

intabs <- function (t, y, parms) {
  ifelse (y[1] <= 0.01, dy1 <- 0, no)
  dy1 <- -51.33 * ((1-y[2]) / y[1])
  dy2 <- 1.54 * y[1] * (1-y[2]) - 2.14 * y[2]
  list(c(dy1, dy2))
}

times <- seq(from = 0, to = 1, by = .002)

out <- ode (times = times, y=yini, func = intabs, parms = NULL, method = "rk4")

head (out, n=50)

I used ifelse statement to denote y[1] less than or equal to 0.001, then keep dy1 as zero. I'm getting results. But it seems I made some error, which brings the results of dy1 in negative values. I'm very new in writing programmes. Please help if you spot the error..

Answers


LutzL February 2016

Such a jump in the ODE function (if correctly encoded) is poison for any step size controller. The ODE function has to be at least as smooth as the error order of the integration method for step size control to work as intended. It still can fail if the ODE is stiff.

A less invasive method to de-singularize an expression 1/y is to use

y/(1e-6+y*y)

which for y>1e-2 will give a good approximation of 1/y and is smooth everywhere. Change the constant to adapt the region where it deviates from 1/y.

Post Status

Asked in February 2016
Viewed 1,607 times
Voted 9
Answered 1 times

Search




Leave an answer