通过净力计算位置，不同的dts得出不同的答案
Calculating position via net force, different dts yield different answers
我正在尝试为带电和有质量的物体编写一个模拟器，该模拟器仅基于计算每个物体上的净力，然后找出用户指定时间段内的位置变化。
然而，我发现当我改变dt时，position的变化是剧烈的，当它不应该有明显变化时，减小dt应该只是让position收敛于正确答案。
例如，对于笛卡尔坐标 (1, 0, 0) 和 (-1, 0, 0) 处的物体，其质量为 9e-31（电子质量）且电荷为 1 库仑（不是一个电子的电荷，我知道），运行 0.1 秒和 0.01 秒的 dt，每个物体的位置总变化为 2048 米。但是，运行0.1秒，0.001秒的dt，位置有1.3e30米左右的变化。这对我来说似乎相当离谱，但我在使用 dt 的部分中找不到任何问题。
我正在使用的代码（c/p 以避免任何可能的更改）
import Data.List
main = print $ mainprog
where
mainprog = runUniverse makeUniverse 1 0.1
type Length = Double
type Mass = Double
type Charge = Double
type Time = Double
type Vector = (Double, Double, Double)
type Position = Vector
type Velocity = Vector
type Acceleration = Vector
type Force = Vector
data Widget = Widget {pos :: Position, mass :: Double, charge :: Double, velocity :: Velocity} deriving (Eq, Show, Read)
--utils
toScalar :: Vector -> Double
toScalar (x, y, z) = sqrt (x ^^ 2 + y ^^ 2 + z ^^ 2)
toUnit :: Vector -> Vector
toUnit (x, y, z) = (x / scalar, y / scalar, z / scalar)
where
scalar = toScalar (x, y, z)
add :: Vector -> Vector -> Vector
add (x1, y1, z1) (x2, y2, z2) = (x1 + x2, y1 + y2, z1 + z2)
mult :: Vector -> Double -> Vector
mult (x, y, z) k = (k * x, k * y, k * z)
diff :: Vector -> Vector -> Vector
diff (x1, y1, z1) (x2, y2, z2) = (x1 - x2, y1 - y2, z1 - z2)
--calcs
gForce :: Widget -> Widget -> Force
gForce (Widget pos1 mass1 _ _) (Widget pos2 mass2 _ _) = mult unitForce scalarForce
where
unitForce = toUnit posdiff
scalarForce = (g * mass1 * mass2) / (radius ^^ 2)
g = 6.674e-11
radius = toScalar posdiff
posdiff = diff pos1 pos2
eForce :: Widget -> Widget -> Force
eForce (Widget pos1 _ charge1 _) (Widget pos2 _ charge2 _) = mult unitForce scalarForce
where
unitForce = (toUnit posdiff)
--necessary to determine attraction vs repulsion, whereas gravitational is always attractive
scalarForce = ((abs (k_C * charge1 * charge2)) / (radius ^^ 2)) * (signum charge1) * (signum charge2)
k_C = 8.988e9
radius = toScalar posdiff
posdiff = diff pos1 pos2
netForce :: [Force] -> Force
netForce = foldl add (0, 0, 0)
toAccel :: Force -> Widget -> Acceleration
toAccel f (Widget _ mass _ _) = mult f (1/mass)
newVeloc :: Velocity -> Acceleration -> Time -> Velocity
newVeloc v a dt = add v (mult a dt)
newPos :: Vector -> Velocity -> Time -> Vector
newPos s v dt = add s (mult v dt)
newWidget :: Widget -> Position -> Velocity -> Widget
newWidget (Widget pos1 mass charge vel1) pos2 vel2 = Widget pos2 mass charge vel2
tUniverse :: [Widget] -> Time -> [Widget]
tUniverse widgets dt = zipWith3 newWidget widgets poses vels
where
netMassForces = map (\w -> gForcePrime w (widgets \ [w])) widgets
gForcePrime w ws = netForce $ map (gForce w) ws
netElectricForces = map (\w -> eForcePrime w (widgets \ [w])) widgets
eForcePrime w ws = netForce $ map (eForce w) ws
volds = map velocity widgets
polds = map pos widgets
accels = zipWith toAccel (map netForce (zipWith (\a b -> a : [b]) netMassForces netElectricForces)) widgets
vels = zipWith (\v a -> newVeloc v a dt) volds accels
poses = zipWith (\s v -> newPos s v dt) polds vels
makeUniverse :: [Widget]
makeUniverse = [(Widget (-1, 0, 0) 1 1 (0, 0, 0)), (Widget (1, 0, 0) 1 1 (0, 0, 0))]
runUniverse :: [Widget] -> Time -> Time -> [Widget]
runUniverse ws t dt
| t <= 0 = ws
| otherwise = runUniverse (tUniverse (inelasticCollide ws) dt) (t-dt) dt
inelasticCollide :: [Widget] -> [Widget]
inelasticCollide [] = []
inelasticCollide (w:[]) = [w]
inelasticCollide (w:ws) = (combine w (sameposes w ws)) : (inelasticCollide $ ws \ (sameposes w ws))
where
sameposes w ws = filter (\w' -> pos w == pos w') ws
combine :: Widget -> [Widget] -> Widget
combine = foldl (\(Widget pos mass1 charge1 veloc1) (Widget _ mass2 charge2 veloc2) -> Widget pos (charge1 + charge2) (mass1 + mass2) (newveloc mass1 mass2 veloc1 veloc2))
--inelastic collision, m1v1 + m2v2 = m3v3 therefore v3 = (m1v1 + m2v2)/(m1 + m2)
newveloc m1 m2 v1 v2 = ((v1 `mult` m1) `add` (v2 `mult` m2)) `mult` (1 / (m1 + m2))
我知道的问题出在 tUniverse 函数中，可能是在计算加速度、速度或位置（加速度、速度或位姿）时。我已经尝试通过将每个乘以 dt 的倒数来更改 toAccel、newVeloc 和 newPos，但它并没有显着改变输出。
请随意忽略 inelasticCollide，我可能会用 id 函数替换它，但我只是将它留在原处，因为它在某些时候是相关的。
编辑：我已经更新了代码以修复加速度的不正确计算、非弹性碰撞中质量和电荷的切换以及 dpos/dvel 的重复计算，但我仍然发现我' m 得到 10 幅度的误差。例如，每个充电 1 C，对于 dt = 0.01，我得到 ~10^8，对于 dt = 0.1，我得到 ~10^7，并且充电 0.01 C每个，dt = 0.01 时约 250，dt = 0.1 时约 65。
似乎 "obvious" 问题是 newWidget 假设 dpos 和 dvel 是增量，但是当它在 tUniverse 中调用时 poses和vels其实已经加完了
为了调试，我重写了要使用的内容 newtypes，我认为可能某处不匹配。 inelasticCollide 中确实存在质量和电荷被调换的问题，但这对我的测试用例无关紧要。我发现这个问题的方法是添加 traces 并看到对象的位置分量在速度分量为 1 时每个刻度都加倍。
我不知道其他计算是否准确。
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
import Data.List
import Debug.Trace (trace)
main = print $ runUniverse makeUniverse 0.1 0.01
newtype Length = Length {unLength::Double}
newtype Mass = Mass {unMass::Double} deriving (Num,Eq,Show)
newtype Charge = Charge {unCharge::Double} deriving (Num,Eq,Show)
newtype Time = Time {unTime::Double} deriving (Num,Eq,Ord,Fractional)
type Vector = (Double,Double,Double)
newtype Position = Position {unPosition::Vector} deriving (Eq,Show)
newtype Velocity = Velocity {unVelocity::Vector} deriving (Eq,Show)
newtype Acceleration = Acceleration {unAcceleration::Vector}
newtype Force = Force {unForce::Vector} deriving (Eq,Show)
data Widget = Widget {pos :: Position, mass :: Mass, charge :: Charge, velocity :: Velocity} deriving (Eq, Show)
--utils
toScalar :: Vector -> Double
toScalar (x, y, z) = sqrt (x ^^ 2 + y ^^ 2 + z ^^ 2)
toUnit :: Vector -> Vector
toUnit (x, y, z) = (x / scalar, y / scalar, z / scalar)
where
scalar = toScalar (x, y, z)
add :: Vector -> Vector -> Vector
add (x1, y1, z1) (x2, y2, z2) = (x1 + x2, y1 + y2, z1 + z2)
mult :: Vector -> Double -> Vector
mult (x, y, z) k = (k * x, k * y, k * z)
diff :: Vector -> Vector -> Vector
diff (x1, y1, z1) (x2, y2, z2) = (x1 - x2, y1 - y2, z1 - z2)
--calcs
gForce :: Widget -> Widget -> Force
gForce (Widget (Position pos1) (Mass mass1) _ _) (Widget (Position pos2) (Mass mass2) _ _) = Force (mult unitForce scalarForce)
where
unitForce = toUnit posdiff
scalarForce = (g * mass1 * mass2) / (radius ^^ 2)
g = 6.674e-11
radius = toScalar posdiff
posdiff = diff pos1 pos2
eForce :: Widget -> Widget -> Force
eForce (Widget (Position pos1) _ (Charge charge1) _) (Widget (Position pos2) _ (Charge charge2) _) = Force (mult unitForce scalarForce)
where
unitForce = (toUnit posdiff)
--necessary to determine attraction vs repulsion, whereas gravitational is always attractive
scalarForce = ((abs (k_C * charge1 * charge2)) / (radius ^^ 2)) * (signum charge1) * (signum charge2)
k_C = 8.988e9
radius = toScalar posdiff
posdiff = diff pos1 pos2
netForce :: [Force] -> Force
netForce = Force . foldl add (0,0,0) . map unForce
toAccel :: Force -> Widget -> Acceleration
toAccel f (Widget _ mass _ _) = Acceleration (mult (unForce f) (unMass mass))
newVeloc :: Velocity -> Acceleration -> Time -> Velocity
newVeloc v a dt = Velocity (add (unVelocity v) (mult (unAcceleration a) (unTime dt)))
newPos :: Position -> Velocity -> Time -> Position
newPos s v dt = Position (add (unPosition s) (mult (unVelocity v) (unTime dt)))
newWidget :: Widget -> Position -> Velocity -> Widget
newWidget w@(Widget pos1 _ _ vel1) dpos dvel = w { pos=Position ((unPosition dpos)),velocity=Velocity ((unVelocity dvel)) }
tUniverse :: [Widget] -> Time -> [Widget]
tUniverse widgets dt = zipWith3 newWidget widgets (trace (show poses) poses) (trace (show vels) vels)
where
netMassForces = map (\w -> gForcePrime w (widgets \ [w])) widgets
gForcePrime w ws = netForce $ map (gForce w) ws
netElectricForces = map (\w -> eForcePrime w (widgets \ [w])) widgets
eForcePrime w ws = netForce $ map (eForce w) ws
volds = map velocity widgets
polds = map pos widgets
accels = zipWith toAccel (map netForce (zipWith (\a b -> a : [b]) netMassForces netElectricForces)) widgets
vels = zipWith (\v a -> newVeloc v a dt) volds accels
poses = zipWith (\s v -> newPos s v dt) polds vels
makeUniverse :: [Widget]
makeUniverse = [Widget (Position (1,0,0)) (Mass 0) (Charge 0) (Velocity (1,0,0))] -- , (Widget (1, 0, 0) 9e-31 1 (0, 0, 0))]
runUniverse :: [Widget] -> Time -> Time -> [Widget]
runUniverse ws t dt
| t < 0 = ws
| otherwise = runUniverse (tUniverse (inelasticCollide ws) dt) (t-dt) dt
inelasticCollide :: [Widget] -> [Widget]
inelasticCollide [] = []
inelasticCollide (w:[]) = [w]
inelasticCollide (w:ws) = (combine w (sameposes w ws)) : (inelasticCollide $ ws \ (sameposes w ws))
where
sameposes w ws = filter (\w' -> pos w == pos w') ws
combine :: Widget -> [Widget] -> Widget
combine = foldl (\(Widget pos mass1 charge1 veloc1) (Widget _ mass2 charge2 veloc2) -> Widget pos (mass1 + mass2) (charge1 + charge2) (Velocity (newveloc (unMass mass1) (unMass mass2) (unVelocity veloc1) (unVelocity veloc2))))
--inelastic collision, m1v1 + m2v2 = m3v3 therefore v3 = (m1v1 + m2v2)/(m1 + m2)
newveloc m1 m2 v1 v2 = ((v1 `mult` m1) `add` (v2 `mult` m2)) `mult` (1 / (m1 + m2))
我正在尝试为带电和有质量的物体编写一个模拟器，该模拟器仅基于计算每个物体上的净力，然后找出用户指定时间段内的位置变化。
然而，我发现当我改变dt时，position的变化是剧烈的，当它不应该有明显变化时，减小dt应该只是让position收敛于正确答案。
例如，对于笛卡尔坐标 (1, 0, 0) 和 (-1, 0, 0) 处的物体，其质量为 9e-31（电子质量）且电荷为 1 库仑（不是一个电子的电荷，我知道），运行 0.1 秒和 0.01 秒的 dt，每个物体的位置总变化为 2048 米。但是，运行0.1秒，0.001秒的dt，位置有1.3e30米左右的变化。这对我来说似乎相当离谱，但我在使用 dt 的部分中找不到任何问题。
我正在使用的代码（c/p 以避免任何可能的更改）
import Data.List
main = print $ mainprog
where
mainprog = runUniverse makeUniverse 1 0.1
type Length = Double
type Mass = Double
type Charge = Double
type Time = Double
type Vector = (Double, Double, Double)
type Position = Vector
type Velocity = Vector
type Acceleration = Vector
type Force = Vector
data Widget = Widget {pos :: Position, mass :: Double, charge :: Double, velocity :: Velocity} deriving (Eq, Show, Read)
--utils
toScalar :: Vector -> Double
toScalar (x, y, z) = sqrt (x ^^ 2 + y ^^ 2 + z ^^ 2)
toUnit :: Vector -> Vector
toUnit (x, y, z) = (x / scalar, y / scalar, z / scalar)
where
scalar = toScalar (x, y, z)
add :: Vector -> Vector -> Vector
add (x1, y1, z1) (x2, y2, z2) = (x1 + x2, y1 + y2, z1 + z2)
mult :: Vector -> Double -> Vector
mult (x, y, z) k = (k * x, k * y, k * z)
diff :: Vector -> Vector -> Vector
diff (x1, y1, z1) (x2, y2, z2) = (x1 - x2, y1 - y2, z1 - z2)
--calcs
gForce :: Widget -> Widget -> Force
gForce (Widget pos1 mass1 _ _) (Widget pos2 mass2 _ _) = mult unitForce scalarForce
where
unitForce = toUnit posdiff
scalarForce = (g * mass1 * mass2) / (radius ^^ 2)
g = 6.674e-11
radius = toScalar posdiff
posdiff = diff pos1 pos2
eForce :: Widget -> Widget -> Force
eForce (Widget pos1 _ charge1 _) (Widget pos2 _ charge2 _) = mult unitForce scalarForce
where
unitForce = (toUnit posdiff)
--necessary to determine attraction vs repulsion, whereas gravitational is always attractive
scalarForce = ((abs (k_C * charge1 * charge2)) / (radius ^^ 2)) * (signum charge1) * (signum charge2)
k_C = 8.988e9
radius = toScalar posdiff
posdiff = diff pos1 pos2
netForce :: [Force] -> Force
netForce = foldl add (0, 0, 0)
toAccel :: Force -> Widget -> Acceleration
toAccel f (Widget _ mass _ _) = mult f (1/mass)
newVeloc :: Velocity -> Acceleration -> Time -> Velocity
newVeloc v a dt = add v (mult a dt)
newPos :: Vector -> Velocity -> Time -> Vector
newPos s v dt = add s (mult v dt)
newWidget :: Widget -> Position -> Velocity -> Widget
newWidget (Widget pos1 mass charge vel1) pos2 vel2 = Widget pos2 mass charge vel2
tUniverse :: [Widget] -> Time -> [Widget]
tUniverse widgets dt = zipWith3 newWidget widgets poses vels
where
netMassForces = map (\w -> gForcePrime w (widgets \ [w])) widgets
gForcePrime w ws = netForce $ map (gForce w) ws
netElectricForces = map (\w -> eForcePrime w (widgets \ [w])) widgets
eForcePrime w ws = netForce $ map (eForce w) ws
volds = map velocity widgets
polds = map pos widgets
accels = zipWith toAccel (map netForce (zipWith (\a b -> a : [b]) netMassForces netElectricForces)) widgets
vels = zipWith (\v a -> newVeloc v a dt) volds accels
poses = zipWith (\s v -> newPos s v dt) polds vels
makeUniverse :: [Widget]
makeUniverse = [(Widget (-1, 0, 0) 1 1 (0, 0, 0)), (Widget (1, 0, 0) 1 1 (0, 0, 0))]
runUniverse :: [Widget] -> Time -> Time -> [Widget]
runUniverse ws t dt
| t <= 0 = ws
| otherwise = runUniverse (tUniverse (inelasticCollide ws) dt) (t-dt) dt
inelasticCollide :: [Widget] -> [Widget]
inelasticCollide [] = []
inelasticCollide (w:[]) = [w]
inelasticCollide (w:ws) = (combine w (sameposes w ws)) : (inelasticCollide $ ws \ (sameposes w ws))
where
sameposes w ws = filter (\w' -> pos w == pos w') ws
combine :: Widget -> [Widget] -> Widget
combine = foldl (\(Widget pos mass1 charge1 veloc1) (Widget _ mass2 charge2 veloc2) -> Widget pos (charge1 + charge2) (mass1 + mass2) (newveloc mass1 mass2 veloc1 veloc2))
--inelastic collision, m1v1 + m2v2 = m3v3 therefore v3 = (m1v1 + m2v2)/(m1 + m2)
newveloc m1 m2 v1 v2 = ((v1 `mult` m1) `add` (v2 `mult` m2)) `mult` (1 / (m1 + m2))
我知道的问题出在 tUniverse 函数中，可能是在计算加速度、速度或位置（加速度、速度或位姿）时。我已经尝试通过将每个乘以 dt 的倒数来更改 toAccel、newVeloc 和 newPos，但它并没有显着改变输出。
请随意忽略 inelasticCollide，我可能会用 id 函数替换它，但我只是将它留在原处，因为它在某些时候是相关的。
编辑：我已经更新了代码以修复加速度的不正确计算、非弹性碰撞中质量和电荷的切换以及 dpos/dvel 的重复计算，但我仍然发现我' m 得到 10 幅度的误差。例如，每个充电 1 C，对于 dt = 0.01，我得到 ~10^8，对于 dt = 0.1，我得到 ~10^7，并且充电 0.01 C每个，dt = 0.01 时约 250，dt = 0.1 时约 65。
似乎 "obvious" 问题是 newWidget 假设 dpos 和 dvel 是增量，但是当它在 tUniverse 中调用时 poses和vels其实已经加完了
为了调试，我重写了要使用的内容 newtypes，我认为可能某处不匹配。 inelasticCollide 中确实存在质量和电荷被调换的问题，但这对我的测试用例无关紧要。我发现这个问题的方法是添加 traces 并看到对象的位置分量在速度分量为 1 时每个刻度都加倍。
我不知道其他计算是否准确。
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
import Data.List
import Debug.Trace (trace)
main = print $ runUniverse makeUniverse 0.1 0.01
newtype Length = Length {unLength::Double}
newtype Mass = Mass {unMass::Double} deriving (Num,Eq,Show)
newtype Charge = Charge {unCharge::Double} deriving (Num,Eq,Show)
newtype Time = Time {unTime::Double} deriving (Num,Eq,Ord,Fractional)
type Vector = (Double,Double,Double)
newtype Position = Position {unPosition::Vector} deriving (Eq,Show)
newtype Velocity = Velocity {unVelocity::Vector} deriving (Eq,Show)
newtype Acceleration = Acceleration {unAcceleration::Vector}
newtype Force = Force {unForce::Vector} deriving (Eq,Show)
data Widget = Widget {pos :: Position, mass :: Mass, charge :: Charge, velocity :: Velocity} deriving (Eq, Show)
--utils
toScalar :: Vector -> Double
toScalar (x, y, z) = sqrt (x ^^ 2 + y ^^ 2 + z ^^ 2)
toUnit :: Vector -> Vector
toUnit (x, y, z) = (x / scalar, y / scalar, z / scalar)
where
scalar = toScalar (x, y, z)
add :: Vector -> Vector -> Vector
add (x1, y1, z1) (x2, y2, z2) = (x1 + x2, y1 + y2, z1 + z2)
mult :: Vector -> Double -> Vector
mult (x, y, z) k = (k * x, k * y, k * z)
diff :: Vector -> Vector -> Vector
diff (x1, y1, z1) (x2, y2, z2) = (x1 - x2, y1 - y2, z1 - z2)
--calcs
gForce :: Widget -> Widget -> Force
gForce (Widget (Position pos1) (Mass mass1) _ _) (Widget (Position pos2) (Mass mass2) _ _) = Force (mult unitForce scalarForce)
where
unitForce = toUnit posdiff
scalarForce = (g * mass1 * mass2) / (radius ^^ 2)
g = 6.674e-11
radius = toScalar posdiff
posdiff = diff pos1 pos2
eForce :: Widget -> Widget -> Force
eForce (Widget (Position pos1) _ (Charge charge1) _) (Widget (Position pos2) _ (Charge charge2) _) = Force (mult unitForce scalarForce)
where
unitForce = (toUnit posdiff)
--necessary to determine attraction vs repulsion, whereas gravitational is always attractive
scalarForce = ((abs (k_C * charge1 * charge2)) / (radius ^^ 2)) * (signum charge1) * (signum charge2)
k_C = 8.988e9
radius = toScalar posdiff
posdiff = diff pos1 pos2
netForce :: [Force] -> Force
netForce = Force . foldl add (0,0,0) . map unForce
toAccel :: Force -> Widget -> Acceleration
toAccel f (Widget _ mass _ _) = Acceleration (mult (unForce f) (unMass mass))
newVeloc :: Velocity -> Acceleration -> Time -> Velocity
newVeloc v a dt = Velocity (add (unVelocity v) (mult (unAcceleration a) (unTime dt)))
newPos :: Position -> Velocity -> Time -> Position
newPos s v dt = Position (add (unPosition s) (mult (unVelocity v) (unTime dt)))
newWidget :: Widget -> Position -> Velocity -> Widget
newWidget w@(Widget pos1 _ _ vel1) dpos dvel = w { pos=Position ((unPosition dpos)),velocity=Velocity ((unVelocity dvel)) }
tUniverse :: [Widget] -> Time -> [Widget]
tUniverse widgets dt = zipWith3 newWidget widgets (trace (show poses) poses) (trace (show vels) vels)
where
netMassForces = map (\w -> gForcePrime w (widgets \ [w])) widgets
gForcePrime w ws = netForce $ map (gForce w) ws
netElectricForces = map (\w -> eForcePrime w (widgets \ [w])) widgets
eForcePrime w ws = netForce $ map (eForce w) ws
volds = map velocity widgets
polds = map pos widgets
accels = zipWith toAccel (map netForce (zipWith (\a b -> a : [b]) netMassForces netElectricForces)) widgets
vels = zipWith (\v a -> newVeloc v a dt) volds accels
poses = zipWith (\s v -> newPos s v dt) polds vels
makeUniverse :: [Widget]
makeUniverse = [Widget (Position (1,0,0)) (Mass 0) (Charge 0) (Velocity (1,0,0))] -- , (Widget (1, 0, 0) 9e-31 1 (0, 0, 0))]
runUniverse :: [Widget] -> Time -> Time -> [Widget]
runUniverse ws t dt
| t < 0 = ws
| otherwise = runUniverse (tUniverse (inelasticCollide ws) dt) (t-dt) dt
inelasticCollide :: [Widget] -> [Widget]
inelasticCollide [] = []
inelasticCollide (w:[]) = [w]
inelasticCollide (w:ws) = (combine w (sameposes w ws)) : (inelasticCollide $ ws \ (sameposes w ws))
where
sameposes w ws = filter (\w' -> pos w == pos w') ws
combine :: Widget -> [Widget] -> Widget
combine = foldl (\(Widget pos mass1 charge1 veloc1) (Widget _ mass2 charge2 veloc2) -> Widget pos (mass1 + mass2) (charge1 + charge2) (Velocity (newveloc (unMass mass1) (unMass mass2) (unVelocity veloc1) (unVelocity veloc2))))
--inelastic collision, m1v1 + m2v2 = m3v3 therefore v3 = (m1v1 + m2v2)/(m1 + m2)
newveloc m1 m2 v1 v2 = ((v1 `mult` m1) `add` (v2 `mult` m2)) `mult` (1 / (m1 + m2))