## 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)
```