I am getting this error, and I’m having trouble finding a way to resize the array to run. is there away to fix this or do I need to have a different approach. (calculate_error
return np.linalg.norm(X_values_fine – X_values_coarse[:len(X_values_fine)]) / np.linalg.norm(X_values_fine)

ValueError: operands could not be broadcast together with shapes (40,6) (20,6))

import matplotlib.pyplot as plt

# Constants
m1 = 12000  # kg
m2 = 10000  # kg
m3 = 8000   # kg
k1 = 1800 * 1000  # N/m
k2 = 2400 * 1000  # N/m
k3 = 3000 * 1000  # N/m
t_end = 20  # seconds

# Time step sizes
h_values = [1, 0.5, 0.25]

def F(X):
    x1, x2, x3, v1, v2, v3 = X
    dx1 = v1
    dx2 = v2
    dx3 = v3
    dv1 = (-k1*x1 + k2*(x2 - x1)) / m1
    dv2 = (k2*(x1 - x2) + k3*(x3 - x2)) / m2
    dv3 = (k3*(x2 - x3)) / m3
    return np.array([dx1, dx2, dx3, dv1, dv2, dv3])

def midpoint_method(h):
    # Initial conditions
    X = np.array([0, 0, 0, 1, 0, 0])  # [x1, x2, x3, v1, v2, v3]
    t = 0
    num_steps = int(t_end / h)
    
    # Storage for plotting
    t_values = []
    X_values = []
    
    for _ in range(num_steps):
        t_values.append(t)
        X_values.append(X)
        
        # Midpoint method steps
        X_half = X + 0.5 * h * F(X)
        X_next = X + h * F(X_half)
        
        # Update
        X = X_next
        t += h
    
    return np.array(t_values), np.array(X_values)

# Solve for different step sizes
results = {}
for h in h_values:
    t_values, X_values = midpoint_method(h)
    results[h] = (t_values, X_values)

# Plot results
fig, axes = plt.subplots(3, 2, figsize=(15, 12))
titles = ['Displacement x1', 'Displacement x2', 'Displacement x3', 'Velocity v1', 'Velocity v2', 'Velocity v3']
for i, (title, idx) in enumerate(zip(titles, [0, 1, 2, 3, 4, 5])):
    ax = axes[i//2, i%2]
    for h in h_values:
        t_values, X_values = results[h]
        ax.plot(t_values, X_values[:, idx], label=f'h={h}')
    ax.set_title(title)
    ax.set_xlabel('Time (s)')
    ax.set_ylabel('Value')
    ax.legend()

plt.tight_layout()
plt.show()

# Analysis of convergence
def calculate_error(X_values_fine, X_values_coarse):
    return np.linalg.norm(X_values_fine - X_values_coarse[:len(X_values_fine)]) / np.linalg.norm(X_values_fine)

errors = {}
for h in h_values[:-1]:
    _, X_values_fine = results[h/2]
    _, X_values_coarse = results[h]
    error = calculate_error(X_values_fine, X_values_coarse)
    errors[h] = error

for h, error in errors.items():
    print(f'Error between h={h/2} and h={h}: {error:.2e}')`