1

我试图在 Haskell 中实现经典的 RK4 算法来解决耦合 ODE 的系统，以简单的钟摆作为示例系统。不幸的是，我无法修复以下错误。

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.

Haskell 推断出 rk4 的类型

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

源代码

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
4

1 回答 1

5

当类似的东西Num [b]出现在推断类型中时，您可以确定在 GHC 需要一个列表的地方出现一个数字文字，或者一个列表作为数字运算符的参数之一出现。事实上：

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

两者ys和k1都是列表，但是+是数字加法。

也许你想在++这里使用，但让我告诉你你的整个方法不太好。为什么您接受[xs, ys]始终包含两个元素的列表？从一个合理的签名开始，类似于

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

先把它写下来，然后从外到内开始实现函数，先留下更多的局部子表达式undefined。频繁地编译，当出现错误时，您会立即看到在哪里。

顺便提一句。在任何级别上使用 Haskell 列表来实现数学向量都不是很好。国际海事组织 Runge-Kutta 的“正确”版本具有沿线的签名

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

具有适当的类型安全向量类。

于 2013-09-23T13:01:26.397 回答