Search code examples
rdifferential-equations

Warning message on Euler Method (Differential Equations) in R


I am studying differential equations and the Euler method and can't find the bug in my program.

The idea is to simulate estimations of the differential equations using Euler's method.

Here's some code we tried in class.

# Define differential equation
y.prime.diff.eqn <- function(p, y) {return(5*y-p)}


initial.condition.x1 <- 0
initial.condition.y1 <- 0

# Define Euler Method's estimations
euler <-  function(x1 = initial.condition.x1,
                   y1 = initial.condition.y1,
                   y.prime = y.prime.diff.eqn(x1, y1),
                   iter = 5,
                   step.size = 1) {

  for (i in 2:iter)

  {

    x1[i]     <-  x1[i-1] + step.size
    y1[i]     <-  y1[i-1] + step.size * (y.prime)
    y.prime   <-  y.prime.diff.eqn(x1[i], y1[i]) 

  }

  return(data.frame(cbind(x1,y1)))

}

output <- euler()

output

It outputs the right result, but with a warning message:

Warning message:
In y1[i] <- y1[i - 1] + step.size * (y.prime) :
number of items to replace is not a multiple of replacement length

Why am I getting this warning?


Solution

  • You are getting odd scoping issues because a function in the arguments of your function is called within the body .... and you think it's called on single values e.g. x1[i] but it looks like it's called on the growing vector x1. Just take the starting parameters out of the iterating function like so:

    y.prime.diff.eqn <- function(p, y) {return(5*y-p)}
    
    
    initial.condition.x1 <- 0
    initial.condition.y1 <- 0
    y.prime = y.prime.diff.eqn(initial.condition.x1, initial.condition.y1)
    
    # Define Euler Method's estimations
    euler <-  function(x1 = initial.condition.x1,
                       y1 = initial.condition.y1,
                       iter = 5,
                       step.size = 1) {
    
      for (i in 2:iter)
    
      {
    
        x1[i]     <-  x1[i-1] + step.size
        y1[i]     <-  y1[i-1] + step.size * (y.prime)
        y.prime   <-  y.prime.diff.eqn(x1[i], y1[i]) 
    
      }
    
      return(data.frame(cbind(x1,y1)))
    
    }
    
    output <- euler()
    
    output
    

    Same answer, no warning.