Search code examples
rfunctionode

What did wrong with my second derivative nonhomogeneous function?


I was trying to write a function for the differential equation

y'' + y = bcos(omega *t)

Here is my code

   library(deSolve)
   yini <- c(y1 = 2, y2 = 0)
   nonvdp1 <- function(t, y, parms) {
      b <- parms['b']
      omega <- parms['omega']
      dy1dt <- y[2]
      dy2dt <- b * cos(omega * t) - y[1]
      list(c(dy1dt, dy2dt))
   }
   output <- as.data.frame(ode(y = yini, func = nonvdp1, times = (0, 30, 0.1), parms = c(2, 2)))

however, the solution came out not quite right

   head(nonvdp1snl)
   time y1 y2
   1    0  2  0
   2    1 NA NA
   3    2 NA NA
   4    3 NA NA
   5    4 NA NA
   6    5 NA NA

Base on the initial condition, and parameters I selected, the solution is

y = (8/3)cos(t) - (2/3)cos(2t)

What did I do wrong with my code?


Solution

  • Your code had several issues, a forgotten seq for the time steps and missing names for the parameters. Compare with the following:

    library(deSolve)
    yini <- c(y1 = 2, y2 = 0)
    nonvdp1 <- function(t, y, parms) {
      b     <- parms['b']
      omega <- parms['omega']
      dy1dt <- y[2]
      dy2dt <- b * cos(omega * t) - y[1]
      list(c(dy1dt, dy2dt))
    }
    
    output <- ode(y = yini, func = nonvdp1, times = seq(0, 30, 0.1), 
      parms = c(b = 2, omega = 2))
    
    output
    
    plot(output)
    

    If you remove as.data.frame()you can use the built-in plot function plot.deSolve.