Solve ODE that depends on x in R?-CodePudding

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)

CodePudding user response：

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(dy) )

}

params  <- c(a = 2, b = 1)
init      <- c(y = params[1]*params[2]) # this is correct
times  <- seq(1, 5, by = 0.1) # start with value > 0
out <- ode(init, times, difeq, params)

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

head(out)

       x     y.a
[1,] 1.0 2.00000
[2,] 1.1 2.20000
[3,] 1.2 2.42000
[4,] 1.3 2.66200
[5,] 1.4 2.92820
[6,] 1.5 3.22102