ODE solver - Numerical method doesn't fit analytical method.

Hey i'm trying to solve some simple ODE numerically and I know the explicit solutions for this ode,
dX = -X*0.01, which have the following explicit solution:
exp(-int_0^t 0.01)

but when i run the following r-code I'm not even close to the explicit solution. Forexample at time=40
the numerical solution gives 0.5999999 and the explicit solution is0.67032.

Why is the numerical method so wrong?

options("scipen"=999, "digits"=20)
Kolmogorov <- function(times, Initial_Value, ...) {
  with(as.list(c(Initial_Value)), {
    dX <- -X*0.01
    list(c(dX))
  })
}

Initial_Value <- c(X <- 1)
times_Prem = seq(0, 40, by = 0.01)


SSH_Prem <- ode(y = Initial_Value, times = times_Prem, func = Kolmogorov, parms = NULL)
SSH_Prem = data.frame(SSH_Prem)
SSH_Prem[100*40+1,c(1,2)]

First, when your code relies on a particular package, please include the code which loads it in your question so no one has to search out what package you are using, in this case,

library(deSolve)

With that out of the way, you really need to read the ode() helpfile, this should be your first, second, and third steps always.

?ode

# Relevant section:
func	  either an R-function that computes the values of the derivatives in the ODE
          system (the model definition) at time t, or a character string giving the name of a 
          compiled function in a dynamically loaded shared library.

          If func is an R-function, it must be defined as: func <- function(t, y, parms,...).
          t is the current time point in the integration, y is the current estimate of the
          variables in the ODE system. If the initial values y has a names attribute, the
          names will be available inside func. parms is a vector or list of parameters; ...
          optional) are any other arguments passed to the function.

          The return value of func should be a list, whose first element is a vector
          containing the derivatives of y with respect to time, and whose next elements
          are global values that are required at each point in times. The derivatives must
          be specified in the same order as the state variables y.

You would see your Kolmogorov() function should be defined as,

Kolmogorov <- function(t, y, ...) {
  list(-y*0.01)
}

NOTE: The names of the function don't matter but you should stick with the convention mentioned in the help file.

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.