I'm trying to write a simulator for charged and massed objects based on just calculating the net force on each object then finding the change in position across the period of time specified by the user.
However, I'm finding that when I change the dt, the change in position is drastic, when it shouldn't change significantly, decreasing the dt should just let the position converge on the correct answer.
For instance, with objects at the Cartesian coordinates (1, 0, 0) and (-1, 0, 0), with masses of 9e-31 (mass of electron) and a charge of 1 Coulomb (not the charge of an electron, I know), run for 0.1 seconds and a dt of 0.01 seconds, there is a total change of position of 2048 meters for each object. However, run for 0.1 seconds and a dt of 0.001 seconds, there is a change in position of about 1.3e30 meters. This seems rather outrageous to me, but I can't find any issues in the parts that use dt.
The code I'm using (c/p'd to avoid any possible changes)
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))
The issue I know is in the tUniverse function, probably in some calculation of either acceleration, velocity, or position (accels, vels, or poses). I've tried changing toAccel, newVeloc, and newPos by multiplying each by the inverse of dt, but it didn't significantly change the outputs.
Feel free to ignore inelasticCollide, I could probably replace it with the id function, but I just left it in because it will be relevant at some point.
EDIT: I've updated the code to fix the incorrect calculation of acceleration, the switching of mass and charge in inelasticallyCollide, and the double counting with dpos/dvel, but I'm still finding that I'm getting an error of by a magnitude of 10. For instance, with a charge of 1 C for each, I got ~10^8 for dt = 0.01 and ~10^7 for dt = 0.1 and with a charge of 0.01 C for each, ~250 for dt = 0.01 and ~65 for dt = 0.1.
It seems the "obvious" issue is that
newWidgetassumesdposanddvelare deltas, but when it's called intUniverseposesandvelshave actually already done the addition.To debug I had rewritten things to use
newtypesthinking that perhaps there was a mismatch somewhere. There did turn out to be an issue of masses and charges being transposed ininelasticCollidebut that didn't matter for my test case. The way I found this issue was by adding thetraces and seeing that the object's position component doubled each tick when the velocity component was 1.I have no idea whether any calculations are accurate otherwise.