Code:
import numpy as np
import matplotlib.pyplot as plt
def harmonic_oscillator(t, y):
return np.array([y[1], -k * y[0]])
def pendulum(t, y):
return np.array([y[1], -g / L * np.sin(y[0])])
k = 1.0
g = 9.81
L = 1.0
initial_conditions = np.array([np.pi / 4, 0.0])
t_start = 0.0
t_end = 10.0
dt = 0.01
num_steps = int((t_end - t_start) / dt)
rk4_results = np.zeros((num_steps + 1, 2))
rk2_results = np.zeros((num_steps + 1, 2))
rk4_results[0] = initial_conditions
rk2_results[0] = initial_conditions
for i in range(num_steps):
t = t_start + i * dt
y = rk4_results[i]
k1 = dt * harmonic_oscillator(t, y)
k2 = dt * harmonic_oscillator(t + dt/2, y + k1/2)
k3 = dt * harmonic_oscillator(t + dt/2, y + k2/2)
k4 = dt * harmonic_oscillator(t + dt, y + k3)
rk4_results[i + 1] = y + (k1 + 2*k2 + 2*k3 + k4) / 6
for i in range(num_steps):
t = t_start + i * dt
y = rk2_results[i]
k1 = dt * harmonic_oscillator(t, y)
k2 = dt * harmonic_oscillator(t + dt, y + k1)
rk2_results[i + 1] = y + (k1 + k2) / 2
analytical_solution = initial_conditions[0] * np.cos(np.sqrt(k) * np.arange(t_start, t_end + dt, dt))
rk4_errors = np.abs(analytical_solution - rk4_results[:, 0])
rk2_errors = np.abs(analytical_solution - rk2_results[:, 0])
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(np.arange(t_start, t_end + dt, dt), rk4_results[:, 0], label='RK4')
plt.plot(np.arange(t_start, t_end + dt, dt), rk2_results[:, 0], label='RK2')
plt.plot(np.arange(t_start, t_end + dt, dt), analytical_solution, label='Analytical')
plt.xlabel('Time')
plt.ylabel('Angle')
plt.legend()
plt.title('Harmonic Oscillator')
plt.subplot(2, 1, 2)
plt.plot(np.arange(t_start, t_end + dt, dt), rk4_errors, label='RK4 Error')
plt.plot(np.arange(t_start, t_end + dt, dt), rk2_errors, label='RK2 Error')
plt.xlabel('Time')
plt.ylabel('Error')
plt.legend()
plt.title('Error Comparison')
plt.tight_layout()
plt.show()
Step-by-step explanation:
Certainly! Let's break down the code and the steps involved in solving the harmonic oscillator using the Runge-Kutta (RK) methods, comparing RK4 and RK2, and analyzing the errors.
1. Import Required Libraries:
We begin by importing the necessary libraries, `numpy` for numerical operations and `matplotlib` for plotting the results.
2. Define Differential Equations:
We define two functions: `harmonic_oscillator` and `pendulum`. These functions represent the differential equations for the harmonic oscillator and pendulum, respectively. They take the current time `t` and the current state vector `y` (angle and angular velocity) as inputs and return the derivatives of the state variables.
3. Set Parameters and Initial Conditions:
We set the values of the physical parameters such as the spring constant `k`, acceleration due to gravity `g`, length of the pendulum `L`, and initial conditions `initial_conditions`.
4. Define Time Parameters:
Set the starting time `t_start`, ending time `t_end`, time step `dt`, and calculate the number of time steps `num_steps` based on the given time range and step size.
5. Initialize Result Arrays:
Create arrays `rk4_results` and `rk2_results` to store the results of the RK4 and RK2 methods. These arrays will hold the angles and angular velocities at each time step.
6. Implement RK4 Method:
Iterate through each time step using a loop. For each step, calculate the intermediate values `k1`, `k2`, `k3`, and `k4` using the RK4 formula for the harmonic oscillator's differential equation. Then, use these intermediate values to update the state vector using a weighted average of `k` values. Store the updated state in the `rk4_results` array.
7. Implement RK2 Method:
Similar to RK4, use a loop to iterate through each time step. For each step, calculate the intermediate values `k1` and `k2` using the RK2 formula for the harmonic oscillator's differential equation. Update the state vector using the weighted average of `k` values and store the results in the `rk2_results` array.
8. Calculate Analytical Solution:
Calculate the analytical solution for the harmonic oscillator using the initial conditions. This solution involves a cosine function based on the initial angle and the square root of the spring constant.
9. Calculate Errors and Plot Results:
Calculate the absolute errors between the analytical solution and the numerical solutions obtained using RK4 and RK2 methods. Then, create a plot to visualize the angles obtained from both methods and the analytical solution. Additionally, create a subplot to visualize the errors of both methods.
10. Display Plots:
Finally, display the generated plots using `matplotlib`.
By comparing the results obtained from the RK4 and RK2 methods to the analytical solution, you can observe how well these numerical methods approximate the actual behavior of the harmonic oscillator. The error comparison plot helps in understanding which method provides more accurate results over the given time range and step size.