I'm writing a program meant to simulate a multiple pendulum with an arbitrary number of bobs, as shown below: enter image description here

Denote the angle of each pendulum from the vertical as $\theta_0,...,\theta_{n-1}$, and let $\mathbf{\Theta} = [\theta_0,...,\theta_{n-1}]^T$. The Lagrangian equations of motion are then of the form $A\ddot{\mathbf{\Theta}}=b$, where $A$ is an $n\times n$ matrix depending on $\mathbf{\Theta}$, and $b$ is an $n$-dimensional vector depending on $\dot{\mathbf{\Theta}}$ and $\mathbf{\Theta}$. My program currently works by keeping track of the current values for $\mathbf{\Theta}$ and $\dot{\mathbf{\Theta}}$, using this information to calculate $A$ and $b$ and solving for $\ddot{\mathbf{\Theta}}$, and then changing the values of each $\theta_i$ to $\theta_i+\dot{\theta}_i\Delta t+\frac 1 2 \ddot{\theta}_i\Delta t^2$, and the value of each $\dot{\theta}_i$ to $\dot{\theta}_i+\ddot{\theta}_i\Delta t$, in that order.

This works somewhat well, but the system tends to gain energy over time, due to how inaccurate this method is. Is there a better way to go about approximating the solution to equations of this form?

  • 2
    $\begingroup$ There is a whole class of methods, actually. $\endgroup$
    – Evgeny
    Jun 10, 2017 at 15:32
  • $\begingroup$ Isn't this system chaotic? I thought there's no hope for a numerical solution as it is extremely sensitive to initial conditions. $\endgroup$
    – user14717
    Jun 11, 2017 at 5:04
  • $\begingroup$ @user14717, You're right, I'm not expecting to be able to accurately predict the state of the pendulum in 10 seconds, but my main problem is that my current method allows the system to gain a substantial amount of energy over time, which leads to odd behavior. $\endgroup$
    – florence
    Jun 11, 2017 at 14:57

1 Answer 1


Let your system of equations be $$A(\Theta)\ddot{\Theta}=\vec{b}(\dot{\Theta},\Theta)$$ You can decompose it into first order ODE: $$\left\{\begin{array}{ll}A(\Theta)\dot{\Phi}=\vec{b}(\Phi,\Theta)\\ \Phi=\dot{\Theta} \end{array}\right.$$ And apply explicit methods (constrained time step through Courant number) such as Runge-Kutta of any order or implicit methods (unconstrained time step) like linear multistep methods. Try to use one that fits into your error criterium.

Since you system is nonlinear apply explicit methods if you want to implement them with ease, they will be direct to solve. If you want a stable and robust code, use implicit methods, but you must solve a nonlinear equation with some iterative scheme (Newton Raphson or fixed point will be adequate).


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .