Friday, 21 June 2013


Now we're going to look at some "classical" numerical analysis. Newton's second law is a differential equation, \[ m \frac{d^2 x}{dt^2} = F, \] summarizing the fact that you need to push heavy things harder than light things, x is a coordinate of as many variables as necessary. For our purposes we are going to write this one simple ordinary differential equation as two even simpler ones. Defining velocity $v = \frac{dx}{dt}$ Newton's equations become: \[ \frac{d}{dt} \left( \begin{array}{c} x \\ v \end{array} \right) = \left( \begin{array}{c} v \\ F/m \end{array} \right), \] calling the vector on the left $\vec{y}(t)$ and the one on the right $\vec{f}(t, \vec{y}(t))$ this becomes: \[ \frac{d\vec{y}(t)}{dt} = \vec{f}(t, \vec{y}(t)). \]

In this form we can make a very simple numerical solution: $\frac{d\vec{y}(t)}{dt}$ means the slope of $\vec{y}(t)$ at $t$. Let $h$ be a very small number, then: \[ \frac{d\vec{y}(t)}{dt} \approx \frac{\vec{y}(t+h) - \vec{y}(t)}{h}. \] Substituting this into the second form of Newton's equation we get, \[ \vec{y}(t+h) = \vec{y}(t) + h \vec{f}(t, \vec{y}(t)) \] so given $\vec{y}$ at $t$ we can calculate it a short time later at $t + h$ from the "force" $\vec{f}$. Let's try it on the standard example in physics, the pendulum.

From the diagram above, the pendulum moves on a circle of length $l$, the pendulum bob has traveled a distance $l \theta$ along the circle when the force on it is $-mg \sin(\theta)$, tangent to the circle. Newton's equations give, \[ \begin{array}{c} m l \frac{d^2 \theta}{dt^2} = -m g \sin(\theta)\\ \frac{d^2 \theta}{dt^2} + \frac{g}{l} \sin(\theta) = 0 \end{array} \] Because we're solving this numerically we don't have to make the small angle approximation $\sin(\theta) \approx \theta$ so the pendulum does indeed move on a circle. Click the animation below and drag the bob to somewhere in the box and let go. The red pendulum is the result of evolving the equation above using the Euler method.

The Euler method is a bad method, its advantage is how easy it is to program. The standard compromise between improvement and complexity is the 4${}^{th}$ order Runge-Kutta method. Runge from the Runge-Lenz vector and Kutta from the Kutta-Joukowski airfoil. This method is basically equivalent to Simpson's Rule. Explicitly, \[ \begin{array}{c} y(t+h) = y(t) + \frac{h}{6} \left( k_1 + 2k_2 + 2k_3 + k_4\right)\\ k_1 = f(t, y(t))\\ k_2 = f(t + h/2, y(t) + k_1h/2)\\ k_3 = f(t + h/2, y(t) + k_2h/2)\\ k_4 = f(t + h, y(t) + k_3h)\\ \end{array} \] This approximates the slope using the values at the beginning, the end and two in the middle, giving greater weight to the values in the middle. The blue bob in the previous animation (behind the red one) evolves according to the Runge-Kutta algorithm.

Given a decent algorithm for solving differential equations we'll use it for something more interesting, the double pendulum. Click the mouse in the box, drag the red bob to its start position click again and drag the blue bob then release. I'll give a derivation of the double (and triple and n-tuple) pendulum equations of motion next week.

No comments:

Post a comment