Double spherical pendulum simulation difficulties

Question

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.

Answer 1

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]