physicsΒΆ
InΒ [Β ]:
#!import ../../../polyglot/lib/fsharp/Plotting.dib
InΒ [Β ]:
//// test
open testing
init_seriesΒΆ
InΒ [Β ]:
//// test
inl x = am'.init_series -3f64 3 0.01
inl y = x |> am'.map_base math.square
"square", "x", "y", ;[ "square", x, y ]
InΒ [Β ]:
//// test
inl x = am'.init_series -10f64 10 0.1
inl y_sin = x |> am'.map_base sin
inl y_cos = x |> am'.map_base cos
"sin cos", "x", "y", ;[ "sin", x, y_sin; "cos", x, y_cos ]
InΒ [Β ]:
//// test
inl y_pos y0 vy0 ay t =
y0 + vy0 * t + ay * (t |> math.square) / 2
inl x = am'.init_series 0f64 5 0.01
inl y = x |> am'.map_base (y_pos 0 20 -9.8)
"projectile motion", "time (s)", "", ;[ "height of projectile (m)", x, y ]
velocity_cfΒΆ
InΒ [Β ]:
type mass = f64
type time = f64
type position = f64
type velocity = f64
type force = f64
type velocity_cf = mass -> velocity -> list force -> (time -> velocity)
inl velocity_cf m v0 fs =
inl f_net = fs |> listm'.sum
inl a0 = f_net / m
inl v t = v0 + a0 * t
v
InΒ [Β ]:
//// test
velocity_cf 0.1f64 0.6 [ 0.04; -0.08 ] 0
|> _assert_eq 0.6
velocity_cf 0.1f64 0.6 [ 0.04; -0.08 ] 1
|> _assert_eq 0.2
{ name = __assert_eq; actual = +0.600000; expected = +0.600000 }
{ name = __assert_eq; actual = +0.200000; expected = +0.200000 }
InΒ [Β ]:
//// test
inl x = am'.init_series 0f64 4 0.1
inl y = x |> am'.map_base (velocity_cf 0.1f64 0.6 [ 0.04; -0.08 ])
"car on an air track", "time (s)", "", ;[ "velocity of car (m/s)", x, y ]
derivativeΒΆ
InΒ [Β ]:
type derivative = (f64 -> f64) -> f64 -> f64
inl derivative dt : derivative =
fun x t =>
(x (t + dt / 2) - x (t - dt / 2)) / dt
InΒ [Β ]:
//// test
derivative 1 (fun x => x ** 4 / 4) 1 - 1
|> _assert_approx_eq None 0.25
derivative 0.001 (fun x => x ** 4 / 4) 1 - 1
|> _assert_approx_eq None 0.0000002499998827953931
derivative 0.000001 (fun x => x ** 4 / 4) 1 - 1
|> _assert_approx_eq None 0.000000000001000088900582341
derivative 0.000000001 (fun x => x ** 4 / 4) 1 - 1
|> _assert_approx_eq None 0.00000008274037099909037
derivative 0.000000000001 (fun x => x ** 4 / 4) 1 - 1
|> _assert_approx_eq None 0.00008890058234101161
derivative 0.000000000000001 (fun x => x ** 4 / 4) 1 - 1
|> _assert_approx_eq None -0.0007992778373592246
derivative 0.000000000000000001 (fun x => x ** 4 / 4) 1 - 1
|> _assert_approx_eq None -1
{ name = __assert_approx_eq; actual = +0.250000; expected = +0.250000 }
{ name = __assert_approx_eq; actual = +0.000000; expected = +0.000000 }
{ name = __assert_approx_eq; actual = +0.000000; expected = +0.000000 }
{ name = __assert_approx_eq; actual = +0.000000; expected = +0.000000 }
{ name = __assert_approx_eq; actual = +0.000089; expected = +0.000089 }
{ name = __assert_approx_eq; actual = -0.000799; expected = -0.000799 }
{ name = __assert_approx_eq; actual = -1.000000; expected = -1.000000 }
integrationΒΆ
InΒ [Β ]:
type integration = (f64 -> f64) -> f64 -> f64 -> f64
inl integral dt : integration =
fun f a b =>
inl rec 루ν t y =
if t < b
then 루ν (t + dt) (y + f t * dt)
else t, y
루ν (a + dt / 2) 0
|> snd
InΒ [Β ]:
//// test
integral 0.01 math.square 0 1
|> _assert_approx_eq None 0.33332500000000004
{ name = __assert_approx_eq; actual = +0.333325; expected = +0.333325 }
InΒ [Β ]:
inl integral' dt : integration =
fun f a b =>
listm'.init_series (a + dt / 2) (b - dt / 2) dt
|> listm.map (f >> (*) dt)
|> listm'.sum
InΒ [Β ]:
//// test
integral' 0.1 math.square 0 1
|> _assert_approx_eq None (integral 0.1 math.square 0 1)
{ name = __assert_approx_eq; actual = +0.332500; expected = +0.332500 }
InΒ [Β ]:
inl integral'' dt : integration =
fun f x y =>
am'.init_series (x + dt / 2) (y - dt / 2) dt
|> fun x => a x : _ int _
|> am.map (f >> (*) dt)
|> am'.sum
InΒ [Β ]:
//// test
integral'' 0.01 math.square 0 1
|> _assert_approx_eq None (integral 0.01 math.square 0 1)
{ name = __assert_approx_eq; actual = +0.333325; expected = +0.333325 }
anti_derivativeΒΆ
InΒ [Β ]:
inl anti_derivative dt v0 a t =
v0 + integral' dt a 0 t
velocity_ftΒΆ
InΒ [Β ]:
type velocity_ft = mass -> velocity -> list (time -> force) -> (time -> velocity)
inl velocity_ft dt : velocity_ft =
fun m v0 fs =>
inl f_net t = fs |> listm.map (fun f => f t) |> listm'.sum
inl a t = f_net t / m
anti_derivative dt v0 a
position_ftΒΆ
InΒ [Β ]:
type position_ft = mass -> position -> velocity -> list (time -> force) -> (time -> position)
inl position_ft dt : position_ft =
fun m x0 v0 fs =>
velocity_ft dt m v0 fs
|> anti_derivative dt x0
InΒ [Β ]:
//// test
inl pedal_coast (t : time) : force =
inl t_cycle = 20
inl n_complete : i32 = t / t_cycle |> conv
inl remainder = t - conv n_complete * t_cycle
if remainder > 0 && remainder < 10
then 10
else 0
inl x = am'.init_series -5f64 45 0.1
inl y = x |> am'.map_base pedal_coast
"child pedaling then coasting", "time (s)", "", ;[ "force on bike (N)", x, y ]
InΒ [Β ]:
//// test
inl x = am'.init_series -5 45 1
inl y = x |> am'.map_base (position_ft 0.1f64 20 0 0 [ pedal_coast ])
"child pedaling then coasting", "time (s)", "", ;[ "position of bike (m)", x, y ]
velocity_fvΒΆ
InΒ [Β ]:
inl newton_second_v m fs v0 =
fs |> listm.map (fun f => f v0) |> listm'.sum |> fun x => x / m
inl update_velocity dt m fs v0 =
v0 + newton_second_v m fs v0 * dt
inl velocity_fv dt m v0 fs t =
stream.iterate (update_velocity dt m fs) v0
|> stream.try_item (t / dt |> math.round |> abs)
|> optionm'.default_value 0
InΒ [Β ]:
inl f_air drag rho area v =
-drag * rho * area * abs v * v / 2
InΒ [Β ]:
//// test
inl x = am'.init_series 0 60 0.5
inl y = x |> am'.map_base (velocity_fv 1 70 0f64 [ fun _ => 100; f_air 2 1.225 0.6 ])
"bike velocity", "time (s)", "", ;[ "velocity of bike (m/s)", x, y ]
velocity_ftvΒΆ
InΒ [Β ]:
inl newton_second_tv m fs (t, v0) =
inl f_net = fs |> listm.map (fun f => f (t, v0)) |> listm'.sum
inl acc = f_net / m
1, acc
inl update_tv dt m fs (t, v0) =
inl dtdt, dvdt = newton_second_tv m fs (t, v0)
t + dtdt * dt, v0 + dvdt * dt
inl velocity_ftv dt m tv0 fs t =
stream.iterate (join update_tv dt m fs) tv0
|> stream.try_item (t / dt |> math.round |> abs)
|> optionm.map snd
|> optionm'.default_value 0
InΒ [Β ]:
//// test
inl x = am'.init_series 0 100 0.1
inl y =
x
|> am'.map_base (
velocity_ftv 0.1 20 (dyn (0, 0)) [ fun (t, _) => pedal_coast t; fun (_, v) => f_air 2 1.225 0.5 v ]
)
"pedaling and coasting with air", "time (s)", "", ;[ "velocity of bike (m/s)", x, y ]
velocity_ftxvΒΆ
InΒ [Β ]:
nominal state_1d = time * position * velocity
nominal rrr = f64 * f64 * f64
inl newton_second_1d m fs (state_1d (t, x0, v0)) =
inl f_net = fs |> listm.map (fun f => f (state_1d (t, x0, v0))) |> listm'.sum
inl acc = f_net / m
rrr (1f64, v0, acc)
inl euler_1d dt deriv (state_1d (t0, x0, v0) as t) =
inl (rrr (_, _, dvdt)) = deriv t
inl t1 = t0 + dt
inl x1 = x0 + v0 * dt
inl v1 = v0 + dvdt * dt
state_1d (t1, x1, v1)
inl update_txv dt m fs =
newton_second_1d m fs |> euler_1d dt
inl states_txv dt m txv0 fs =
seq.iterate_ (update_txv dt m fs) txv0
inl velocity_1d sts t =
inl (state_1d (t0, _, _)) = sts 0
inl (state_1d (t1, _, _)) = sts 1
inl dt = t1 - t0
inl num_steps = t / dt |> math.round |> abs
inl (state_1d (_, _, v0)) = sts num_steps
v0
inl velocity_ftxv dt m txv0 fs =
states_txv dt m txv0 fs |> velocity_1d
inl position_1d sts t =
inl (state_1d (t0, _, _)) = sts 0
inl (state_1d (t1, _, _)) = sts 1
inl dt = t1 - t0
inl num_steps = t / dt |> math.round |> abs
inl (state_1d (_, x0, _)) = sts num_steps
x0
inl position_ftxv dt m txv0 fs =
states_txv dt m txv0 fs |> position_1d
inl spring_force k (state_1d (_, x0, _)) =
-k * x0
InΒ [Β ]:
//// test
inl damped_ho_forces () =
[
spring_force 0.8
fun (state_1d (_, _, v0)) => f_air 2 1.225 (pi * math.square 0.02) v0
fun _ => -0.0027 * 9.80665
]
inl damped_ho_states () =
states_txv 0.001 0.0027 (state_1d (0, 0.1, 0)) (damped_ho_forces ())
inl pingpong_position t =
position_ftxv 0.001 0.0027 (state_1d (0, 0.1, 0)) (damped_ho_forces ()) t
inl x = am'.init_series 0 3 0.01
inl y = x |> am'.map_base pingpong_position
"ping pong ball on a slinky", "time (s)", "", ;[ "position (m)", x, y ]
InΒ [Β ]:
//// test
inl pingpong_velocity t =
velocity_ftxv 0.001 0.0027 (state_1d (0, 0.1, 0)) (damped_ho_forces ()) t
inl x = am'.init_series 0 3 0.01
inl y = x |> am'.map_base pingpong_velocity
"ping pong ball on a slinky", "time (s)", "", ;[ "velocity (m/s)", x, y ]
shiftΒΆ
InΒ [Β ]:
type update_function s = s -> s
type differential_equation s ds = s -> ds
type numerical_method s ds = differential_equation s ds -> update_function s
inl solver method =
method >> seq.iterate
inl solver' method =
method >> seq.iterate'
inl solver_ method =
method >> seq.iterate_
inl euler_cromer_1d dt deriv (state_1d (t0, x0, v0) as t) =
inl (rrr (_, _, dvdt)) = deriv t
inl t1 = t0 + dt
inl v1 = v0 + dvdt * dt
inl x1 = x0 + v1 * dt
state_1d (t1, x1, v1)
inl update_txv_ec dt m fs =
euler_cromer_1d dt (newton_second_1d m fs)
prototype (+++) ds : ds -> ds -> ds
prototype scale ds : f64 -> ds -> ds
instance (+++) rrr = fun (rrr (dtdt0, dxdt0, dvdt0)) (rrr (dtdt1, dxdt1, dvdt1)) =>
rrr (dtdt0 + dtdt1, dxdt0 + dxdt1, dvdt0 + dvdt1)
instance scale rrr = fun w (rrr (dtdt0, dxdt0, dvdt0)) =>
rrr (w * dtdt0, w * dxdt0, w * dvdt0)
prototype shift s : forall ds. f64 -> ds -> s -> s
instance shift state_1d = fun dt ds (state_1d (t, x, v)) =>
inl dtdt, dxdt, dvdt =
real
match ds with
| rrr x => x
| state_1d x => x
state_1d (t + dtdt * dt, x + dxdt * dt, v + dvdt * dt)
inl euler dt deriv st0 =
shift dt (deriv st0) st0
inl runge_kutta_4 dt deriv st0 =
inl m0 = deriv st0
inl m1 = deriv (shift (dt / 2) m0 st0)
inl m2 = deriv (shift (dt / 2) m1 st0)
inl m3 = deriv (shift dt m2 st0)
shift (dt / 6) (m0 +++ m1 +++ m1 +++ m2 +++ m2 +++ m3) st0
inl exponential (_, x0, v0) =
1f64, v0, x0
inl of_state_1d (state_1d (t, x, v)) =
t, x, v
InΒ [Β ]:
//// test
solver (euler 0.01) (of_state_1d >> exponential >> state_1d) (state_1d (0, 1, 1)) 800i32
|> _assert_eq (state_1d (7.999999999999874, 2864.8311229272326, 2864.8311229272326))
solver (euler_cromer_1d 0.1) (of_state_1d >> exponential >> rrr) (state_1d (0, 1, 1)) 80i32
|> _assert_eq (state_1d (7.999999999999988, 3043.379244966009, 2895.0121485099035))
solver (runge_kutta_4 1) (of_state_1d >> exponential >> rrr) (state_1d (0, 1, 1)) 8i32
|> _assert_eq (state_1d (8.0, 2894.789038540849, 2894.789038540849))
{ name = __assert_eq; actual = struct (8.0, 2864.831123, 2864.831123); expected = struct (8.0, 2864.831123, 2864.831123) }
{ name = __assert_eq; actual = struct (8.0, 3043.379245, 2895.012149); expected = struct (8.0, 3043.379245, 2895.012149) }
{ name = __assert_eq; actual = struct (8.0, 2894.789039, 2894.789039); expected = struct (8.0, 2894.789039, 2894.789039) }
vecΒΆ
InΒ [Β ]:
type vec =
{
x : f64
y : f64
z : f64
}
inl vec x y z : vec =
{ x y z }
InΒ [Β ]:
//// test
vec 1 2 3 .z
|> _assert_eq 3
{ name = __assert_eq; actual = +3.000000; expected = +3.000000 }
constsΒΆ
InΒ [Β ]:
inl i_hat () = vec 1 0 0
inl j_hat () = vec 0 1 0
inl k_hat () = vec 0 0 1
inl zero_vec () = vec 0 0 0
^+^ΒΆ
InΒ [Β ]:
inl (^+^) (a : vec) (b : vec) =
vec (a.x + b.x) (a.y + b.y) (a.z + b.z)
InΒ [Β ]:
//// test
vec 1 2 3 ^+^ vec 4 5 6
|> _assert_eq (vec 5 7 9)
{ name = __assert_eq; actual = { x = +5.000000; y = +7.000000; z = +9.000000 }; expected = { x = +5.000000; y = +7.000000; z = +9.000000 } }
sum_vecΒΆ
InΒ [Β ]:
inl sum_vec vs =
vs |> listm.fold (^+^) (zero_vec ())
InΒ [Β ]:
//// test
[ vec 1 2 3; vec 4 5 6 ]
|> sum_vec
|> _assert_eq (vec 5 7 9)
{ name = __assert_eq; actual = { x = +5.000000; y = +7.000000; z = +9.000000 }; expected = { x = +5.000000; y = +7.000000; z = +9.000000 } }
*^ΒΆ
InΒ [Β ]:
inl (*^) c { x y z } =
vec (c * x) (c * y) (c * z)
InΒ [Β ]:
//// test
5 *^ vec 1 2 3
|> _assert_eq (vec 5 10 15)
{ name = __assert_eq; actual = { x = +5.000000; y = +10.000000; z = +15.000000 }; expected = { x = +5.000000; y = +10.000000; z = +15.000000 } }
InΒ [Β ]:
//// test
3 *^ i_hat () ^+^ 4 *^ k_hat ()
|> _assert_eq (vec 3 0 4)
{ name = __assert_eq; actual = { x = +3.000000; y = +0.000000; z = +4.000000 }; expected = { x = +3.000000; y = +0.000000; z = +4.000000 } }
^*ΒΆ
InΒ [Β ]:
inl (^*) v c =
(*^) c v
InΒ [Β ]:
//// test
vec 1 2 3 ^* 5
|> _assert_eq (vec 5 10 15)
{ name = __assert_eq; actual = { x = +5.000000; y = +10.000000; z = +15.000000 }; expected = { x = +5.000000; y = +10.000000; z = +15.000000 } }
^/ΒΆ
InΒ [Β ]:
inl (^/) { x y z } c =
vec (x / c) (y / c) (z / c)
InΒ [Β ]:
//// test
vec 1 2 3 ^/ 5
|> _assert_eq (vec 0.2 0.4 0.6)
{ name = __assert_eq; actual = { x = +0.200000; y = +0.400000; z = +0.600000 }; expected = { x = +0.200000; y = +0.400000; z = +0.600000 } }
negate_vecΒΆ
InΒ [Β ]:
inl negate_vec v =
v ^* -1
InΒ [Β ]:
//// test
vec 1 2 3
|> negate_vec
|> _assert_eq (vec -1 -2 -3)
{ name = __assert_eq; actual = { x = -1.000000; y = -2.000000; z = -3.000000 }; expected = { x = -1.000000; y = -2.000000; z = -3.000000 } }
^-^ΒΆ
InΒ [Β ]:
inl (^-^) a b =
a ^+^ (negate_vec b)
InΒ [Β ]:
//// test
vec 1 2 3 ^-^ vec 4 5 6
|> _assert_eq (vec -3 -3 -3)
{ name = __assert_eq; actual = { x = -3.000000; y = -3.000000; z = -3.000000 }; expected = { x = -3.000000; y = -3.000000; z = -3.000000 } }
<.>ΒΆ
InΒ [Β ]:
inl (<.>) { x = ax y = ay z = az } { x = bx y = by z = bz } =
ax * bx + ay * by + az * bz
InΒ [Β ]:
//// test
vec 1 2 3 <.> vec 4 5 6
|> _assert_eq 32
{ name = __assert_eq; actual = +32.000000; expected = +32.000000 }
><ΒΆ
InΒ [Β ]:
inl (><) (a : vec) (b : vec) =
vec
(a.y * b.z - a.z * b.y)
(a.z * b.x - a.x * b.z)
(a.x * b.y - a.y * b.x)
InΒ [Β ]:
//// test
vec 1 2 3 >< vec 4 5 6
|> _assert_eq (vec -3 6 -3)
{ name = __assert_eq; actual = { x = -3.000000; y = +6.000000; z = -3.000000 }; expected = { x = -3.000000; y = +6.000000; z = -3.000000 } }
magnitudeΒΆ
InΒ [Β ]:
inl magnitude v =
v <.> v |> sqrt
InΒ [Β ]:
//// test
vec 1 2 3
|> magnitude
|> _assert_approx_eq None 3.7416573867739413
{ name = __assert_approx_eq; actual = +3.741657; expected = +3.741657 }
v1ΒΆ
InΒ [Β ]:
inl v1 t =
2 *^ (t ** 2 *^ i_hat () ^+^ 3 *^ (t ** 3 *^ j_hat () ^+^ t ** 4 *^ k_hat ()))
InΒ [Β ]:
//// test
v1 1
|> _assert_eq (vec 2 6 6)
{ name = __assert_eq; actual = { x = +2.000000; y = +6.000000; z = +6.000000 }; expected = { x = +2.000000; y = +6.000000; z = +6.000000 } }
vec_derivativeΒΆ
InΒ [Β ]:
type vec_derivative = (f64 -> vec) -> f64 -> vec
inl vec_derivative dt : vec_derivative =
fun v t =>
(v (t + dt / 2) ^-^ v (t - dt / 2)) ^/ dt
InΒ [Β ]:
//// test
vec_derivative 0.01 v1 3 .x
|> _assert_approx_eq None (derivative 0.01 (v1 >> fun v => v.x) 3)
{ name = __assert_approx_eq; actual = +12.000000; expected = +12.000000 }
states_psΒΆ
InΒ [Β ]:
nominal particle_state =
{
mass : f64
charge : f64
time : f64
pos_vec : vec
velocity : vec
}
inl default_particle_state () : particle_state =
particle_state {
mass = 1
charge = 0
time = 0
pos_vec = zero_vec ()
velocity = zero_vec ()
}
type one_body_force = particle_state -> vec
nominal d_particle_state =
{
dmdt : f64
dqdt : f64
dtdt : f64
drdt : vec
dvdt : vec
}
inl newton_second_ps (fs : list one_body_force) (st : particle_state) : d_particle_state =
inl f_net = fs |> listm.map (fun f => f st) |> sum_vec
d_particle_state {
dmdt = 0
dqdt = 0
dtdt = 1
drdt = st.velocity
dvdt = f_net ^/ st.mass
}
inl earth_surface_gravity (st : particle_state) =
inl g = 9.80665
-st.mass * g *^ k_hat ()
inl air_resistance drag rho area (st : particle_state) =
-0.5 * drag * rho * area * magnitude st.velocity *^ st.velocity
inl euler_cromer_ps dt (deriv : particle_state -> d_particle_state) (particle_state st) =
inl dst : d_particle_state = deriv (particle_state st)
inl v' = st.velocity ^+^ dst.dvdt ^* dt
particle_state { st with
time = st.time + dt
pos_vec = st.pos_vec ^+^ v' ^* dt
velocity = st.velocity ^+^ dst.dvdt ^* dt
}
instance (+++) d_particle_state = fun (dps : d_particle_state) (dps' : d_particle_state) =>
d_particle_state {
dmdt = dps.dmdt + dps'.dmdt
dqdt = dps.dqdt + dps'.dqdt
dtdt = dps.dtdt + dps'.dtdt
drdt = dps.drdt ^+^ dps'.drdt
dvdt = dps.dvdt ^+^ dps'.dvdt
}
instance scale d_particle_state = fun w (dps : d_particle_state) =>
d_particle_state {
dmdt = w * dps.dmdt
dqdt = w * dps.dqdt
dtdt = w * dps.dtdt
drdt = w *^ dps.drdt
dvdt = w *^ dps.dvdt
}
instance shift particle_state = fun dt dps (particle_state st) =>
inl (d_particle_state dps) =
real
match dps with
| d_particle_state _ => dps
particle_state { st with
time = st.time + dps.dtdt * dt
pos_vec = st.pos_vec ^+^ dps.drdt ^* dt
velocity = st.velocity ^+^ dps.dvdt ^* dt
}
inl states_ps (method : numerical_method particle_state d_particle_state) : _ -> _ -> i32 -> particle_state =
newton_second_ps >> method >> seq.iterate_
inl z_ge0 sts =
sts
|> seq.take_while_ (fun (particle_state st) _ => st.pos_vec.z >= 0)
inl trajectory sts =
sts |> listm.map (fun (particle_state st) => st.pos_vec.y, st.pos_vec.z)
InΒ [Β ]:
//// test
inl update_ps (method : numerical_method particle_state d_particle_state) =
newton_second_ps >> method
inl position_ps (method : numerical_method particle_state d_particle_state) fs st t =
inl states : i32 -> particle_state = states_ps method fs st
inl dt = (states 1).time - (states 0).time
inl num_steps = t / dt |> math.round |> abs
inl st1 = solver' method (newton_second_ps fs) st num_steps
st1.pos_vec
inl sun_gravity (st : particle_state) : vec =
inl big_g = 0.0000000000667408
inl sun_mass = 1988480000000000000000000000000
-big_g * sun_mass * st.mass *^ st.pos_vec ^/ magnitude st.pos_vec ** 3
inl wind_force v_wind drag rho area (st : particle_state) =
inl v_rel = st.velocity ^-^ v_wind
-0.5 * drag * rho * area * magnitude v_rel *^ v_rel
inl rock_state () =
inl (particle_state default_particle_state') = default_particle_state ()
particle_state { default_particle_state' with
mass = 2
velocity = vec 3 0 4
}
inl halley_update dt =
update_ps (euler_cromer_ps dt) [ sun_gravity ]
inl halley_initial () =
inl (particle_state default_particle_state') = default_particle_state ()
particle_state { default_particle_state' with
mass = 220000000000000
pos_vec = 87660000000 *^ i_hat ()
velocity = 54569 *^ j_hat ()
}
InΒ [Β ]:
//// test
inl baseball_forces () =
inl area = pi * (0.074 / 2) ** 2
[
earth_surface_gravity
air_resistance 0.3 1.225 area
]
inl baseball_trajectory dt v0 theta_deg =
inl theta_rad = theta_deg * pi / 180
inl vy0 = v0 * cos theta_rad
inl vz0 = v0 * sin theta_rad
inl initial_state =
particle_state {
mass = 0.145
charge = 0
time = 0
pos_vec = zero_vec ()
velocity = vec 0 vy0 vz0
}
states_ps (euler_cromer_ps dt) (baseball_forces ()) initial_state
>> Some
|> z_ge0
|> trajectory
inl baseball_range dt v0 theta_deg =
baseball_trajectory dt v0 theta_deg
|> listm.fold (fun _ (y, _) => y) 0
inl x = am'.init_series 10 80 1
inl y = x |> am'.map_base (baseball_range 0.01 45)
"range for a baseball hit at 45 m/s",
"angle above horizontal (degrees)",
"",
;[ "horizontal range (m)", x, y ]
InΒ [Β ]:
//// test
inl best_angle (min, max) =
let rec 루ν theta_deg (best_range, best_theta_deg) =
if theta_deg > max
then best_range, best_theta_deg
else
inl range = baseball_range 0.01 45 theta_deg
루ν
(theta_deg + 1)
(if range > best_range
then range, theta_deg
else best_range, best_theta_deg)
루ν min (0f64, min)
best_angle (30f64, 60f64)
|> _assert_eq (116.77499158246208, 41)
{ name = __assert_eq; actual = +116.774992, +41.000000; expected = +116.774992, +41.000000 }
relativity_psΒΆ
InΒ [Β ]:
inl relativity_ps fs (st : particle_state) =
inl f_net = fs |> listm.map (fun f => f st) |> sum_vec
inl c = 299792458
inl u = st.velocity ^/ c
inl acc = sqrt (1 - (u <.> u)) *^ (f_net ^-^ (f_net <.> u) *^ u) ^/ st.mass
d_particle_state {
dmdt = 0
dqdt = 0
dtdt = 1
drdt = st.velocity
dvdt = acc
}
InΒ [Β ]:
//// test
inl year = 365.25 * 24 * 60 * 60
inl c = 299792458
inl ~method = runge_kutta_4 100000
inl forces = [ fun _ => 10 *^ i_hat () ]
inl (particle_state default_particle_state') = default_particle_state ()
inl initial_state =
particle_state { default_particle_state' with
mass = 1
}
inl newton_states = solver_ method (newton_second_ps forces) initial_state
inl relativity_states = solver_ method (relativity_ps forces) initial_state
inl newton_x, newton_y =
newton_states
>> Some
|> seq.take_while_ (fun (particle_state st) (_ : i32) => st.time <= year)
|> listm.map (fun (particle_state st) => st.time / year, st.velocity.x / c)
|> listm'.unzip
inl _, relativity_y =
relativity_states
>> Some
|> seq.take_while_ (fun (particle_state st) (_ : i32) => st.time <= year)
|> listm.map (fun (particle_state st) => st.time / year, st.velocity.x / c)
|> listm'.unzip
inl newton_x = newton_x |> listm'.box |> listm'.to_array'
inl newton_y = newton_y |> listm'.box |> listm'.to_array'
inl relativity_y = relativity_y |> listm'.box |> listm'.to_array'
"response to a constant force",
"time (years)",
"velocity (multiples of c)",
;[
"newtonian", newton_x, newton_y
"relativistic", newton_x, relativity_y
]
InΒ [Β ]:
inl uniform_lorentz_force v_e v_b (st : particle_state) =
st.charge *^ (v_e ^+^ st.velocity >< v_b)
InΒ [Β ]:
//// test
inl c : f64 = 299792458
inl ~method = runge_kutta_4 0.000000001
inl forces = [ uniform_lorentz_force (zero_vec ()) (k_hat ()) ]
inl (particle_state default_particle_state') = default_particle_state ()
inl initial_state =
particle_state { default_particle_state' with
mass = 0.000000000000000000000000001672621898
charge = 0.0000000000000000001602176621
velocity = 0.8 *^ (c *^ j_hat ())
}
inl newton_states = solver_ method (newton_second_ps forces) initial_state
inl relativity_states = solver_ method (relativity_ps forces) initial_state
inl newton_x, newton_y =
newton_states
>> Some
|> seq.take_while_ (fun (particle_state st) i => i < 100i32)
|> listm.map (fun (particle_state st) => st.pos_vec.x, st.pos_vec.y)
|> listm'.unzip
inl relativity_x, relativity_y =
relativity_states
>> Some
|> seq.take_while_ (fun (particle_state st) i => i < 165i32)
|> listm.map (fun (particle_state st) => st.pos_vec.x, st.pos_vec.y)
|> listm'.unzip
inl newton_x = newton_x |> listm'.box |> listm'.to_array'
inl newton_y = newton_y |> listm'.box |> listm'.to_array'
inl relativity_x = relativity_x |> listm'.box |> listm'.to_array'
inl relativity_y = relativity_y |> listm'.box |> listm'.to_array'
"proton in a 1-t magnetic field",
"x (m)",
"y (m)",
;[
"newtonian", newton_x, newton_y
"relativistic", relativity_x, relativity_y
]
system kinetic energy versus time 1ΒΆ
InΒ [Β ]:
//// test
inl central_force f (particle_state st1) (particle_state st2) =
inl r1 = st1.pos_vec
inl r2 = st2.pos_vec
inl r21 = r2 ^-^ r1
inl r21mag = magnitude r21
f r21mag *^ r21 ^/ r21mag
inl billiard_force k re =
inl f r =
if r >= re
then 0
else -k * (r - re)
central_force f
type force_vector = vec
type two_body_force = particle_state -> particle_state -> force_vector
union force =
| ExternalForce : i32 * one_body_force
| InternalForce : i32 * i32 * two_body_force
nominal multi_particle_state = list particle_state
nominal d_multi_particle_state = list d_particle_state
inl force_on n sts force =
match force with
| ExternalForce (n0, f_one_body) =>
if n = n0
then f_one_body
else fun _ => zero_vec ()
| InternalForce (n0, n1, f_two_body) =>
if n = n0
then f_two_body (sts |> listm'.item n1)
elif n = n1
then f_two_body (sts |> listm'.item n0)
else fun _ => zero_vec ()
inl forces_on n (multi_particle_state sts) fs =
fs |> listm.map (force_on n sts)
inl newton_second_mps fs (multi_particle_state sts) : d_multi_particle_state =
inl deriv (n, st) =
newton_second_ps (forces_on n (multi_particle_state sts) fs) st
sts |> listm'.indexed |> listm.map deriv |> d_multi_particle_state
instance (+++) d_multi_particle_state = fun (d_multi_particle_state dsts1) (d_multi_particle_state dsts2) =>
d_multi_particle_state (listm'.zip_with_ (+++) dsts1 dsts2)
instance scale d_multi_particle_state = fun w (d_multi_particle_state dsts) =>
d_multi_particle_state (dsts |> listm.map (scale w))
instance shift multi_particle_state = fun dt dsts (multi_particle_state sts) =>
inl (d_multi_particle_state dsts) =
real
match dsts with
| d_multi_particle_state _ => dsts
listm'.zip_with_ (shift dt) dsts sts |> multi_particle_state
inl euler_cromer_mps dt : numerical_method multi_particle_state d_multi_particle_state =
fun deriv mpst0 =>
inl mpst1 = euler dt deriv mpst0
inl (multi_particle_state sts0) = mpst0
inl (multi_particle_state sts1) = mpst1
sts1
|> listm'.zip_ sts0
|> listm.map (fun ((particle_state st0), (particle_state st1)) =>
particle_state {
st1 with
pos_vec = st0.pos_vec ^+^ st1.velocity ^* dt
}
)
|> multi_particle_state
inl update_mps (method : numerical_method multi_particle_state d_multi_particle_state) =
newton_second_mps >> method
inl states_mps (method : numerical_method multi_particle_state d_multi_particle_state) =
newton_second_mps >> method >> seq.iterate_
inl kinetic_energy (particle_state st) =
inl m = st.mass
inl v = magnitude st.velocity
0.5 * m * v ** 2
inl system_ke (multi_particle_state sts) =
sts |> listm.map kinetic_energy |> listm'.sum
inl linear_spring_pe k re (particle_state st1) (particle_state st2) =
inl r1 = st1.pos_vec
inl r2 = st2.pos_vec
inl r21 = r2 ^-^ r1
inl r21mag = magnitude r21
k * (r21mag - re) ** 2 / 2
inl earth_surface_gravity_pe (particle_state st) =
inl g = 9.80665
inl m = st.mass
inl z = st.pos_vec.z
m * g * z
inl two_springs_pe (multi_particle_state sts) =
inl st0 = sts |> listm'.item 0i32
inl st1 = sts |> listm'.item 1i32
linear_spring_pe 100 0.5 (default_particle_state ()) st0
+ linear_spring_pe 100 0.5 st0 st1
+ earth_surface_gravity_pe st0
+ earth_surface_gravity_pe st1
inl two_springs_me mpst =
system_ke mpst + two_springs_pe mpst
inl ball_radius () = 0.03
inl billiard_forces k =
[ InternalForce (0, 1, billiard_force k (2 * ball_radius ())) ]
inl billiard_update n_method k dt =
update_mps (n_method dt) (billiard_forces k)
inl billiard_initial () =
inl ball_mass = 0.160
inl (particle_state default_particle_state') = default_particle_state ()
multi_particle_state [
particle_state {
default_particle_state' with
mass = ball_mass
pos_vec = zero_vec ()
velocity = 0.2 *^ i_hat ()
}
particle_state {
default_particle_state' with
mass = ball_mass
pos_vec = i_hat () ^+^ 0.02 *^ j_hat ()
velocity = zero_vec ()
}
]
inl billiard_states ~n_method k dt =
states_mps (n_method dt) (billiard_forces k) (billiard_initial ())
inl billiard_states_finite n_method k dt =
billiard_states n_method k dt
>> Some
|> seq.take_while_ (fun (multi_particle_state mpst) (_ : i32) =>
(mpst |> listm'.item 0i32).time <= 10
)
inl momentum (particle_state st) =
inl m = st.mass
inl v = st.velocity
m *^ v
inl system_p (multi_particle_state sts) =
sts |> listm.map momentum |> sum_vec
inl time_ke_ec_x, time_ke_ec_y =
billiard_states_finite euler_cromer_mps 30 0.03
|> listm.map (fun (multi_particle_state mpst) =>
(mpst |> listm'.item 0i32).time, system_ke (multi_particle_state mpst)
)
|> listm'.unzip
inl time_ke_rk4_x, time_ke_rk4_y =
billiard_states_finite runge_kutta_4 30 0.03
|> listm.map (fun (multi_particle_state mpst) =>
(mpst |> listm'.item 0i32).time, system_ke (multi_particle_state mpst)
)
|> listm'.unzip
inl time_ke_ec_x = time_ke_ec_x |> listm'.box |> listm'.to_array'
inl time_ke_ec_y = time_ke_ec_y |> listm'.box |> listm'.to_array'
inl time_ke_rk4_x = time_ke_rk4_x |> listm'.box |> listm'.to_array'
inl time_ke_rk4_y = time_ke_rk4_y |> listm'.box |> listm'.to_array'
"system kinetic energy versus time",
"time (s)",
"system kinetic energy (j)",
;[
"euler-cromer", time_ke_ec_x, time_ke_ec_y
"runge-kutta 4", time_ke_rk4_x, time_ke_rk4_y
]
wave 1ΒΆ
InΒ [Β ]:
//// test
inl linear_spring k re (particle_state st1) (particle_state st2) =
inl r1 = st1.pos_vec
inl r2 = st2.pos_vec
inl r21 = r2 ^-^ r1
inl r21mag = magnitude r21
-k * (r21mag - re) *^ r21 ^/ r21mag
inl fixed_linear_spring k re r1 =
inl (particle_state default_particle_state') = default_particle_state ()
linear_spring k re (particle_state { default_particle_state' with pos_vec = r1 })
inl forces_string () =
[
ExternalForce (0, fixed_linear_spring 5384 0 (zero_vec ()))
ExternalForce (63, fixed_linear_spring 5384 0 (0.65 *^ i_hat ()))
] ++ (
listm'.init_series 0 59 1
|> listm.map (fun n => InternalForce (n, n + 1, linear_spring 5384 0))
)
inl string_update dt =
update_mps (runge_kutta_4 dt) (forces_string ())
inl string_initial_overtone n =
inl ball_mass = 0.0008293 * 0.65 / 64
inl (particle_state default_particle_state') = default_particle_state ()
listm'.init_series 0.01 0.64 0.01
|> listm.map (fun x =>
inl y = 0.005 * sin (conv n * pi * x / 0.65)
particle_state {
default_particle_state' with
mass = ball_mass
pos_vec = x *^ i_hat () ^+^ y *^ j_hat ()
velocity = zero_vec ()
}
)
|> multi_particle_state
inl string_initial_pluck () =
inl ball_mass = 0.0008293 * 0.65 / 64
inl (particle_state default_particle_state') = default_particle_state ()
listm'.init_series 0.01 0.64 0.01
|> listm.map (fun x =>
inl y =
inl n = if x <= 0.51 then 0 else 0.65
0.005 / (0.51 - n) * (x - n)
particle_state {
default_particle_state' with
mass = ball_mass
pos_vec = x *^ i_hat () ^+^ y *^ j_hat ()
velocity = zero_vec ()
}
)
|> multi_particle_state
let main () =
inl ~frames = listm'.init_series 0 9 1f64
inl initial_state = string_initial_overtone 3i32
inl frames =
frames
|> listm.map (fun n =>
inl (multi_particle_state sts) =
seq.iterate' (string_update 0.000025) initial_state |> fun f => f 0f64
inl rs =
[ zero_vec () ]
++ (sts |> listm.map (fun (particle_state st) => st.pos_vec))
++ [ 0.65 *^ i_hat () ]
inl x, y =
rs
|> listm.map (fun r => r.x, r.y)
|> listm'.unzip
inl x = x |> listm'.box |> listm'.to_array'
inl y = y |> listm'.box |> listm'.to_array'
x, y
)
|> listm'.box |> listm'.to_array'
inl n = 0i32
inl x, y = a frames |> am'.index n
"wave",
"position (m)",
"displacement (m)",
;[
($'$"{!n}"' : string), x, y
]
system kinetic energy versus time 2ΒΆ
InΒ [Β ]:
//// test
inl central_force f (particle_state st1) (particle_state st2) =
inl r1 = st1.pos_vec
inl r2 = st2.pos_vec
inl r21 = r2 ^-^ r1
inl r21mag = magnitude r21
f r21mag *^ r21 ^/ r21mag
inl billiard_force k re =
inl f r =
if r >= re
then 0
else -k * (r - re)
central_force f
type force_vector = vec
type two_body_force = particle_state -> particle_state -> force_vector
union force t =
| ExternalForce : t * one_body_force
| InternalForce : t * t * two_body_force
nominal multi_particle_state = stream.stream particle_state
nominal d_multi_particle_state = stream.stream d_particle_state
inl force_on n s force =
match force with
| ExternalForce (n0, f_one_body) =>
if n = n0
then f_one_body
else fun _ => zero_vec ()
| InternalForce (n0, n1, f_two_body) =>
if n = n0
then s |> stream.try_item n1 |> optionm.map f_two_body
elif n = n1
then s |> stream.try_item n0 |> optionm.map f_two_body
else None
|> optionm'.default_value (fun _ => zero_vec ())
inl forces_on n (multi_particle_state sts) fs =
fs
|> listm.map (force_on n sts)
inl newton_second_mps fs ((multi_particle_state sts) as mpst) =
inl deriv (n, st) =
newton_second_ps (forces_on n mpst fs) st
sts |> stream.indexed |> stream.map deriv |> d_multi_particle_state
instance (+++) d_multi_particle_state =
fun (d_multi_particle_state dsts1) (d_multi_particle_state dsts2) =>
(dsts1, dsts2)
||> stream.zip_with (+++)
|> d_multi_particle_state
instance scale d_multi_particle_state = fun w (d_multi_particle_state dsts) =>
dsts
|> stream.map (scale w)
|> d_multi_particle_state
instance shift multi_particle_state = fun dt dsts (multi_particle_state sts) =>
inl (d_multi_particle_state dsts) =
real
match dsts with
| d_multi_particle_state _ => dsts
(dsts, sts)
||> stream.zip_with (shift dt)
|> stream.memoize
|> fun x => x ()
|> multi_particle_state
inl euler_cromer_mps dt : numerical_method multi_particle_state d_multi_particle_state =
fun deriv ((multi_particle_state sts0) as mpst0) =>
inl (multi_particle_state sts1) = euler dt deriv mpst0
(sts0, sts1)
||> stream.zip
|> stream.map (fun ((particle_state st0), (particle_state st1)) =>
particle_state {
st1 with
pos_vec = st0.pos_vec ^+^ st1.velocity ^* dt
}
)
|> multi_particle_state
inl update_mps (method : numerical_method multi_particle_state d_multi_particle_state) =
newton_second_mps >> method
inl states_mps (method : numerical_method multi_particle_state d_multi_particle_state) =
newton_second_mps
>> method
>> (fun x (multi_particle_state y) =>
y
|> stream.memoize
|> (fun x => x ())
|> multi_particle_state |> x
)
// >> stream.iterate
>> seq.iterate'
inl kinetic_energy (particle_state st) =
inl m = st.mass
inl v = magnitude st.velocity
0.5 * m * v ** 2
inl system_ke (multi_particle_state sts) =
sts
|> stream.map kinetic_energy
|> stream.sum
inl linear_spring_pe k re (particle_state st1) (particle_state st2) =
inl r1 = st1.pos_vec
inl r2 = st2.pos_vec
inl r21 = r2 ^-^ r1
inl r21mag = magnitude r21
k * (r21mag - re) ** 2 / 2
inl earth_surface_gravity_pe (particle_state st) =
inl g = 9.80665
inl m = st.mass
inl z = st.pos_vec.z
m * g * z
inl ball_radius () = 0.03
inl billiard_forces k =
[ InternalForce (0i32, 1, billiard_force k (2 * ball_radius ())) ]
inl billiard_initial () =
inl ball_mass = 0.160
inl (particle_state default_particle_state') = default_particle_state ()
[
particle_state {
default_particle_state' with
mass = ball_mass
pos_vec = zero_vec ()
velocity = 0.2 *^ i_hat ()
}
particle_state {
default_particle_state' with
mass = ball_mass
pos_vec = i_hat () ^+^ 0.02 *^ j_hat ()
velocity = zero_vec ()
}
]
|> stream.from_list
|> multi_particle_state
inl billiard_states ~n_method k dt =
states_mps (n_method dt) (billiard_forces k) (billiard_initial ())
inl billiard_states_finite n_method k dt =
billiard_states n_method k dt
>> Some
|> seq.take_while_ (fun (multi_particle_state mpst) (_ : i32) =>
match mpst |> stream.try_item 0i32 with
| Some st =>
st.time <= 10
| None => false
)
inl momentum (particle_state st) =
inl m = st.mass
inl v = st.velocity
m *^ v
inl system_p (multi_particle_state sts) =
sts
|> stream.map momentum
|> stream.fold (^+^) (zero_vec ())
inl time_ke_ec_x, time_ke_ec_y =
billiard_states_finite euler_cromer_mps 30 0.03
|> listm.map (fun (multi_particle_state mpst) =>
mpst |> stream.try_item 0i32
|> optionm.map (fun st =>
st.time, system_ke (multi_particle_state mpst)
)
)
// |> stream.to_list
|> listm'.choose id
|> listm'.unzip
inl time_ke_rk4_x, time_ke_rk4_y =
billiard_states_finite runge_kutta_4 30 0.03
|> listm.map (fun (multi_particle_state mpst) =>
mpst |> stream.try_item 0i32
|> optionm.map (fun st =>
st.time, system_ke (multi_particle_state mpst)
)
)
// |> stream.to_list
|> listm'.choose id
|> listm'.unzip
inl time_ke_ec_x = time_ke_ec_x |> listm'.box |> listm'.to_array'
inl time_ke_ec_y = time_ke_ec_y |> listm'.box |> listm'.to_array'
inl time_ke_rk4_x = time_ke_rk4_x |> listm'.box |> listm'.to_array'
inl time_ke_rk4_y = time_ke_rk4_y |> listm'.box |> listm'.to_array'
"system kinetic energy versus time",
"time (s)",
"system kinetic energy (j)",
;[
"euler-cromer", time_ke_ec_x, time_ke_ec_y
"runge-kutta 4", time_ke_rk4_x, time_ke_rk4_y
]
wave 2ΒΆ
InΒ [Β ]:
//// test
inl linear_spring k re (particle_state st1) (particle_state st2) =
inl r1 = st1.pos_vec
inl r2 = st2.pos_vec
inl r21 = r2 ^-^ r1
inl r21mag = magnitude r21
-k * (r21mag - re) *^ r21 ^/ r21mag
inl fixed_linear_spring k re r1 =
inl (particle_state default_particle_state') = default_particle_state ()
linear_spring k re (particle_state { default_particle_state' with pos_vec = r1 })
inl forces_string () =
[
ExternalForce (0i32, fixed_linear_spring 5384 0 (zero_vec ()))
ExternalForce (63, fixed_linear_spring 5384 0 (0.65 *^ i_hat ()))
] ++ (
listm'.init_series 0 59 1
|> listm.map (fun n => InternalForce (n, n + 1, linear_spring 5384 0))
)
inl string_update dt =
update_mps (join runge_kutta_4 dt) (join forces_string ())
inl string_initial_overtone n =
inl ball_mass = 0.0008293 * 0.65 / 64
inl (particle_state default_particle_state') = default_particle_state ()
listm'.init_series 0.01 0.64 0.01
|> listm.map (fun x =>
inl y = 0.005 * sin (conv n * pi * x / 0.65)
particle_state {
default_particle_state' with
mass = ball_mass
pos_vec = x *^ i_hat () ^+^ y *^ j_hat ()
velocity = zero_vec ()
}
)
|> stream.from_list
|> multi_particle_state
let main () =
inl ~frames = listm'.init_series 0 65 1f64 |> stream.from_list
inl ~initial_state = string_initial_overtone 3i32
inl frames =
frames
|> stream.map (fun n =>
inl (multi_particle_state sts) =
stream.iterate (string_update 0.000025) initial_state |> stream.item n
inl x, y =
[ zero_vec () ]
++ (sts |> stream.map (fun (particle_state st) => st.pos_vec) |> stream.to_list)
++ [ 0.65 *^ i_hat () ]
|> listm.map (fun r => r.x, r.y)
|> stream.from_list
|> stream.unzip
inl x = x |> stream.to_list |> listm'.box |> listm'.to_array'
inl y = y |> stream.to_list |> listm'.box |> listm'.to_array'
x, y
)
inl plots =
frames
|> stream.indexed
|> stream.map (fun ((n : i32), (x, y)) =>
"wave",
"position (m)",
"displacement (m)",
;[
($'$"{!n}"' : string), x, y
]
)
plots |> stream.to_list |> listm'.box |> listm'.to_array'
| index | value |
|---|---|
| 0 | |
| 1 | |
| 2 | |
| 3 | |
| 4 | |
| 5 | |
| 6 | |
| 7 | |
| 8 | |
| 9 | |
| 10 | |
| 11 | |
| 12 | |
| 13 | |
| 14 | |
| 15 | |
| 16 | |
| 17 | |
| 18 | |
| 19 | |
| 20 | |
| 21 | |
| 22 | |
| 23 | |
| 24 | |
| 25 | |
| 26 | |
| 27 | |
| 28 | |
| 29 | |
| 30 | |
| 31 | |
| 32 | |
| 33 | |
| 34 | |
| 35 | |
| 36 | |
| 37 | |
| 38 | |
| 39 | |
| 40 | |
| 41 | |
| 42 | |
| 43 | |
| 44 | |
| 45 | |
| 46 | |
| 47 | |
| 48 | |
| 49 | |
| 50 | |
| 51 | |
| 52 | |
| 53 | |
| 54 | |
| 55 | |
| 56 | |
| 57 | |
| 58 | |
| 59 | |
| 60 | |
| 61 | |
| 62 | |
| 63 | |
| 64 | |
| 65 |