# Different Results for Double Pendulum

In this study, (Hidden Fractals in the Dynamics of the Compound Double Pendulum) the authors provide various fliptime fractals (of a double pendulum) for different length combinations. However, when I tried to replicate their results, the shapes (for $$3l_1 = l_2$$ in particular) did not match up. I used a fourth order Runge Kutta with adaptive step size and a max-step of 0.01 or 0.05 and simulated the pendulum for 100 seconds. Here are my results.

$$l_1 = l_2$$, $$2l_1 = l_2$$, and $$3l_1 = l_2$$ However, the article provides this image for $$3l_1 = l_2$$:

My $$2l_1 = l_2$$ and the article's $$3_l1 = l_2$$ resembles a square-ish shape. My $$2l_1$$ and $$3l_1$$ do not match the article, only $$l_1 = l_2$$ bares a similar resemblance. The article does not include a $$2l_1 = l_2$$ image.

Which image is right? I find it fishy that my $$2_l1$$ resembles their $$3_l1$$. Thanks :)

(some other parameters: masses are equal, gravity is 9.81, used solve_ivp() from scipy. I am happy to provide more if needed)

Here is my code:

from matplotlib.colors import LogNorm, Normalize
from scipy.integrate import solve_ivp, odeint
import matplotlib.pyplot as plt
from tqdm import tqdm
import seaborn as sns
import numpy as np

def derivs(t, y, length_1, length_2, mass_1, mass_2, gravity, angle_2_init):
angle_1, angle_1_d, angle_2, angle_2_d = y
sin1=np.sin(angle_1 - angle_2)
cos1=np.cos(angle_1 - angle_2)
cos2=(2 * mass_1 + mass_2 - mass_2 * np.cos(2 * angle_1 - 2 * angle_2))
angle_1_dd = (-gravity * (2 * mass_1 + mass_2) * np.sin(angle_1) - mass_2 * gravity * np.sin(angle_1 - 2 * angle_2) - 2 * sin1 * mass_2 * (angle_2_d ** 2 * length_2 + angle_1_d ** 2 * length_1 * cos1)) / (length_1 * cos2)
angle_2_dd = (2 * sin1 * (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 * cos1)) / (length_2 * cos2)
return [angle_1_d, angle_1_dd, angle_2_d, angle_2_dd]

def event(t, y, length_1, length_2, mass_1, mass_2, gravity, angle_2_init):
angle2 = y[2]
angle_2_init = np.deg2rad(angle_2_init)
return abs(angle2 - angle_2_init) - 2*np.pi
event.terminal=True
event.direction=1

def fliptime(length_1, length_2, mass_1, mass_2, angle_1_init, angle_2_init, angle_1_d_init, angle_2_d_init, gravity, time):
event.terminal=True
time_span = (0, time)
y0 = [np.deg2rad(angle_1_init), np.deg2rad(angle_1_d_init), np.deg2rad(angle_2_init), np.deg2rad(angle_2_d_init)]
sol = solve_ivp(derivs, time_span, y0, args=(length_1, length_2, mass_1, mass_2, gravity, angle_2_init), max_step=0.05, events=event)
if sol.status == 0:
return 100
else:
return sol.t_events[0][0]

angle_1_range = np.arange(-180, 180, 1)
angle_2_range = np.arange(-180, 180, 1)
fliptime_matrix = np.zeros((len(angle_1_range), len(angle_2_range)))
for i, angle_1 in tqdm(enumerate(angle_1_range), desc='angle_1', leave=False):
for j, angle_2 in tqdm(enumerate(angle_2_range), desc="angle_2", leave=False):
time = fliptime(1, 1, 1, 1, angle_1, angle_2, 0, 0, 9.81, 10000)
fliptime_matrix[i, j] = time

sns.heatmap(fliptime_matrix, square=True, cbar_kws={'label': 'Time to Flip'}, norm=LogNorm())
plt.xlabel('Angle 2 (degrees)')
plt.ylabel('Angle 1 (degrees)')
plt.title('Fliptime Heatmap')
plt.gca().invert_yaxis()
plt.show()
$$$$
`
• The default RK45 method is a 7-stage order 5 method, with an embedded 4th order method that is only used for step size control. It was developed by Dormand and Prince. Commented Mar 30 at 21:17
• Why don't you ask the authors of the paper in question? Commented Apr 1 at 16:27

## 1 Answer

You should make sure your integration is precise enough. I suggest you use the DOP853 algorithm of solve_ivp, which of higher order, and set tight error tolerances, e.g. rtol=atol=1e-12. Then assess if that changes your results.

EDIT: note that the authors have a slide dedicated to the effect of integration accuracy, which seems relatively low... So maybe is it indeed a mistake in the legend?