Introduction

When delving into the simulation of physical systems, a common approach involves solving differential equations that describe the evolution of the system over time. In this blog post, we’ll take a closer look at a specific example: the motion of an object subject to gravity and a time-dependent force. We’ll represent the governing differential equation in matrix form, implement a numerical solution using the Runge-Kutta method in Python, and visualize the results using Matplotlib.

Understanding the Differential Equation

The differential equation that models the system is expressed in matrix form as:

\[\frac{d}{dt} \begin{bmatrix} h \\ v \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} h \\ v \end{bmatrix} + \begin{bmatrix} 0 \\ -g + F(t) \end{bmatrix}\]

Here:

\( h \) represents the height of the object.
\( v \) represents its velocity.
\( g \) is the gravitational constant (9.81 m/s²).
\( F(t) \) is a time-dependent force, introduced here as \( 100 \sin(\pi t) \).

Implementing the Runge-Kutta Method

The Python code provided implements the fourth-order Runge-Kutta method. This numerical technique involves iteratively updating the solution based on weighted averages of function evaluations at different points within each time step. Let’s break down the key components of the code:

The `f` Function

def f(t, w):
    g = 9.81  # m/s^2, gravitational constant
    
    # Introduce a time-dependent force
    force = 100 * np.sin(np.pi * t)
    
    b = np.array([[0.], [-g + force]])

    L = np.array([[0., 1.], [0., 0.]])
    return L @ w + b

This function defines the right-hand side of the differential equation. It calculates the derivative of the state vector \( \begin{bmatrix} h \ v \end{bmatrix} \) with respect to time.

The `runge_kutta` Function

def runge_kutta(f, a, b, N, alpha):
    h = (b - a) / N
    t = a
    w = alpha
    
    time_steps = [t]
    heights = [w[0, 0]]
    velocities = [w[1, 0]]

    for i in range(1, N+1):
        K1 = h * f(t, w)
        K2 = h * f(t + h/2, w + K1/2)
        K3 = h * f(t + h/2, w + K2/2)
        K4 = h * f(t + h, w + K3)

        w = w + (K1 + 2*K2 + 2*K3 + K4) / 6
        t = a + i * h

        time_steps.append(t)
        heights.append(w[0, 0])
        velocities.append(w[1, 0])

    return time_steps, heights, velocities

This function performs the actual numerical integration using the Runge-Kutta method. It takes the differential equation function f, the time interval [a, b], the number of steps N, and the initial condition alpha.

Initial Conditions and Time Parameters

h0 = 100.0  # initial height
v0 = 0.0    # initial velocity
alpha = np.asarray([[h0], [v0]])
a = 0.0
b = 8.0
N = 100  # 

These variables set up the initial conditions of the system and define the time interval over which the simulation will run.

Solving the Differential Equation and Plotting

# Solve the differential equation using Runge-Kutta
time_steps, heights, velocities = runge_kutta(f, a, b, N, alpha)

# Plotting
fig, ax1 = plt.subplots()

color = 'tab:red'
ax1.set_xlabel('Time (s)')
ax1.set_ylabel('Height (m)', color=color)
ax1.plot(time_steps, heights, color=color)
ax1.tick_params(axis='y', labelcolor=color)

ax2 = ax1.twinx()
color = 'tab:blue'
ax2.set_ylabel('Velocity (m/s)', color=color)
ax2.plot(time_steps, velocities, color=color)
ax2.tick_params(axis='y', labelcolor=color)

fig.tight_layout()
plt.show()

This part of the code calls the runge_kutta function to obtain the solution and then visualizes the results using Matplotlib.

Conclusion

In summary, this blog post has covered the matrix representation of a time-dependent differential equation, the implementation of the Runge-Kutta method in Python, and the visualization of the simulation results. Understanding and applying numerical methods to solve such equations are crucial skills in various scientific and engineering disciplines. Experiment with different parameters and initial conditions to observe how the system’s behavior changes in response to these modifications.