在 Haskell 中共享力的计算

Question

我正在 Haskell 中实施 N 体模拟。 https://github.com/thorlucas/N-Body-Simulation

现在，每个粒子计算它的力，然后计算相对于其他粒子的加速度。换句话说，O(n²) 力的计算。如果我要计算每个组合一次，我可以将其减少到 O(n 选择 2)。

let combs = [(a, b) | (a:bs) <- tails ps, b <- bs ]
    force = map (\comb -> gravitate (fst comb) (snd comb)) combs

但我不知道如何在不使用状态的情况下将这些应用于粒子。在上面的例子中，ps 是 [Particle] 其中

data Particle = Particle Mass Pos Vel Acc deriving (Eq, Show)

理论上，在有状态的语言中，我可以简单地遍历组合，计算来自每个 a 和 b 的力的相关加速度，然后更新每个 Particle 在 ps 加速，因为我这样做。

我考虑过 foldr f ps combs 之类的事情。起始累加器将是当前的 ps，f 将是某个函数，它接受每个 comb 并更新 ps 中的相关 Particle，以及 [=41] =] 那个累加器。对于这样一个简单的过程，这似乎确实占用大量内存并且相当复杂。

有什么想法吗？

Answer 1

从 https://github.com/thorlucas/N-Body-Simulation

中获取代码

updateParticle :: Model -> Particle -> Particle
updateParticle ps p@(Particle m pos vel acc) =
    let accs = map (gravitate p) ps
        acc' = foldr (\(accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs
        vel' = (fst vel + fst acc, snd vel + snd acc)
        pos' = (fst pos + fst vel, snd pos + snd vel)
    in  Particle m pos' vel' acc'

step :: ViewPort -> Float -> Model -> Model
step _ _ ps = map (updateParticle ps) ps

并对其进行修改，以便在矩阵（好吧，列表的列表...）中计算加速度，而不是更新每个粒子，我们得到...

updateParticle :: Model -> (Particle, [Acc]) -> Particle
updateParticle ps (p@(Particle m pos vel acc), accs) =
    let acc' = foldr (\(accx, accy) (x, y) -> (accx + x, accy + y)) (0, 0) accs
        vel' = (fst vel + fst acc, snd vel + snd acc)
        pos' = (fst pos + fst vel, snd pos + snd vel)
    in  Particle m pos' vel' acc'

step :: ViewPort -> Float -> Model -> Model
step _ _ ps = map (updateParticle ps) $ zip ps accsMatrix where
    accsMatrix = [map (gravitate p) ps | p <- ps]

... 所以问题本质上是如何在计算 accsMatrix 时减少对 gravitate 的调用次数，利用 gravitate a b = -1 * gravitate b a 这一事实.

如果我们要打印出来 accsMatrix，它看起来像...

[[( 0.0,  0.0), ( 1.0,  2.3), (-1.0,  0.0), ...
[[(-1.0, -2.3), ( 0.0,  0.0), (-1.2,  5.3), ...
[[( 1.0,  0.0), ( 1.2, -5.3), ( 0.0,  0.0), ...
...

...所以我们看到 accsMatrix !! i !! j == -1 * accsMatrix !! j !! i.

因此，要使用上述事实，我们需要访问一些索引。首先，我们索引外部列表...

accsMatrix = [map (gravitate p) ps | (i,p) <- zip [0..] ps]

...并用列表理解替换内部列表...

accsMatrix = [[ gravitate p p' | p' <- ps] | (i,p) <- zip [0..] ps]

...通过 zip 获取更多索引...

accsMatrix = [[ gravitate p p' | (j, p') <- zip [0..] ps] | (i,p) <- zip [0..] ps]

... 然后，关键是让 accsMatrix 依赖于矩阵的一半...

accsMatrix = [[ if i == j then 0 else if i < j then gravitate p p' else -1 * accsMatrix !! j !! i | (j, p') <- zip [0..] ps] | (i, p) <- zip [0..] ps]

我们也可以拆分一下，如下所示...

accsMatrix = [[ accs (j, p') (i, p) | (j, p') <- zip [0..] ps] | (i, p) <- zip [0..] ps]
accs (j, p') (i, p) 
    | i == j    = 0
    | i < j     = gravitate p p'
    | otherwise = -1 * accsMatrix !! j !! i

... 或使用 map

避免列表理解

accsMatrix = map (flip map indexedPs) $ map accs indexedPs  
indexedPs = zip [0..] ps
accs (i, p) (j, p')
    | i == j    = 0
    | i < j     = gravitate p p'
    | otherwise = -1 * accsMatrix !! j !! i

...或使用列表 monad...

accsMatrix = map accs indexedPs >>= (:[]) . flip map indexedPs

...虽然（对我而言）更难看出其中发生了什么。

这种列表方法可能存在一些可怕的性能问题，尤其是使用 !!，而且您仍然运行 O(n²) 操作由于遍历，以及 O (n · (n – 1)) ≡ O (n²) 正如@leftaroundabout 提到的，但每次迭代应该调用 gravitate n * (n-1) / 2 次.

在 Haskell 中共享力的计算

Sharing Computation of Force in Haskell

simulation

haskell

physics