Is this a viable way to calculate the Lyapunov Exponent?

Question

Is this a viable way to compute the Lyapunov Exponent for a double pendulum?

Here is the code for the double pendulum (You don't have to look at this part, it just returns the angles and angular velocity over time. I just placed it here just in case. I used the euler method to numerically solve the ODE):

import matplotlib.pyplot as plt
import numpy as np


def run_simulation(length_1, length_2, mass_1, mass_2, angle_1, angle_2, angle_1_d, angle_2_d, gravity, dt, num_steps):

    angle_1 = np.deg2rad(angle_1)
    angle_2 = np.deg2rad(angle_2)

    # Angle lists
    angle_1_lst = [angle_1]
    angle_2_lst = [angle_2]

    # Angular velocity lists
    angle_1_d_lst = [angle_1_d]
    angle_2_d_lst = [angle_2_d]

    # Angular acceleration lists
    angle_1_dd_lst = [0]
    angle_2_dd_lst = [0]

    # Trajectory lists (second pendulum)
    p1_x_trace = []
    p1_y_trace = []

    p2_x_trace = []
    p2_y_trace = []

    for _ in range(num_steps):
        angle_1_dd = (-gravity * (2 * mass_1 + mass_2) * np.sin(angle_1) - mass_2 * gravity * np.sin(angle_1 - 2 * angle_2) - 2 * np.sin(angle_1 - angle_2) * mass_2 * (angle_2_d ** 2 * length_2 + angle_1_d ** 2 * length_1 * np.cos(angle_1 - angle_2))) / (length_1 * (2 * mass_1 + mass_2 - mass_2 * np.cos(2 * angle_1 - 2 * angle_2)))
        angle_2_dd = (2 * np.sin(angle_1 - angle_2) * (angle_1_d ** 2 * length_1 * (mass_1 + mass_2) + gravity * (mass_1 + mass_2) * np.cos(angle_1) + angle_2_d ** 2 * length_2 * mass_2 * np.cos(angle_1 - angle_2))) / (length_2 * (2 * mass_1 + mass_2 - mass_2 * np.cos(2 * angle_1 - 2 * angle_2)))

        angle_1_dd_lst.append(angle_1_dd)
        angle_2_dd_lst.append(angle_2_dd)

        # Updating angular velocities
        angle_1_d += angle_1_dd * dt
        angle_2_d += angle_2_dd * dt

        angle_1_d_lst.append(angle_1_d)
        angle_2_d_lst.append(angle_2_d)

        # Updating angles
        angle_1 += angle_1_d * dt
        angle_2 += angle_2_d * dt

        angle_1_lst.append(angle_1)
        angle_2_lst.append(angle_2)

        # Trajectory
        x_1 = length_1 * np.sin(angle_1)
        x_2 = x_1 + length_2 * np.sin(angle_2)
        y_1 = - length_1 * np.cos(angle_1)
        y_2 = y_1 - length_2 * np.cos(angle_2)

        p1_x_trace.append(x_1)
        p1_y_trace.append(y_1)

        p2_x_trace.append(x_2)
        p2_y_trace.append(y_2)

    return angle_1_lst, angle_2_lst, angle_1_d_lst, angle_2_d_lst

And here is how I calculate the Lyapunov Exponent:

epsilon = 0.001 # Initial separation in angle
output_1 = run_simulation(1, 1, 1, 1, 30, 25, 0, 0, 9.81, 0.02, 50)
output_2 = run_simulation(1, 1, 1, 1, 30 + epsilon, 25, 0, 0, 9.81, 0.02, 50)

distance_lst = []

...

According to this study, the separation of the two trajectories is

distance = np.sqrt((output_2[0][index] - output_1[0][index])**2 + (output_2[1][index] - output_1[1][index])**2 + (output_2[2][index]*0.02 - output_1[2][index]*0.02)**2 + (output_2[3][index]*0.02 - output_1[3][index]*0.02)**2)

But this study states that the separation is the Euclidean distance of (theta_1, theta_2, theta_dot_1, theta_dot_2) of the angles and angular velocity, which does not include multiplying the velocities by dt (time step). For my code, I did not use * 0.02 for the angular velocities (which study is correct? or am I misinterpreting something?).

Then finally,

distance_array = np.array(distance_lst)
lyapunov_exponent = np.mean(np.log(distance_array[1:] / distance_array[:-1])) / 0.02
print("Lyapunov Exponent:", lyapunov_exponent)

Does this give the Lyapunov Exponent for the initial conditions? Thanks in advance. I was unsure which stack exchange I should post this on, so apologies if this post does not belong here. (I was directed here from stack overflow)

**Without code: **

I have the angles and angular velocity for two trajectories separated by a small initial condition over time. I found two studies (linked above) that used different methods of calculating the "separation" (Euclidean distance of the angles and angular acceleration, the other one did the same thing but multiplied the angular velocities by some time increment). I used the Euclidean distance (without multiplying the angular velocities by the time increment) for my distance calculation. Then, I calculate the mean of the natural log of the difference between two consecutive elements in my distance list and divide all that by the time increment (based off of $\lambda = \frac{1}{t}ln(\frac{\delta(t)}{\delta_0})$). Here is an example:

• distance_array = [1, 2, 3, 4, 5]
• distance_array[1:] / distance_array[:-1] = [2/1, 3/2, 4/3, 5/4]
• Then I take mean of the absolute value of the natural logarithm for every value in the list, and dividing by the time step I used for this experiment, 0.02 (or do I have to divide by the time of the entire simulation), should (hopefully) result in a (good enough) Lyapunov Exponent.