Skip to content

Update Rust implementations for Verlet algorithms #297

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Aug 30, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion contents/verlet_integration/code/julia/verlet.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
27 changes: 18 additions & 9 deletions contents/verlet_integration/code/rust/verlet.rs
Original file line number Diff line number Diff line change
@@ -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;

Expand All @@ -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;

Expand All @@ -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);
}
26 changes: 11 additions & 15 deletions contents/verlet_integration/verlet_integration.md
Original file line number Diff line number Diff line change
@@ -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:

Expand All @@ -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 %}
Expand Down Expand Up @@ -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)
Expand All @@ -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" %}
Expand Down Expand Up @@ -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 %}
Expand All @@ -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}
Expand All @@ -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}
Expand Down Expand Up @@ -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

Expand Down