# Robot Dynamics: Numerical simulation in Python

30 Dec 2018Recall from last episode, we were studying the dynamics of a simple pendulum in 2D. We derived the equations of motion of the pendulum

And the position of the pendulum is fully described by the angle

Our goal now is to solve the ODE given by the equations of motion, so we can compute the position $r(t)$ of the pendulum at any point in time starting from a know arbitrary position $r_0$.

# Analytic solution

The equation of motion we’re dealing with is not linear because of the $\sin{\theta}$ component. In practice, you never solve this kind of system analytically unless the system is simple (like our pendulum). Even for the is simple problem, the general solution is very tedious to derive analytically.

To make the problem easier to solve analytically, let’s assume the oscillations of the pendulum will be small. Which translates mathematically as

Our ODE becomes linear

Intuitively, the pendulum will oscillate over time. One way to model oscillation mathematically is with a $\sin{}$ function. We can write a candidate solution $\tilde{\theta}$

where $\alpha$, $\beta$, and $\gamma$ are unknown parameters. Computing the first and second derivatives of the candidate $\tilde{\theta}$ we get

Plugging this back in our ODE, we can solve for unknown parameters. We’ll also need an initial condition constraint: $\theta(t=0) = \theta_0$

So we have solved our equation of motion and can express $\theta$ explicity as function of time $t$

That’s nice, and easy to use to simulate our pendulum. But it only holds for small angles ($\sin{\theta} \approx \theta$).

# Numerical solution

Another popular approach that scales better to solving complex differential equations is numerical integration. The basic idea behind it is that the future position $x(t + \Delta t)$ of a particle can be approximated as the sum of the current position $x(t)$ and a small displacement proportional to the velocity $v(t)$ and the time increment $\Delta t$.

This is called forward integration, and it’s one of the most basic ways we can solve an ODE by numerical integration. We can solve first order ODEs this way, but the equation of motion of our pendulum is a second order ODE. Now I’m going to attempt a magic trick to transform the second order ODE into first order: a change of variable.

Let’s write our state $s = \begin{bmatrix} \theta \ \dot{\theta} \end{bmatrix}$. The first order time derivative of $s$ is $\dot{s} = \begin{bmatrix} \dot{\theta} \ \ddot{\theta} \end{bmatrix}$. How convenient! Using our beloved equation of motion, we can rewrite our problem as first order ODE in $s$

where $f$ is again an non arbitrary function of the state $s$. Great success! Let’s simulate it!

This is the Jupyter notebook used for simulation. With the numerical approach, we replicated the expected pendulum oscillation also obtained analytically (although with small angle approximation).

That’s all folks! In the next episode, we’ll look at the Lagrangian formalism and how it makes it easier to derive the dynamics of a system from energies instead of forces.