Modelling a pendulum in Haskell using RK4

Question

I am trying to implement the classical RK4 algorithm in Haskell to solve systems of coupled ODEs, taking the simple pendulum as a sample system. Unfortunately I haven’t been able to fix the following error.

Pendulum.hs:43:27:
    No instance for (Num [Double]) arising from a use of `rk4'
    Possible fix: add an instance declaration for (Num [Double])
    In the first argument of `iterate', namely `(rk4 fs h)'
    In the expression: iterate (rk4 fs h) initial
    In an equation for `result': result = iterate (rk4 fs h) initial
Failed, modules loaded: none.

Where Haskell has inferred the type of rk4 to be

rk4  :: (Fractional b, Num [b]) => [[b] -> [b] -> b] -> b -> [[b]] -> [[b]]

The source code

module Pendulum where
import Data.List

{-
  The simple pendulum

  This is system which obeys the following equation.

  $ \frac{d^2 \theta}{dt^2} + \frac{g}{L}\theta = 0 $

  where g is the acceleration due to gravity, L the length of the pendulum and theta
  the angular displacement of the pendulum

  By observation you could see thar for small angles the solution for simple harmonic osscilator could be used,
  but for large angles a numerical solution is required.

  This simulation uses the Runge Kutta 4th Order Method to solve the ODE.

-}

data PhysicalParam = PhysicalParam { pMass   :: Double -- kg   Pendulum mass
                                   , pLength :: Double -- m    Pendulum length
                                   }
data PendulumState = PendulumState { omega   :: Double -- rad/s Angular Velocity
                                   , theta   :: Double -- rad   Angular Displacement
                                   , time    :: Double -- s     Time
                                   }
type Time = Double
type TIncrement = Double

list2State :: [[Double]] -> PendulumState
list2State [xs, ys]  = PendulumState { theta = (ys !! 0)
                                     , omega = (ys !! 1)
                                     , time  = (xs !! 0)
                                     }

state2List :: PendulumState -> [[Double]]
state2List state = [[time state],[theta state, omega state ]]


--pendulum :: PhysicalParam -> PendulumState -> Time -> TIncrement -> [PendulumState]
pendulum physicalParam initState time dt = result -- map list2State result
  where result = iterate (rk4 fs h) initial
        omegaDerivative pp state = adg * sin (theta state) / (pLength pp)
        thetaDerivative pp state = omega state
        fs   = [(\xs ys -> omegaDerivative physicalParam (list2State [xs, ys]))
               ,(\xs ys -> thetaDerivative physicalParam (list2State [xs, ys]))
               ]
        initial  = state2List initState
        h   = dt
        adg = 9.81 -- m/s

--rk4 :: Floating a => [ [a] -> [a] -> a ] -> a -> [[a]] -> [[a]]
rk4 fs h [xs, ys] = [xs', ys']
  where xs' = map (+h) xs
        ys' = ys + map (*(h/6.0)) (k1 + (map (*2.0) k2) + (map (*2.0) k3) + k4)
          where k1 = [f (xs)              (ys)                     | f <- fs]
                k2 = [f (map (+0.5*h) xs) (ys + map (*(0.5*h)) k1) | f <- fs]
                k3 = [f (map (+0.5*h) xs) (ys + map (*(0.5*h)) k2) | f <- fs]
                k4 = [f (map (+1.0*h) xs) (ys + map (*(1.0*h)) k3) | f <- fs]

physParam = PhysicalParam { pMass = 1.0, pLength = 0.1 }

-- Initial Conditions
initState = PendulumState { omega = 1.0, theta   = 1.0, time = 0.0 }

dt = 1.0/100   -- time increment
simTime = 30.0 -- simulation time

simulatedStates = pendulum physParam initState simTime dt

Answer 1

When something like Num [b] turns up in an inferred type, you can be quite sure there appears a numerical literal where GHC expects a list, or a list turns up as one of the arguments to a numerical operator. And indeed:

... ys + map (*(0.5*h)) k1 ...

both ys and k1 are lists, but + is numerical addition.

Perhaps you're trying to use ++ here, but let me tell you your whole approach is less than nice. Why do you accept [xs, ys] as a always–two-element list? Start with a reasonable signature, something similar to

rk4 :: (Floating a, Floating t) => [a -> t -> a] -> a -> [t] -> [a]`

Write that down first of all, and then start to implement the function from the outside to the inside, first leaving more local subexpressions undefined. Compile in frequent intervals, and when an error turns up you'll immediately see where.

---

BTW. It's not really great to use Haskell lists for implementing mathematical vectors, on any level. The "right" version of Runge-Kutta, IMO, has a signature along the lines

rk4 :: (VectorSpace v, t ~ Scalar v, Floating t)
     => (v -> t -> v) -> v -> [t] -> [v]

with a proper type-safe vector class.