From fa89a9b2369c00eb8920feec624e5586bde314e1 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 24 Jul 2018 12:23:14 -0400 Subject: [PATCH 1/4] Adjust Rust code for Verlet according to #257 --- .../verlet_integration/code/rust/verlet.rs | 27 ++++++++++++------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/contents/verlet_integration/code/rust/verlet.rs b/contents/verlet_integration/code/rust/verlet.rs index 4793afda7..e1518a223 100644 --- a/contents/verlet_integration/code/rust/verlet.rs +++ b/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); } From b988523212ba6a1b88c7da5954c6b95f49ce1ac3 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 24 Jul 2018 12:23:57 -0400 Subject: [PATCH 2/4] Operator whitespace in Julia Verlet implementation --- contents/verlet_integration/code/julia/verlet.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/contents/verlet_integration/code/julia/verlet.jl b/contents/verlet_integration/code/julia/verlet.jl index 0b0dd6a68..b55b85a9e 100644 --- a/contents/verlet_integration/code/julia/verlet.jl +++ b/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 From caf336c5851f3d5f19d2a4f732c4a8a5be0a212a Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 24 Jul 2018 12:30:06 -0400 Subject: [PATCH 3/4] Changed boundaries for Rust code in Verlet integration --- contents/verlet_integration/verlet_integration.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/contents/verlet_integration/verlet_integration.md b/contents/verlet_integration/verlet_integration.md index e1b28c7c9..613907087 100644 --- a/contents/verlet_integration/verlet_integration.md +++ b/contents/verlet_integration/verlet_integration.md @@ -101,7 +101,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) {% endmethod %} @@ -158,7 +158,7 @@ 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) {% endmethod %} From 73e4458569bb3b08e9a6b947b790f0cfe1056a57 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Tue, 24 Jul 2018 12:30:40 -0400 Subject: [PATCH 4/4] Minor grammar changes in Verlet integration text --- .../verlet_integration/verlet_integration.md | 22 ++++++++----------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/contents/verlet_integration/verlet_integration.md b/contents/verlet_integration/verlet_integration.md index 613907087..dfe5fbfd4 100644 --- a/contents/verlet_integration/verlet_integration.md +++ b/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 %} @@ -58,9 +58,7 @@ Unfortunately, this has not yet been implemented in LabVIEW, so here's Julia cod [import:1-13, lang:"rust"](code/rust/verlet.rs) {% 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) @@ -72,8 +70,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" %} @@ -109,7 +106,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} @@ -119,7 +116,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} @@ -161,8 +158,7 @@ Unfortunately, this has not yet been implemented in LabVIEW, so here's Julia cod [import:34-45, lang:"rust"](code/rust/verlet.rs) {% 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