1
$\begingroup$

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 3l_1 = l_2

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

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()
```
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Mar 30 at 21:17
  • $\begingroup$ Why don't you ask the authors of the paper in question? $\endgroup$ Apr 1 at 16:27

1 Answer 1

2
$\begingroup$

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?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.