我正在 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 的力计算相关加速度，然后更新每个正如我所做的那样，ps 中的 Particle 加速。
我考虑过做类似 foldr f ps Coms 的事情。起始累加器将是当前的ps，而f将是一些函数，它接受每个comb并更新相关的Particle > 在 ps 中，并返回该累加器。对于这样一个简单的过程来说，这似乎确实是内存密集型的并且相当复杂。
有什么想法吗？
最佳答案
从 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]
...所以问题本质上是如何减少gravitate的调用次数锻炼时accsMatrix ，利用 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]
...然后，关键是 make 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²) 操作，并且正如 @leftaroundabout 提到的 O (n · (n – 1)) ≡ O (n²) ，但每次迭代都应该调用 gravitate n * (n-1) / 2次。
关于haskell - 在 Haskell 中共享力的计算，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/44622757/