I think this post will wrap up the series on Vectors and Simple Mechanics and we’ll look at Simple Harmonic Motion (SHM) and compare the numerical solutions to SHM using the naive step function from the previous post – aka the ‘*Euler*‘ step and the more accurate ‘*Euler-Cromer*‘ method.

Here’s the Euler Step from last time.

1 2 3 4 5 |
step :: Accnf -> Float -> State -> State step f dt st@(t, r, v) = (t', r', v') where t' = t + dt r' = r ^+^ v ^* dt v' = v ^+^ f st ^* dt |

a very simple change to the above yields the ‘*Euler-Cromer*‘ step where the ‘*new velocity*‘ rather than the ‘*old*‘ is used to determine the ‘*new*‘ position.

1 2 3 4 5 |
ecStep :: Accnf -> Float -> State -> State ecStep f dt st@(t, r, v) = (t', r', v') where t' = t + dt r' = r ^+^ v' ^* dt v' = v ^+^ f st ^* dt |

These two functions have the same signature, Accnf -> Float -> State -> State which allows us to generalise a solution based on a step function:

1 2 |
solutionWithStep :: ( Accnf -> Float -> State -> State )-> Accnf -> Float -> State -> [State] solutionWithStep stp a dt = iterate (stp a dt) |

which simply iterates the supplied step function to produce a list of *State*. This will now allow us to calculate solutions for both ‘*Euler*‘ and ‘*Euler-Cromer*‘.

A detailed discussion of SHM can be found here but essentially the force is proportional to the displacement and is directed towards the origin. This translates quite easily to one of our acceleration functions as

1 2 3 |
-- Harmonic oscillator where the force is -kx hosc :: Scalar -> (Time, Displacement, Velocity) -> V.Vector hosc k (_, r , _) = (-1)*k *^ r |

And, if we are to plot velocity (or displacement) against time for a given list of *State* then we need to extract the velocity (or displacement) from the *States*. This is achieved by these two simple functions:

1 2 3 4 5 |
timeDistanceX :: State -> Scalar timeDistanceX (_, V (x, _, _), _) = x timeDistanceXS :: [State] -> Path timeDistanceXS sts = zip [1..] (map timeDistanceX sts) |

The *timeDistanceX *function simply pattern matches on the x-component of the velocity vector within the state and the *timeDistanceXS* function zips up values, while there still are values, from the result of mapping *timeDistanceX* over all *States*. (Isn’t Haskell satisfying?).

Now we create the function

1 2 3 4 |
oneHosc :: Color -> (Accnf -> Float -> State -> State) -> Picture oneHosc c stp = pathPlot c . timeVelocityXS . take 50000 $ solutionWithStep stp (hosc 2.0) 0.01 (0, p0, v0) where p0 = V (0, 0, 0) v0 = V (50, 0, 0) |

which plots a path, in a given color, of a harmonic oscillator using the supplied step function. And if we are to compare the results of step functions then we need to produce a combined picture of each step function as shown in the *spring *function.

1 2 3 4 |
spring :: IO () spring = drawPics [oneHosc red step, oneHosc blue ecStep, color black (line [ (0, 50), (50000, 50) ])] |

In other words we’ll plot the ‘Euler’ solution in red and the ‘Euler-Cromer’ output in blue along with a horizontal line showing the start state. Here are such results!

The first graph of the solutions shows the first couple of oscillations and already the ‘*Euler*‘ solution (in red) is diverging – as can be seen by the gap between it and the horizontal line.

After a few more oscillations the difference becomes more pronounced whilst the ‘*Euler-Cromer*‘ solution looks to be stable.

and finally we can see that the ‘*Euler*‘ solution diverges completely whilst the ‘*Euler-Cromer*‘ solution remains stable!

This behaviour – i.e. the ‘*Euler*‘ solution diverging and the ‘*Euler-Cromer*‘ solution remaining stable is to be expected and there are many references and discussions about this behaviour e.g.

All the code is on Github and includes functions to plot a damped harmonic oscillator for both ‘Euler’ and ‘Euler-Cromer’. I would have liked to write the ‘definitive’ numerical method – Runge-Kutta but time etc. but if anyone wants to contribute this then please fork from Github and let me know how it goes!

Thanks for reading!