# Developers Planet

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")

``````

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

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