This is an attempt to simulate a double spring model using 4th order Runge-Kutta simulation in Haskell. The model consists of two springs where the first is fixed at one location. At the bottom of the first spring a point mass is located. Attached to the first mass is a second spring which ends in a second point mass. At first my model seemed to be unstable when the mass was above the spring. Finally I realized that this was due to an erroneous calculation of the angle (`phi`) from the location of the masses (`x1, y1, x2, y2`). `atan` has two solutions (separated by `pi`) where only one is valid. I had to add pi (via `addPi` function) to handle that case.

The haskell code for Runge-Kutta is largely taken from Runge-Kutta with Haskell, but I have changed the algorithm to match the following algorithm instead, 16.1 Runge-Kutta Method. gloss is used to output graphics.

The haskell code is shown below.

``````module Main (main) where
import Graphics.Gloss
import System.Random

{- "(Double) Spring Pendulum"
See: http://www.myphysicslab.com/dbl_spring2d.html
Compile: ghc -O2 SpringPendulum.hs
Optimize: strip SpringPendulum -o SpringPendulum
Run: ./SpringPendulum -}

instance (Num a) => Num [a] where
(+) = zipWith (+)       -- e.g. [a,b]+[c,d]=[a+c,b+d]
negate = fmap negate
fromInteger = undefined
(*) = undefined
abs = undefined
signum = undefined

-- One iteration of runge kutta. Calculate y(t+h) from y(t)
-- f'(y) is system derivate, h step length and y current state, y(t)
-- See: http://tnt.phys.uniroma1.it/twiki/pub/TNTgroup/AngeloVulpiani/runge.pdf
rk4 :: (Functor f, Floating a, Num (f a)) => (f a -> f a) -> a -> (f a) -> (f a)
rk4 f' h y =
let
(.*) n = fmap (*n)
shf = ((1/2).*)
k0 = h .* (f'   y)
k1 = h .* (f' \$ y + shf k0)
k2 = h .* (f' \$ y + shf k1)
k3 = h .* (f' \$ y + k2)
in
y + (1/6) .* (k0 + 2 .* k1 + 2 .* k2 + k3)

-- Integrate system using runge kutta approximation
integrate :: (Functor f, Floating a, Num (f a)) => (f a -> f a) -> f a -> a -> [f a]
integrate f' y0 h = iterate (rk4 f' h) y0

-- System derivate, f'(y), for two masses (m1, m2) where m1 is connected to origo through spring
-- with k1 coefficient. m2 is connected to m1 through spring with k2 coefficient
-- See: http://www.myphysicslab.com/dbl_spring2d.html
double_spring m1 m2 k1 k2 r1 r2 g (x1:y1:x2:y2:v1x:v1y:v2x:v2y:_) =
[ v1x, v1y, v2x, v2y,
(-1)*(k1/m1)*s1*sin(phi1) +      (k2/m1)*s2*sin(phi2),
(-1)*(k1/m1)*s1*cos(phi1) +      (k2/m1)*s2*cos(phi2) + g ,
(-1)*(k2/m2)*s2*sin(phi2)     ,
(-1)*(k2/m2)*s2*cos(phi2) + g   ]
where
l1 = sqrt(x1^2+y1^2)
l2 = sqrt((x2-x1)^2+(y2-y1)^2)
s1 = l1-r1
s2 = l2-r2
phi1 = atan(x1/y1) + addPi y1
phi2 = atan((x2-x1)/(y2-y1)) + addPi (y2-y1)

-- Helper converting radians to the degrees used by gloss

-- Helper used when determining the angle from atan(a/b)
addPi b = if (b>0) then 0 else pi

-- Show a mass m on a spring from p1 to p2
showSpring :: Float -> Float -> Point -> Point -> Picture
showSpring m r0 (x1,y1) (x2,y2) =
pictures \$ [ spring (round r0) (x1,-y1) (x2,-y2), Translate x2 (-y2) \$ Color c \$ ThickCircle r 1 ]
where
c = light \$ light black

-- Show spring from p1 to p2 n is rest length
spring :: Int -> Point -> Point -> Picture
spring n (x1,y1) (x2,y2) = Translate x1 y1 \$ Rotate (radiansToDegrees phi) \$ Scale (l/r) 1 \$ Line ps where
(dx,dy) = (x2-x1,y2-y1)
r = fromIntegral n              -- rest len
l = sqrt(dx^2+dy^2)             -- actual len
phi = atan(dy/dx) + addPi dx    -- angle
xs = [0, 0.25 .. r]
ys = take (length xs) \$ concat \$ repeat [0,0.5,0,-0.5]
ps = zipWith (\x y -> (x, y)) xs ys

main = do
let h = 0.001
let total_time = 20*60
let (r1, r2, m1, m2, k1, k2) = (5, 5, 0.1, 0.1, 1,1)
let (x1, y1, x2, y2) = (5, 0, 10, 0)
let series = take (floor (total_time/h)) (integrate (double_spring m1 m2 k1 k2 5 5 9.81)  [x1,y1,x2,y2,0,0,0,0] h)
c <- getChar
animate
(InWindow "Double Spring" (400, 400) (150, 150))
white
(doubleSpring m1 m2 r1 r2 h series)   -- function of t that produce the picture

-- Produce picture from constants and pre-calculated model-data (series)
doubleSpring m1 m2 r1 r2 h series t = Translate 0 trans \$ Scale scale scale \$ picture
where
size = 400
scale = 15
trans = 0.4*(fromIntegral size)
picture = showSystem (series !! (floor (t/h)))
showSystem (x1:y1:x2:y2:v1x:v1y:v2x:v2y:_) = pictures \$ [ ThickCircle 0.2 0.1,
showSpring m1 r1 ( 0, 0) (x1,y1),
showSpring m2 r2 (x1,y1) (x2,y2) ]
``````

A movie of the simulation can be found below.

## Try out the hamilton library

Try the hamilton library for this purpose.

## Try out the ad library

Try the ad library for this purpose.