# Accurate way to numerically solve a system of second-order ODEs

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

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?

• There is a whole class of methods, actually. Jun 10, 2017 at 15:32
• Isn't this system chaotic? I thought there's no hope for a numerical solution as it is extremely sensitive to initial conditions. Jun 11, 2017 at 5:04
• @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. Jun 11, 2017 at 14:57

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.