2. Solving differential equations: Numerical solution#
Differential equations can often be solved using a mathematical approach, which is called “analytical solution”. This requires good skills in math, as a minimum. Some differential equations are not even possible to solve analytically. Therefore it is very useful to be able to solve differential equations using an approximation and a computer. That is called a “numerical solution” or solving the equation numerically.
2.1. The forward Euler’s method: The simple but imprecise method.#
Let us consider a one-dimensional differential equation written in a general form:
where:
\(x\) is some variable, e.g. the firing rate of a neuron,
\(t\) is time
\(dx/dt\), is the rate of change of \(x\) and is equal to \(f()\) is some function of both \(x\) and \(t\).
To solve a differential equation, there must be an Initial Value that is used to calculate the change. This introduces an interesting property of differential equations: The values of the variables (in this case it is \(x\)) are different depending on the initial value, but the equation that governs the values is the same. Let us call the initial value of \(x\), i.e. the value at \(t=0\)
Although time and \(x\) are continuous variables, we will approximate their values at discrete times, \(t_n = n\Delta t\), where \(n\) indicates the \(n\)th point and \(n\Delta t\) is the associated time point.
Now we approximate the differential equation above as the rate of change in the time window of consideration, \(\Delta t\):
Hence, we can calculate the approximated updated value of \(x\) as
We can rewrite as
where the slope \(\frac{dx}{dt}\) is corresponding to the time point \(t_n\). These iterative equationas are approximation, which is more correct if \(f(x_n,t_n)\) is more linear. Also, choosing a smaller \(\Delta t\) increase the accuracy of the update, but this requires more updates for the same amount of time. The accuracy can be estimated to an order of:
Where \(T\) is the timescale of change of \(f\). So ideally \(\Delta t\) is small compared with the rate of change of \(T\). If not, then the approximation quickly breaks down. As we will discuss below, there is a much better estimation method called 4th order Runge-Kutta, which has an accuracy of \((\frac{\Delta t}{T})^4\), i.e. 10 000 time better than the forward Euler method.
2.2. Using for-loop to solve numerically#
Since the forward Euler method is an iterative process, i.e. calculating the updates in \(x\) as \(x_{n+1} = x_{n}+f(x_n,t_n)\Delta t\), it makes sense to use a for loop. A for loop in Python is a control structure used to iterate over a sequence (like a list, tuple, string, or range) and repeat a block of code for each item in that sequence. The basic syntax of a for loop is the following:
for item in sequence:
# Do something with item
Cell In[1], line 2
# Do something with item
^
SyntaxError: incomplete input
Here is an example of a for loop using a list, where the script is going through a loop printing each element in the list:
fruits = ["apple", "banana", "cherry"]
for fruit in fruits:
print(fruit)
Here it is using range() to repeat an action:
for i in range(5):
print(i)
Key Points:
The for loop automatically gets each item from the sequence.
It runs the indented block once for each item.
No need to manage index counters manually like in some other languages (though you can if needed).
2.3. Exercise: Integrate the membrane equation using the forward Euler’s method#
This is a classic leaky integrator model from neuroscience. Let’s solve it with:
\(V_{rest}=−65 mV\)
\(\tau = 10\) ms
Initial condition: \(V(0)=-55\) mV
Time step: \(\Delta t = 1\) ms
Simulate for 100 ms
2.4. Euler’s method formula:#
import numpy as np
import matplotlib.pyplot as plt
# Parameters
V_rest = -65 # Resting potential (mV)
tau = 10 # Time constant (ms)
dt = 1 # Time step (ms)
T = 100 # Total simulation time (ms)
# Time array
time = np.arange(0, T + dt, dt)
# Initialize voltage array
V = np.zeros(len(time))
V[0] = -55 # Initial condition
# Euler integration using a for loop
for i in range(1, len(time)):
dVdt = -(V[i-1] - V_rest) / tau
V[i] = V[i-1] + dt * dVdt
# Plot the result
plt.plot(time, V)
plt.xlabel('Time (ms)')
plt.ylabel('Membrane potential V (mV)')
plt.title("Euler's Method: Leaky Integrator")
plt.grid(True)
plt.show()
Exercise: Plot the solution to the following differential equations:
\( \frac{dx}{dt} = exp(t/\tau)\), \(x(t=0)=1\). \(\tau = 0.1\). Plot from t=-5 to 5. Try varying the \(\Delta t\), and other parameters
\( \frac{dx}{dt} = 2sin(2x)\), \(x(t=0)=0\). Plot from t=-5 to 5. Try varying the \(\Delta t\), and other parameters
\( \frac{dx}{dt} = 2x^2\), Plot from t=-2 to 2. Try varying the \(\Delta t\), and other parameters
2.4.1. Euler’s method of a 2-dimensional system#
Consider a system of 2-dimensional coupled linear differential equations:
To integrate this using Euler’s method we do the following:
Here is the python code for numerically integrating this system:
import numpy as np
import matplotlib.pyplot as plt
# Define the system parameters
a, b, c, d = -2, 1, 1, -1
# Initial conditions
x0, y0 = 1, 0
# Time parameters
t_start, t_end, dt = 0, 10, 0.1
time_steps = int((t_end - t_start) / dt)
# Arrays to store the solutions
x_values = np.zeros(time_steps + 1)
y_values = np.zeros(time_steps + 1)
t_values = np.linspace(t_start, t_end, time_steps + 1)
# Set initial conditions
x_values[0], y_values[0] = x0, y0
# Euler's method iteration
for n in range(time_steps):
x_values[n + 1] = x_values[n] + dt * (-2 * x_values[n] + y_values[n])
y_values[n + 1] = y_values[n] + dt * (x_values[n] - y_values[n])
# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(t_values, x_values, label='x(t)')
plt.plot(t_values, y_values, label='y(t)')
plt.title("Solution of the 2D Linear System using Euler's Method")
plt.xlabel('Time t')
plt.ylabel('Values')
plt.legend()
plt.grid(True)
plt.show()
2.5. A better method: 4th order Runge-Kutta#
A very popular numerical method for integrating ordinary differential equations is the 4th order Runge-Kutta (RK4). It improves on Euler’s method by taking multiple estimates of the slope at each step and combining them in a weighted average. Let us again look at the differential equation from above:
The RK4 computes the values of \(x_{n+1}\) from the current value of \(x_n\) using the following set of estimations:
Then:
Where
\(k_1\) is the slope at the beginning of the interval.
\(k_2\) and \(k_3\) are slopes at the midpoint (with different estimates).
\(k_4\) is the slope at the end of the interval.
The weighted average gives better accuracy than a single slope.
In fact, the accuracy is
which is 10.000 times better than when using Forward Euler method.
2.6. Lorenz system: example of a 3-dimenstional system#
The Lorenz system is a famous set of three coupled, nonlinear differential equations. Using the 4th Order Runge-Kutta (RK4) method to solve it numerically is a classic example in nonlinear dynamics and chaos theory.
The Lorenz equations are:
Where typical parameters are:
\(\sigma = 10\)
\(\rho =28\)
\(\beta = \frac{8}{3}\)
Goal is to numerically integrate these equations using the RK4. Here is a suggestion of how to do it:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
# Parameters
sigma = 10
rho = 28
beta = 8 / 3
dt = 0.01
T = 50
num_steps = int(T / dt)
# Arrays to hold the results
x = np.zeros(num_steps)
y = np.zeros(num_steps)
z = np.zeros(num_steps)
# Initial conditions
x[0], y[0], z[0] = 1, 1, 1
# Lorenz system
def lorenz(x, y, z):
dx = sigma * (y - x)
dy = x * (rho - z) - y
dz = x * y - beta * z
return dx, dy, dz
# RK4 Integration
for i in range(num_steps - 1):
dx1, dy1, dz1 = lorenz(x[i], y[i], z[i])
dx2, dy2, dz2 = lorenz(
x[i] + 0.5 * dt * dx1,
y[i] + 0.5 * dt * dy1,
z[i] + 0.5 * dt * dz1
)
dx3, dy3, dz3 = lorenz(
x[i] + 0.5 * dt * dx2,
y[i] + 0.5 * dt * dy2,
z[i] + 0.5 * dt * dz2
)
dx4, dy4, dz4 = lorenz(
x[i] + dt * dx3,
y[i] + dt * dy3,
z[i] + dt * dz3
)
x[i+1] = x[i] + (dt / 6) * (dx1 + 2*dx2 + 2*dx3 + dx4)
y[i+1] = y[i] + (dt / 6) * (dy1 + 2*dy2 + 2*dy3 + dy4)
z[i+1] = z[i] + (dt / 6) * (dz1 + 2*dz2 + 2*dz3 + dz4)
# Plotting the Lorenz attractor
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, lw=0.5)
ax.set_title("Lorenz Attractor (RK4)")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
plt.tight_layout()
plt.show()
To animate the 3D properties of the Lorenz attractor:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
# Define the Lorenz system
def lorenz(state, t, sigma, rho, beta):
x, y, z = state
dx = sigma * (y - x)
dy = x * (rho - z) - y
dz = x * y - beta * z
return [dx, dy, dz]
# Parameters
sigma = 10
rho = 28
beta = 8 / 3
dt = 0.01
T = 40
t = np.arange(0, T, dt)
# Initial condition
initial_state = [1.0, 1.0, 1.0]
# Integrate the system
solution = odeint(lorenz, initial_state, t, args=(sigma, rho, beta))
x, y, z = solution.T # Transpose to get individual arrays
# Set up 3D plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
line, = ax.plot([], [], [], lw=0.8)
# Set axis limits
ax.set_xlim(np.min(x), np.max(x))
ax.set_ylim(np.min(y), np.max(y))
ax.set_zlim(np.min(z), np.max(z))
ax.set_title("Lorenz Attractor (Animated with odeint)")
# Animation update function
def update(num):
line.set_data(x[:num], y[:num])
line.set_3d_properties(z[:num])
return line,
# Create animation (no blitting for 3D!)
ani = FuncAnimation(fig, update, frames=len(t), interval=1, blit=False)
plt.style.use('dark_background')
plt.show()
If you want to save the the plot as a animation mp4 file use the following:
ani.save('lorenz_animation.mp4', fps=30, dpi=150)
# Use plt.style.use('dark_background') before plotting to make it visually striking.