Solve ODE that depends on x in R?

Question

Consider the following ODE with boundary condition

y'[x] = (a - 1)*y[x]/x

y[1] = a*b

where a > 0, b > 0, and x > 0.

I want to solve this ODE numerically in R using the command ode from the R package deSolve.

This is, I want to calculate the solution y[x] for x>0

Questions: (i) I do not know how to include "x" in the denominator in the model equation. (ii) I also don't know how to specify the initial condition.

Disclaimer: I know this equation has an analytical solution, but I am trying to compare it with the numerical solution.

So far, I have unsuccessfully tried:

library(deSolve)

difeq <- function(x, y, parms) {

  a <- parms[1]
  b <- parms[2]
  
  # model equations
  dy <-  y*(a - 1)/***x*** # HOW TO SPECIFY x?
  
  # result  
  return( list(dy) )

}

params  <- c(a = 2, b = 1)
init      <- c(y = params[1]*params[2]) # HOW TO SPECIFY THE INITIAL CONDITION?
times  <- seq(0,5, by = 0.1)
out <- ode(init, times, difeq, params)

plot(out)

Answer 1

Another part of the question was, how to specify the initial condition at x>0. As this is a somewhat additional question to the technical part how to include x in the denominator, I will give it as separate answer.

The approach uses a backward simulation. As negative time steps are not supported in deSolve a workaround can be used:

• Instead of running the simulation from 1 to 0 backwards, we run it forward from -1 to 0.
• Because of a pole at x=0, we use a small value close to zero.
• In the model, we change the sign of the derivative, and as x appears also in the model, we multiply it with the sign too.
• Note: this is explicitly included here to show the general approach for didactical reasons. In the special case here, signs will of course cancel out.

In the script below, we first start integration from the known initial value at x=1. Then we do the "backward" integration to estimate the initial value close to zero. Finally we run the model with the new initial value for the whole range of x and compare the results.

library(deSolve)

difeq <- function(x, y, parms, sign = 1) {
  with(as.list(parms),{
    dy <-  sign * y * (a - 1) / (sign * x)
    return(list(y=dy))
  })
}

## initial simulation starting at x = 1
params  <- c(a = 2, b = 1, sign = 1)
init    <- c(y = params[["a"]] * params[["b"]])
times   <- seq(1, 5, by = 0.1)
out1    <- ode(init, times, difeq, params)

## backwards simulation
close_to_zero <- 1e-16 # numerical precision is limited, do not decrease this more
times  <-  c(-1, -close_to_zero)
out0 <- ode(init, times, difeq, params, sign = -1)

## forward simulation starting close to zero
times  <- c(close_to_zero, seq(0.1, 10, 0.1))
init   <- c(y = out0[[nrow(out0), 2]])
cat("y[", close_to_zero,"] =", init, "\n")
out2   <- ode(init, times, difeq, params)#, hmax=0.01)

## comparison of the initial simulation with an extended time period
plot(out2, out1, lwd=c(1, 5), lty=c("solid", "dotted"), xlab="x")

Answer 2

Looking at the equation, I would say that x is equivalent to time t. Package deSolve uses time in the denominator because this is quite common, but it is not limited to time dependent systems. It can also be something else, e.g. a spatial component, just set x=t, that works exactly the same.

Note also to avoid 0 in the denominator.

library(deSolve)

difeq <- function(x, y, parms) {

  a <- parms[1]
  b <- parms[2]
  
  # model equations
  dy <-  y * (a - 1) / x
  
  # result  
  return(list(y=dy))
}

params <- c(a = 2, b = 1)
init   <- c(y = unname(params[1]*params[2]))
times  <- seq(1, 5, by = 0.1)
out    <- ode(init, times, difeq, params)

## just rename "time" to "x"
colnames(out)[1] <- "x"
head(out)

To get the initial value at x -> 0 one may use an optimizer, run the system backwards (see separate answer) or use another solver (see CRAN Task view).

   time    y
1   1.0  2.0
2   1.1  2.2
3   1.2  2.4
4   1.3  2.6
5   1.4  2.8
6   1.5  3.0
7   1.6  3.2
8   1.7  3.4