Update Rust implementations for Verlet algorithms (algorithm-archivists#297)

Author: berquist

contents/verlet_integration/code/julia/verlet.jl

@@ -24,7 +24,7 @@ function stormer_verlet(pos::Float64, acc::Float64, dt::Float64)
         prev_pos = temp_pos
 
         # Because acceleration is constant, velocity is straightforward
-        vel += acc*dt
+        vel += acc * dt
     end
 
     return time, vel

contents/verlet_integration/code/rust/verlet.rs

@@ -1,4 +1,4 @@
-fn verlet(mut pos: f64, acc: f64, dt: f64) {
+fn verlet(mut pos: f64, acc: f64, dt: f64) -> f64 {
     let mut prev_pos = pos;
     let mut time = 0.0;
 
@@ -9,24 +9,29 @@ fn verlet(mut pos: f64, acc: f64, dt: f64) {
         prev_pos = temp_pos;
     }
 
-    println!("{}", time);
+    time
 }
 
-fn stormer_verlet(mut pos: f64, acc: f64, dt: f64) {
+fn stormer_verlet(mut pos: f64, acc: f64, dt: f64) -> (f64, f64) {
     let mut prev_pos = pos;
     let mut time = 0.0;
+    let mut vel = 0.0;
 
     while pos > 0.0 {
         time += dt;
         let temp_pos = pos;
         pos = pos * 2.0 - prev_pos + acc * dt * dt;
         prev_pos = temp_pos;
+
+        // Because acceleration is constant, velocity is
+        // straightforward
+        vel += acc * dt;
     }
 
-    println!("{}", time);
+    (time, vel)
 }
 
-fn velocity_verlet(mut pos: f64, acc: f64, dt: f64) {
+fn velocity_verlet(mut pos: f64, acc: f64, dt: f64) -> (f64, f64) {
     let mut time = 0.0;
     let mut vel = 0.0;
 
@@ -36,11 +41,15 @@ fn velocity_verlet(mut pos: f64, acc: f64, dt: f64) {
         vel += acc * dt;
     }
 
-    println!("{}", time);
+    (time, vel)
 }
 
 fn main() {
-    verlet(5.0, -10.0, 0.01);
-    stormer_verlet(5.0, -10.0, 0.01);
-    velocity_verlet(5.0, -10.0, 0.01);
+    let time_v = verlet(5.0, -10.0, 0.01);
+    let (time_sv, vel_sv) = stormer_verlet(5.0, -10.0, 0.01);
+    let (time_vv, vel_vv) = velocity_verlet(5.0, -10.0, 0.01);
+
+    println!("Time for original Verlet integration: {}", time_v);
+    println!("Time and velocity for Stormer Verlet integration: {}, {}", time_sv, vel_sv);
+    println!("Time and velocity for velocity Verlet integration: {}, {}", time_vv, vel_vv);
 }

contents/verlet_integration/verlet_integration.md

@@ -1,12 +1,12 @@
 # Verlet Integration
 
-Verlet Integration is essentially a solution to the kinematic equation for the motion of any object,
+Verlet integration is essentially a solution to the kinematic equation for the motion of any object,
 
 $$
 x = x_0 + v_0t + \frac{1}{2}at^2 + \frac{1}{6}bt^3 + \cdots
 $$
 
-Where $$x$$ is the position, $$v$$ is the velocity, $$a$$ is the acceleration, $$b$$ is the often forgotten jerk term, and $$t$$ is time. This equation is a central equation to almost every Newtonian physics solver and brings up a class of algorithms known as *force integrators*. One of the first force integrators to work with is *Verlet Integration*.
+where $$x$$ is the position, $$v$$ is the velocity, $$a$$ is the acceleration, $$b$$ is the often forgotten jerk term, and $$t$$ is time. This equation is a central equation to almost every Newtonian physics solver and brings up a class of algorithms known as *force integrators*. One of the first force integrators to work with is *Verlet Integration*.
 
 So, let's say we want to solve for the next timestep in $$x$$. To a close approximation (actually performing a Taylor Series Expansion about $$x(t\pm \Delta t)$$), that might look like this:
 
@@ -20,13 +20,13 @@ $$
 x(t-\Delta t) = x(t) - v(t)\Delta t + \frac{1}{2}a(t)\Delta t^2 - \frac{1}{6}b(t) \Delta t^3 + \mathcal{O}(\Delta t^4)
 $$
 
-Now, we have two equations to solve for two different timesteps in x, one of which we already have. If we add the two equations together and solve for $$x(t+\Delta t)$$, we find
+Now, we have two equations to solve for two different timesteps in x, one of which we already have. If we add the two equations together and solve for $$x(t+\Delta t)$$, we find that
 
 $$
 x(t+ \Delta t) = 2x(t) - x(t-\Delta t) + a(t)\Delta t^2 + \mathcal{O}(\Delta t^4)
 $$
 
-So, this means, we can find our next $$x$$ simply by knowing our current $$x$$, the $$x$$ before that, and the acceleration! No velocity necessary! In addition, this drops the error to $$\mathcal{O}(\Delta t^4)$$, which is great!
+So, this means we can find our next $$x$$ simply by knowing our current $$x$$, the $$x$$ before that, and the acceleration! No velocity necessary! In addition, this drops the error to $$\mathcal{O}(\Delta t^4)$$, which is great!
 Here is what it looks like in code:
 
 {% method %}
@@ -60,9 +60,7 @@ Unfortunately, this has not yet been implemented in LabVIEW, so here's Julia cod
 [import:1-15, lang:"swift"](code/swift/verlet.swift)
 {% endmethod %}
 
-
-
-Now, obviously this poses a problem, what if we want to calculate a term that requires velocity, like the kinetic energy, $$\frac{1}{2}mv^2$$? In this case, we certainly cannot get rid of the velocity! Well, we can find the velocity to $$\mathcal{O}(\Delta t^2)$$ accuracy by using the Stormer-Verlet method, which is the same as before, but we calculate velocity like so
+Now, obviously this poses a problem; what if we want to calculate a term that requires velocity, like the kinetic energy, $$\frac{1}{2}mv^2$$? In this case, we certainly cannot get rid of the velocity! Well, we can find the velocity to $$\mathcal{O}(\Delta t^2)$$ accuracy by using the Stormer-Verlet method, which is the same as before, but we calculate velocity like so
 
 $$
 v(t) = \frac{x(t+\Delta t) - x(t-\Delta t)}{2\Delta t} + \mathcal{O}(\Delta t^2)
@@ -74,8 +72,7 @@ $$
 v(t+\Delta t) = \frac{x(t+\Delta t) - x(t)}{\Delta t} + \mathcal{O}(\Delta t)
 $$
 
-However, the error for this is $$\mathcal{O}(\Delta t)$$, which is quite poor, but get's the job done in a pinch.
-Here's what it looks like in code:
+However, the error for this is $$\mathcal{O}(\Delta t)$$, which is quite poor, but it gets the job done in a pinch.  Here's what it looks like in code:
 
 {% method %}
 {% sample lang="jl" %}
@@ -103,7 +100,7 @@ Unfortunately, this has not yet been implemented in LabVIEW, so here's Julia cod
 {% sample lang="javascript" %}
 [import:18-35, lang:"javascript"](code/javascript/verlet.js)
 {% sample lang="rs" %}
-[import:15-27, lang:"rust"](code/rust/verlet.rs)
+[import:15-32, lang:"rust"](code/rust/verlet.rs)
 {% sample lang="swift" %}
 [import:17-34, lang:"swift"](code/swift/verlet.swift)
 {% endmethod %}
@@ -113,7 +110,7 @@ Now, let's say we actually need the velocity to calculate out next timestep. Wel
 
 # Velocity Verlet
 
-In some ways, this algorithm is even simpler than above. We can calculate everything like so
+In some ways, this algorithm is even simpler than above. We can calculate everything like
 
 $$
 \begin{align}
@@ -123,7 +120,7 @@ v(t+\Delta t) &= v(t) + \frac{1}{2}(a(t) + a(t+\Delta t))\Delta t
 \end{align}
 $$
 
-Which is literally the kinematic equation above, solving for $$x$$, $$v$$, and $$a$$ every timestep. You can also split up the equations like so
+which is literally the kinematic equation above, solving for $$x$$, $$v$$, and $$a$$ every timestep. You can also split up the equations like so
 
 $$
 \begin{align}
@@ -162,13 +159,12 @@ Unfortunately, this has not yet been implemented in LabVIEW, so here's Julia cod
 {% sample lang="javascript" %}
 [import:37-50, lang:"javascript"](code/javascript/verlet.js)
 {% sample lang="rs" %}
-[import:29-40, lang:"rust"](code/rust/verlet.rs)
+[import:34-45, lang:"rust"](code/rust/verlet.rs)
 {% sample lang="swift" %}
 [import:36-49, lang:"swift"](code/swift/verlet.swift)
 {% endmethod %}
 
-
-Even though this method is more used than the simple Verlet method mentioned above, it unfortunately has an error term of $$\mathcal{O} \Delta t^2$$, which is two orders of magnitude worse. That said, if you want to have a simulation with many objects that depend on one another --- like a gravity simulation --- the Velocity Verlet algorithm is a handy choice; however, you may have to play further tricks to allow everything to scale appropriately. These types of simulations are sometimes called *n-body* simulations and one such trick is the Barnes-Hut algorithm, which cuts the complexity of n-body simulations from $$\sim \mathcal{O}(n^2)$$ to $$\sim \mathcal{O}(n\log(n))$$
+Even though this method is more widely used than the simple Verlet method mentioned above, it unforunately has an error term of $$\mathcal{O}(\Delta t^2)$$, which is two orders of magnitude worse. That said, if you want to have a simulaton with many objects that depend on one another --- like a gravity simulation --- the Velocity Verlet algorithm is a handy choice; however, you may have to play further tricks to allow everything to scale appropriately. These types of simulatons are sometimes called *n-body* simulations and one such trick is the Barnes-Hut algorithm, which cuts the complexity of n-body simulations from $$\sim \mathcal{O}(n^2)$$ to $$\sim \mathcal{O}(n\log(n))$$.
 
 ## Example Code