Forward Euler method from scratch

A differential equation gives a rate of change, not all future values directly. The forward Euler method uses the current rate to estimate the next value, then repeats the calculation to construct an approximate trajectory.

Biological scenario

Suppose 100 people are infectious and no new infections occur. Infectious people recover at rate \(\gamma=0.2\) per day.

We want to estimate the remaining infectious population over 10 days using only the differential equation and simple Python calculations.

The recovery differential equation

The recovery-only model is

\[\frac{dI}{dt}=-\gamma I.\]

The minus sign means that the infectious population decreases. At \(I=100\),

\[\frac{dI}{dt}=-0.2\times100=-20\quad\text{people per day}.\]

This is the instantaneous rate at the current state. As \(I\) decreases, the recovery rate also becomes smaller.

From a derivative to an approximate change

Over a short interval \(\Delta t\), Euler’s method approximates

\[\text{change}\approx\text{current rate}\times\Delta t.\]

\[I_{n+1}=I_n+\left(-\gamma I_n\right)\Delta t.\]

The subscript \(n\) identifies the current stored time point. The subscript \(n+1\) identifies the next one.

The four Euler operations

1 CURRENT VALUERead \(I_n\).

2 CURRENT RATECalculate \(-\gamma I_n\).

3 APPROXIMATE CHANGEMultiply the rate by \(\Delta t\).

4 NEXT VALUEAdd the change to \(I_n\).

Calculate the first two steps by hand

Use \(I_0=100\) and \(\Delta t=1\) day.

Current\(I_0=100\)

Rate\(-0.2\times100=-20\)

Change\(-20\times1=-20\)

Next\(I_1=100-20=80\)

Current\(I_1=80\)

Rate\(-0.2\times80=-16\)

Change\(-16\times1=-16\)

Next\(I_2=80-16=64\)

Translate one step into Python

rate = -gamma * I[n]
change = rate * dt
I[n + 1] = I[n] + change

Python expression	Meaning
`I[n]`	The current infectious value
`-gamma * I[n]`	The rate evaluated at the current value
`rate * dt`	The approximate change during one step
`I[n + 1]`	The storage position for the next value

Repeat the step with a loop

for n in range(number_of_steps):
    rate = -gamma * I[n]
    I[n + 1] = I[n] + rate * dt

for begins a repetition.
range(number_of_steps) produces the indices \(0,1,2,\ldots\) needed for the updates.
The indented statements run once for each index.
The next value must be calculated from the current value, I[n].

Run the complete Euler program

The program also calculates the exact solution \(I(t)=I_0e^{-\gamma t}\) so that the numerical approximation can be checked.

Interactive PythonEuler recovery model

import numpy as np
import matplotlib.pyplot as plt

# Biological and numerical settings
I0 = 100
gamma = 0.2
start_time = 0
end_time = 10
dt = 1

# Check the numerical settings
if dt <= 0:
    raise ValueError("dt must be positive.")
if gamma * dt > 1:
    raise ValueError("Choose dt so that gamma * dt is at most 1.")

# Prepare the time grid and storage
number_of_steps = int(round((end_time - start_time) / dt))
time = np.linspace(start_time, end_time, number_of_steps + 1)
I_euler = np.zeros(number_of_steps + 1)
I_euler[0] = I0

# Forward Euler method
for n in range(number_of_steps):
    rate = -gamma * I_euler[n]
    change = rate * dt
    I_euler[n + 1] = I_euler[n] + change

# Exact solution for comparison
I_exact = I0 * np.exp(-gamma * time)
absolute_error = np.abs(I_exact - I_euler)

# Display a results table
print("Day   Euler     Exact     Absolute error")
print("-----------------------------------------")
for n in range(len(time)):
    print(
        f"{time[n]:3.0f}  {I_euler[n]:8.3f}  "
        f"{I_exact[n]:8.3f}  {absolute_error[n]:14.3f}"
    )

print("-----------------------------------------")
print(f"Maximum error: {absolute_error.max():.3f}")

# Plot the outcome
plt.figure(figsize=(7, 4.5))
plt.plot(time, I_euler, "o-", label="Forward Euler")
plt.plot(time, I_exact, "--", label="Exact solution")
plt.xlabel("Time (days)")
plt.ylabel("Infectious population")
plt.title("Recovery-only model")
plt.grid(alpha=0.25)
plt.legend()
plt.tight_layout()

Output

Run the code to see the result.

Expected early values

Day	Euler value	Exact value	Absolute error
0	100.000	100.000	0.000
1	80.000	81.873	1.873
2	64.000	67.032	3.032
3	51.200	54.881	3.681

Biological and numerical interpretation

The infectious population declines because the model contains recovery but no new infection. Each day, Euler’s method removes 20% of the current infectious value. The amount removed becomes smaller as the infectious population declines.

The Euler values differ from the exact curve because Euler assumes the current rate remains fixed across each whole step. A smaller time step updates the rate more frequently and normally improves the approximation.

Important conditions and common mistakes

\(\Delta t>0\) is required.
The time grid and storage array must have compatible lengths.
The loop stops at the final update; otherwise I[n + 1] would exceed the array.
Use the current value on the right-hand side before storing the next value.
For this decay equation, \(0\leq\gamma\Delta t\leq1\) prevents a single Euler step from producing a negative population.
A time step can satisfy non-negativity and still be too large for good accuracy.

Try a smaller time step

Change dt = 1 to dt = 0.25.

Predict what happens to the number of Euler steps.
Run the program.
Compare the maximum error with the original result.
Inspect how closely the Euler graph follows the exact curve.