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