This page gives some simple code for implementing Euler's method in Python and R for two particular examples. See here for a brief introduction to Euler's method.

Python code to implement Euler's method for a simple example

Here's some simple code that implements Euler's method for the ordinary differential equation $f(x)=f^{\prime }(x)$ with $f(0)=1$ in Python:

# Python-code

import numpy as np

# x0, initial value of x
# xn, final value of x
# y0, initial value of y
# h, step size
# f, function

def euler(x0, xn, y0, h, f):
    x = np.arange(x0, xn, h)
    y = np.zeros(x.shape)
    y[0] = y0
    for i in range(1, x.shape[0]):
         y[i] = y[i-1] + h * f(y[i-1])

   return {'x':x, 'y':y}

def dxdy(x):   # function for f'(x) = f(x)
    return x      # dx/dy = x

solution = euler(x0=0, xn=2.5, h=0.5, f=dxdy, y0=1)

print(solution['y'])

R code to implement Euler's method for a simple example

Here's some simple code that implements Euler's method for the ordinary differential  equation $f(x)=f^{\prime }(x)$ with $f(0)=1$ in R:

# R-code

# x0, initial value of x
# xn, final value of x
# y0, initial value of y
# h, step size
# f, function
euler = function(x0, xn, y0, h, f){
    x = seq(x0, xn, h)
    y = array(data=0, dim=length(x))
    y[1] = y0
    for (i in 2:length(x)){
       y[i] = y[i-1] + h * f(y[i-1])
    }
    return(list(x=x, y=y))
}

dxdy = function(x)   # function for f'(x) = f(x)
    return(x)               # dx/dy = x

solution = euler(x0=0, xn=2, y0=1, h=0.5, f=dxdy)

print(solution$y)

Python code to implement Euler's method for the SIR model

The SIR model is another interesting example involving ordinary differential equations. It's one of the simplest models to describe the spread of an infectious disease — you can find out more about the model in this article.

The equations defining the model are:

$$dS/dt=-\beta S(t)I(t)$$

$$dI/dt=\beta S(t)I(t)-\gamma I(t)$$

$$dR/dt=\gamma I(t)$$

where $S$ denotes the proportion of people in a population that are susceptible to the disease, $I$ denotes the proportion of infected people, and $R$ denotes the proportion of people who have recovered and are immune.  The parameter $\beta$ is the rate at which people pass from the susceptible to the infected class and  $\gamma$ is the rate at which people pass from the infected to the recovered class. Here is some simple code for implementing Euler's method for this model in Python:

# Python-code

import numpy as np
import matplotlib.pyplot as plt

# t0, initial value of time
# tn, final value of time
# y0, initial values
# h, step size
# f, ODE model equations
def euler(t0, tn, y0, h, f):                                  # vectorised version of Euler's method
    t = np.arange(t0, tn, h)
    y = np.zeros((t.shape[0], y0.shape[0]))
    y[0] = y0
    for i in range(1, t.shape[0]):
       y[i] = y[i-1] + h * f(y[i-1])

   return {'t':t, 'y':y}

def sir_model(x):                                           # SIR model ODE equations
    beta = 0.6                                                   # rate parameter S -> I
    gamma = 0.2                                              # rate parameter I -> R
    S = - beta * x[0] * x[1]                              # dS/dt = -beta * S(t) * I(t)
    I =   beta * x[0] * x[1] - gamma * x[1]     # dI/dt = beta * S(t) * I(t) - gamma * I(t).                                                   R = - S - I                                                   # dR/dt = gamma * I(t)
   return np.array([S, I, R])

y0 = np.array([0.9, 0.1, 0]) # initial values of S, I and R as a proportion of the total population
solution = euler(t0=0, tn=40, y0=y0, h=0.1, f=sir_model)

# plot
labels = ['S(t), susecptible', 'I(t), infectious', 'R(t), recovered']
plt.plot(solution['t'], solution['y'], label=labels)
plt.legend(loc='center right')
plt.show()
 

R code to implement Euler's method for the SIR model

Here is some simple code for implementing Euler's method for the SIR model in R:

# R-code

# t0, initial value of time
# tn, final value of time
# y0, initial values
# h, step size
# f, ODE model equations
euler = function(t0, tn, y0, h, f){                   # vectorised version of Euler's method
    t = seq(t0, tn, h)
    y = array(data=0, dim=c(length(t),length(y0)))
    y[1,] = y0
    for (i in 2:length(t)){
      y[i,] = y[i-1,] + h * f(y[i-1,])
    }
    return(list(t=t, y=y))
}

sir_model = function(x){                            # SIR model ODE equations
   beta = 0.6                                                 # rate parameter S -> I
   gamma = 0.2                                            # rate parameter I -> R
   S = - beta * x[1] * x[2]                            # dS/dt = -beta * S(t) * I(t)
   I =   beta * x[1] * x[2] - gamma * x[2]   # dI/dt =  beta * S(t) * I(t) - gamma * I(t)
   R = - S - I                                                 # dR/dt = gamma * I(t)
   return(c(S, I, R))
}

y0 = c(0.9, 0.1, 0) # initial values of S, I and R as a proportion of the total population
solution = euler(t0=0, tn=40, y0=y0, h=0.1, f=sir_model)

# plot
labels = c('S(t), susecptible', 'I(t), infectious', 'R(t), recovered')
colours = c('blue', 'red', 'green')
linetypes = c(1,2,3)
matplot(solution$t, solution$y, type='l', lty=linetypes, xlab='time', ylab='proportion', col=colours)
legend(14,0.6, legend=labels, cex=0.7, col=colours, lty=linetypes)