# Double spherical pendulum simulation difficulties

I've been working on a simulation project of mine that I kind of need help with. So, as made obvious by the title, I am attempting to simulate the motion of a double spherical pendulum. I am writing my code in Java; Processing to be exact.

I was inspired by The Coding Train's video on this subject and it really caught my attention as something to challenge myself with. I have followed the basic "flow" of the program in the afformentioned video. Here is the basic idea of how his code runs and how I want my code to run.

• Step 1: define variables such as angles, lengths, masses, gravity etc.
• Step 2: Draw the "bobs" of the pendulum and the lines connecting the and change this value every time the angles are changed.
• Step 3: I solved for the second time derivative (acceleration) in the equations of motion, I ran those equations. (I also solved for the equations of motion in Mathmatica and I'm pretty sure they are the same).
• Step 4: Based on the logic of the Coding Train video, each of the angles will change by it's derivative and each of the derivatives will change by their derivatives (This is a bad explanation and I don't have the firmest grip on this type of thing since i'm only 16). But based on this logic, velocity changes by velocity += acceleration and angle += velocity.
• Step 5: Repeat step 2-4 every frame. (this can easily be done in Processing)

My issue is that the acceleration is constantly increasing and making the sim go CRAZY. Is it how I have implemented the equations of motion? is it how I change the values of the angles?

I'm trying not to solve using Euler's method because it seems kind of hard to write code for. If it is necessary I should be able to write one but I don't want to deal with is unless absolutely necessary.

If I am wrong about how I implemented the equations or changed the values of the angles, how can I make this sim work?

P.S. Re-reading this, I may come off as kind of lazy, by not wanting to use an Euler's equation solver or anything, it's just that I want to make the code as efficient as possible so that the computer doesn't have to do some cooky ODE solving for each frame. I'm sorry.

• Did you transfer your equations to first order differential equations?. If you use Euler integration method, you need very small integration step size to have stable results, it is better to use implicit Euler Methode for the numerical simulation, or RK4 for example
– Eli
Commented Jul 13, 2020 at 17:41
• In the paper the equations of motion are already first order. How would I go around solving using the Runge-Kutta method? Commented Jul 14, 2020 at 1:38
• @I will write you the implementation code for that.
– Eli
Commented Jul 14, 2020 at 5:20
• Oh thanks. That's great Commented Jul 14, 2020 at 7:38

Theory

We want to solve this first order vector differential equations numerically.

$$\vec{\dot{y}}(t)=\vec{f}(t,\vec{y}(t))$$

with the initial condition

$$\vec{y}(t=t_0)=\vec{y}_0$$

Implementation

Explicit Integration method

\begin{align*} &t_a=t_0 \\ &\vec{y}_a=\vec{y}_0\\\\ &\text{while t_a \le t_f do}\\ &\vec{y}_n=\vec{y}_a+\vec{F}(t_a,\vec{y}_a,h) \\ &t_a\mapsto t_a+h\\ &\vec{y}_a\mapsto \vec{y}_n\\ &\text{end do} \end{align*}

Implicit Integration method

\begin{align*} &t_a=t_0 \\ &\vec{y}_a=\vec{y}_0\\\\ &\text{while t_a\le t_f do}\\ &\vec{y}_n=\vec{y}_a+\left[I_n+{F'}\right]^{-1}\vec{F}(t_a,\vec{y}_a,h) \\ &t_a\mapsto t_a+h\\ &\vec{y}_a\mapsto \vec{y}_n\\ &\text{end do} \end{align*}

where:

• $$t_0\quad$$ Start simulation time
• $$t_f\quad$$ Final simulation time
• $$h\quad$$ Simulation step time
• $$t_a=(t_0\,,t_0+h\,,t_0+2h+\ldots)$$
• $$I_n$$ unitiy $$n\times n$$ matrix
• $$F'=\frac{\partial \vec{F}}{\partial \vec{y}}\bigg|_{t_a,\vec{y}_a,h}\quad n\times n\,\,$$matrix
• $$\vec{y}_a(t_a)$$ The solution

The function $$F()$$ is depends on the Integration method that you want to implement:

RK4 Method

\begin{align*} &\vec F=\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\\ &k_1=h\,\vec{f}(t_a,\vec y_a)\\ &k_2=h\,\vec f(t_a+h/2,\vec y_a+k_1/2)\\ &k_3=h\,\vec f(t_a+h/2,\vec y_a+k_2/2)\\ &k_4=h\,\vec f(t_a+h,\vec y_a+k_3) \end{align*}

Euler: \begin{align*} &\vec F=h\,f(t_a,\vec y_a) \end{align*}

Euler implicit \begin{align*} &\vec F=h\,f(t_a,\vec y_a)\\ &F'=h\,f'(t_a,\vec y_a)=h\,\frac{\partial \vec f}{\partial \vec y}\bigg|_{t_a,\vec y_a} \end{align*}

Remark :

for your simulation choose $$h=0.001$$[s]