## Numerical integration

In order to ‘solve’ a differential equation for continuous time using a method of numerical integration, one could code it like in the spreadsheet assignment below. For R and Matlab there are so-called solvers available, functions that will do the integration for you. For R look at the Examples in package deSolve.

### Euler’s method and more…

The result of applying a method of numerical integration is called a numerical solution of the differential equation. The analytical solution is the equation which will give you a value of $$Y$$ for any point in time, given an initial value $$Y_0$$. Systems which have an analytical solution can be used to test the accuracy of numerical solutions.

#### Analytical solution

Remember that the analytical solution for the logistic equation is:

$Y(t) = \frac{K}{1 + \left(\frac{K}{Y_{0} - 1}\right) * e^{-r*t} }$

If we want to know the growth level $$Y_t$$ at $$t=10$$, with $$Y_0=.0001$$, $$r=1.1$$ and $$K=4$$, we can just fill it in:

# Define a function for the solution
logSol <- function(Y0, r, K, t){K/(1+(K/Y0-1)*exp(-r*t))}

# Call the function
logSol(Y0=.0001, r=1.1, K=4, t=10)
## [1] 2.398008

We can pas a vector of timepoints to create the exact solution, the same we would get if we were to iterate the differential/difference equation.

# Plot from t=1 to t=100
plot(logSol(Y0=.0001, r=1.1, K=4, t=seq(1,20)), type = "b",
ylab = expression(Y[t]), xlab = "t")
# Plot t=10 in red
points(10,logSol(Y0=.0001, r=1.1, K=4, t=10), col="red", pch=16)

#### Numerical solution (discrete)

If we would iterate the differential equation …

$\frac{dY}{dt} = Y_t * (1 + r - r * \frac{Y_t}{K})$

… as if it were a difference equation, that is, not simulating continuous time.

logIter <-  function(Y0,r,K,t){
N <- length(t)
Y <- as.numeric(c(Y0, rep(NA,N-2)))
sapply(seq_along(Y), function(t){ Y[[t+1]] <<- Y[t] * (1 + r - r * Y[t] / K)})
}

# Plot from t=1 to t=100
plot(logIter(Y0=.0001, r=1.1, K=4, t=seq(1,20)), type = "b",
ylab = expression(Y[t]), xlab = "t")
# Plot t=10 in red
points(10,logSol(Y0=.0001, r=1.1, K=4, t=10), col="red", pch=16)