Verlet Integration

Verlet integration is essentially a solution to the kinematic equation for the motion of any object,

where is the position, is the velocity, is the acceleration, is the often forgotten jerk term, and 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 . To a close approximation (actually performing a Taylor Series Expansion about ), that might look like this:

This means that if we need to find the next , we need the current , , , etc. However, because few people calculate the jerk term, our error is typically . That said, we can calculate with less knowledge and higher accuracy if we play a trick! Let's say we want to calculate of the previous timestep. Again, to a close approximation, that might look like this:

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 , we find that

So, this means we can find our next simply by knowing our current , the before that, and the acceleration! No velocity necessary! In addition, this drops the error to , which is great! Here is what it looks like in code:

function verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0

while (pos > 0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos
end

return time
end
double verlet(double pos, double acc, double dt) {

double prev_pos = pos;
double time = 0;

while (pos > 0) {
time += dt;
double next_pos = pos * 2 - prev_pos + acc * dt * dt;
prev_pos = pos;
pos = next_pos;
}

return time;
}
void verlet(double *time, double pos, double acc, double dt) {
double prev_pos, temp_pos;
prev_pos = pos;
*time = 0.0;

while (pos > 0) {
*time += dt;
temp_pos = pos;
pos = pos * 2 - prev_pos + acc * dt * dt;
prev_pos = temp_pos;
}
}
static double verlet(double pos, double acc, double dt) {

// Note that we are using a temp variable for the previous position
double prev_pos, temp_pos, time;
prev_pos = pos;
time = 0;

while (pos > 0) {
time += dt;
temp_pos = pos;
pos = pos*2 - prev_pos + acc * dt * dt;
prev_pos = temp_pos;
}

return time;
}
def verlet(pos, acc, dt):
prev_pos = pos
time = 0

while pos > 0:
time += dt
next_pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos, pos = pos, next_pos

return time
type Method = Model -> Time -> Particle -> Particle -> Particle

verlet :: Method
verlet acc dt (xOld, _, _, _) (x, _, a, t) = (x', v', a', t + dt)
where
x' = 2 * x - xOld + a * dt ^ 2
v' = 0
a' = acc (x', v', a, t + dt)

Unfortunately, this has not yet been implemented in scratch, so here's Julia code:

function verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0

while (pos > 0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos
end

return time
end

Unfortunately, this has not yet been implemented in matlab, so here's Julia code:

function verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0

while (pos > 0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos
end

return time
end

Unfortunately, this has not yet been implemented in LabVIEW, so here's Julia code:

function verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0

while (pos > 0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos
end

return time
end
function verlet(pos, acc, dt) {
let prevPos = pos;
let time = 0;
let tempPos;

while (pos > 0) {
time += dt;
tempPos = pos;
pos = pos * 2 - prevPos + acc * dt * dt;
prevPos = tempPos;
}

return time;
}
fn verlet(mut pos: f64, acc: f64, dt: f64) -> f64 {
let mut prev_pos = pos;
let mut time = 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;
}

time
}
func verlet(pos: Double, acc: Double, dt: Double) -> Double {
var pos = pos
var temp_pos, time: Double
var prev_pos = pos
time = 0.0

while (pos > 0) {
time += dt
temp_pos = pos
pos = pos*2 - prev_pos + acc * dt * dt
prev_pos = temp_pos
}

return time
}
SUBROUTINE verlet(pos, acc, dt, time)
IMPLICIT NONE
REAL(8), INTENT(INOUT) :: pos, acc, dt, time
REAL(8)                :: prev_pos, next_pos

prev_pos = pos
time     = 0d0

DO
IF (pos > 0d0) THEN
time     = time + dt
next_pos = pos * 2d0 - prev_pos + acc * dt ** 2
prev_pos = pos
pos      = next_pos
ELSE
EXIT
END IF
END DO
END SUBROUTINE verlet
def verlet(pos, acc, dt)

prev_pos = pos
time = 0
while pos > 0 do
time += dt
temp_pos = pos
pos = pos*2 - prev_pos + acc * dt * dt
prev_pos = temp_pos
end

return time

end
func verlet(pos, acc, dt float64) (time float64) {
prevPos := pos
time = 0

for pos > 0 {
time += dt
nextPos := pos*2 - prevPos + acc*dt*dt
prevPos, pos = pos, nextPos
}

return
}
# xmm0 - pos
# xmm1 - acc
# xmm2 - dt
# RET xmm0 - time
verlet:
pxor   xmm7, xmm7                  # Holds 0 for comparisons
pxor   xmm3, xmm3                  # Holds time value
comisd xmm0, xmm7                  # Check if pos is greater then 0.0
jbe    verlet_return
movsd  xmm6, xmm1                  # xmm6 = acc * dt * dt
mulsd  xmm6, xmm2
mulsd  xmm6, xmm2
movsd  xmm5, xmm0                  # Holds previous position
verlet_loop:
addsd  xmm3, xmm2                  # Adding dt to time
movsd  xmm4, xmm0                  # Hold old value of posistion
addsd  xmm0, xmm0                  # Calculating new position
subsd  xmm0, xmm5
movsd  xmm5, xmm4
comisd xmm0, xmm7                  # Check if position is greater then 0.0
ja     verlet_loop
verlet_return:
movsd  xmm0, xmm3                  # Saving time value
ret
fun verlet(_pos: Double, acc: Double, dt: Double): Double {
var pos = _pos  // Since function parameter are val and can't be modified
var prevPos = pos
var time = 0.0

while (pos > 0) {
time += dt
val nextPos = pos * 2 - prevPos + acc * dt * dt
prevPos = pos
pos = nextPos
}
return time
}
proc verlet(pos_in, acc, dt: float): float =
var
pos: float = pos_in
prevPos: float = pos
time: float = 0.0
tempPos: float

while pos > 0.0:
time += dt
tempPos = pos
pos = pos * 2 - prevPos + acc * dt * dt
prevPos = tempPos

return time
(defun verlet (pos acc dt)
"Integrates Newton's equation for motion while pos > 0 using Verlet integration."
(loop
with prev-pos = pos
for time = 0 then (incf time dt)
while (> pos 0)
;; The starting speed is assumed to be zero.
do (psetf
pos (+ (* pos 2) (- prev-pos) (* acc dt dt))
prev-pos pos)
finally (return time)))

Now, obviously this poses a problem; what if we want to calculate a term that requires velocity, like the kinetic energy, ? In this case, we certainly cannot get rid of the velocity! Well, we can find the velocity to accuracy by using the Stormer-Verlet method, which is the same as before, but we calculate velocity like so

Note that the 2 in the denominator appears because we are going over 2 timesteps. It's essentially solving . In addition, we can calculate the velocity of the next timestep like so

However, the error for this is , which is quite poor, but it gets the job done in a pinch. Here's what it looks like in code:

function stormer_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos

# Because acceleration is constant, velocity is straightforward
vel += acc * dt
end

return time, vel
end
timestep stormer_verlet(double pos, double acc, double dt) {

double prev_pos = pos;
double time = 0;
double vel = 0;
while (pos > 0) {
time += dt;
double next_pos = pos * 2 - prev_pos + acc * dt * dt;
prev_pos = pos;
pos = next_pos;

// The acceleration is constant, so the velocity is
// straightforward
vel += acc * dt;
}

return timestep { time, vel };
}
void stormer_verlet(double *time, double *vel,
double pos, double acc, double dt) {
double prev_pos, temp_pos;
prev_pos = pos;
*vel = 0.0;
*time = 0.0;

while (pos > 0) {
*time += dt;
temp_pos = pos;
pos = pos * 2 - prev_pos + acc * dt * dt;
prev_pos = temp_pos;

*vel += acc * dt;
}
}
static VerletValues stormer_verlet(double pos, double acc, double dt) {

// Note that we are using a temp variable for the previous position
double prev_pos, temp_pos, time, vel;
prev_pos = pos;
vel = 0;
time = 0;
while (pos > 0) {
time += dt;
temp_pos = pos;
pos = pos*2 - prev_pos + acc * dt * dt;
prev_pos = temp_pos;

// The acceleration is constant, so the velocity is straightforward
vel += acc*dt;
}

return new VerletValues(time, vel);
}
def stormer_verlet(pos, acc, dt):
prev_pos = pos
time = 0
vel = 0

while pos > 0:
time += dt
next_pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos, pos = pos, next_pos
vel += acc * dt

return time, vel
stormerVerlet :: Method
stormerVerlet acc dt (xOld, _, _, _) (x, _, a, t) = (x', v', a', t + dt)
where
x' = 2 * x - xOld + a * dt ^ 2
v' = (x' - x) / dt
a' = acc (x', v', a, t + dt)

Unfortunately, this has not yet been implemented in scratch, so here's Julia code:

function stormer_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos

# Because acceleration is constant, velocity is straightforward
vel += acc * dt
end

return time, vel
end

Unfortunately, this has not yet been implemented in matlab, so here's Julia code:

function stormer_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos

# Because acceleration is constant, velocity is straightforward
vel += acc * dt
end

return time, vel
end

Unfortunately, this has not yet been implemented in LabVIEW, so here's Julia code:

function stormer_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos

# Because acceleration is constant, velocity is straightforward
vel += acc * dt
end

return time, vel
end
function stormerVerlet(pos, acc, dt) {
let prevPos = pos;
let time = 0;
let vel = 0;
let tempPos;

while (pos > 0) {
time += dt;
tempPos = pos;
pos = pos * 2 - prevPos + acc * dt * dt;
prevPos = tempPos;

vel += acc * dt;
}

return { time, vel };
}
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;
}

(time, vel)
}
func stormerVerlet(pos: Double, acc: Double, dt: Double) -> (time: Double, vel: Double) {
var pos = pos
var temp_pos, time, vel: Double
var prev_pos = pos
vel = 0
time = 0

while (pos > 0) {
time += dt
temp_pos = pos
pos = pos*2 - prev_pos + acc * dt * dt
prev_pos = temp_pos

vel += acc*dt
}

return (time:time, vel:vel)
}
SUBROUTINE stormer_verlet(pos, acc, dt, time, vel)
IMPLICIT NONE
REAL(8), INTENT(INOUT) :: pos, acc, dt, time, vel
REAL(8)                :: prev_pos, next_pos

prev_pos = pos
time     = 0d0
vel      = 0d0

DO
IF (pos > 0d0) THEN
time     = time + dt
next_pos = pos * 2 - prev_pos + acc * dt ** 2
prev_pos = pos
pos      = next_pos
vel      = vel + acc * dt
ELSE
EXIT
END IF
END DO
END SUBROUTINE stormer_verlet
def stormer_verlet(pos, acc, dt)

prev_pos = pos
vel = 0
time = 0
while pos > 0 do
time += dt
temp_pos = pos
pos = pos*2 - prev_pos + acc * dt * dt
prev_pos = temp_pos

vel += acc*dt
end

return time

end
func stormerVerlet(pos, acc, dt float64) (time, vel float64) {
prevPos := pos
time, vel = 0, 0

for pos > 0 {
time += dt
vel += acc * dt
nextPos := pos*2 - prevPos + acc*dt*dt
prevPos, pos = pos, nextPos
}

return
}
# xmm0 - pos
# xmm1 - acc
# xmm2 - dt
# RET xmm0 - time
# RET xmm1 - velocity
stormer_verlet:
pxor   xmm7, xmm7                  # Holds 0 for comparisons
pxor   xmm3, xmm3                  # Holds time value
comisd xmm0, xmm7                  # Check if pos is greater then 0.0
jbe    stormer_verlet_return
movsd  xmm6, xmm1                  # xmm6 = acc * dt * dt
mulsd  xmm6, xmm2
mulsd  xmm6, xmm2
movsd  xmm5, xmm0                  # Holds previous position
stormer_verlet_loop:
addsd  xmm3, xmm2                  # Adding dt to time
movsd  xmm4, xmm0                  # Hold old value of posistion
addsd  xmm0, xmm0                  # Calculating new position
subsd  xmm0, xmm5
movsd  xmm5, xmm4
comisd xmm0, xmm7                  # Check if position is greater then 0.0
ja     stormer_verlet_loop
stormer_verlet_return:
movsd  xmm0, xmm3                  # Saving time and velocity
mulsd  xmm3, xmm1
movsd  xmm1, xmm3
ret
fun stormerVerlet(_pos: Double, acc: Double, dt: Double): VerletValues {
var pos = _pos
var prevPos = pos
var time = 0.0
var vel = 0.0
while (pos > 0) {
time += dt
val nextPos = pos * 2 - prevPos + acc * dt * dt
prevPos = pos
pos = nextPos
vel += acc * dt
}
return VerletValues(time, vel)
}
proc stormerVerlet(pos_in, acc, dt: float): (float, float) =
var
pos: float = pos_in
prevPos: float = pos
time: float = 0.0
vel: float = 0.0
tempPos: float

while pos > 0.0:
time += dt
tempPos = pos
pos = pos * 2 - prevPos + acc * dt * dt
prevPos = tempPos

vel += acc * dt

return (time, vel)
(defun stormer-verlet (pos acc dt)
"Integrates Newton's equation for motion while pos > 0 using the Stormer-Verlet method."
(loop
with prev-pos = pos
for time = 0 then (incf time dt)
for vel = 0 then (incf vel (* acc dt))
while (> pos 0)
;; Variables are changed simultaneously by 'psetf', so there's no need for a temporary variable.
do (psetf
pos (+ (* pos 2) (- prev-pos) (* acc dt dt))
prev-pos pos)
finally (return (list time vel))))

Now, let's say we actually need the velocity to calculate out next timestep. Well, in this case, we simply cannot use the above approximation and instead need to use the Velocity Verlet algorithm.

Velocity Verlet

In some ways, this algorithm is even simpler than above. We can calculate everything like

which is literally the kinematic equation above, solving for , , and every timestep. You can also split up the equations like so

Here is the velocity Verlet method in code:

function velocity_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
end

return time, vel
end
timestep velocity_verlet(double pos, double acc, double dt) {

double time = 0;
double vel = 0;
while (pos > 0) {
time += dt;
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
}

return timestep { time, vel };
}
void velocity_verlet(double *time, double *vel,
double pos, double acc, double dt) {
*vel = 0.0;
*time = 0.0;

while (pos > 0) {
*time += dt;
pos += (*vel) * dt + 0.5 * acc * dt * dt;
*vel += acc * dt;
}
}
static VerletValues velocity_verlet(double pos, double acc, double dt) {

// Note that we are using a temp variable for the previous position
double time, vel;
vel = 0;
time = 0;
while (pos > 0) {
time += dt;
pos += vel*dt + 0.5*acc * dt * dt;
vel += acc*dt;
}
return new VerletValues(time, vel);
}
def velocity_verlet(pos, acc, dt):
time = 0
vel = 0

while pos > 0:
time += dt
pos += vel * dt + 0.5 * acc * dt * dt
vel += acc * dt

return time, vel
velocityVerlet :: Method
velocityVerlet acc dt (xOld, _, aOld, _) (x, v, a, t) = (x', v', a', t + dt)
where
x' = 2 * x - xOld + a * dt ^ 2
v' = v + 0.5 * (aOld + a) * dt
a' = acc (x', v', a, t + dt)

Unfortunately, this has not yet been implemented in scratch, so here's Julia code:

function velocity_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
end

return time, vel
end

Unfortunately, this has not yet been implemented in matlab, so here's Julia code:

function velocity_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
end

return time, vel
end

Unfortunately, this has not yet been implemented in LabVIEW, so here's Julia code:

function velocity_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
end

return time, vel
end
function velocityVerlet(pos, acc, dt) {
let time = 0;
let vel = 0;

while (pos > 0) {
time += dt;
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
}

return { time, vel };
}
fn velocity_verlet(mut pos: f64, acc: f64, dt: f64) -> (f64, f64) {
let mut time = 0.0;
let mut vel = 0.0;

while pos > 0.0 {
time += dt;
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
}

(time, vel)
}
func velocityVerlet(pos: Double, acc: Double, dt: Double) -> (time: Double, vel: Double) {
var pos = pos
var time, vel : Double
vel = 0
time = 0

while (pos > 0) {
time += dt
pos += vel*dt + 0.5*acc * dt * dt
vel += acc*dt
}

return (time:time, vel:vel)
}
SUBROUTINE velocity_verlet(pos, acc, dt, time, vel)
IMPLICIT NONE
REAL(8), INTENT(INOUT) :: pos, acc, dt, time, vel

time     = 0d0
vel      = 0d0

DO
IF (pos > 0d0) THEN
time = time + dt
pos  = pos + vel * dt + 0.5d0 * acc * dt ** 2
vel  = vel + acc * dt
ELSE
EXIT
END IF
END DO
END SUBROUTINE velocity_verlet
def velocity_verlet(pos, acc, dt)

vel = 0
time = 0
while pos > 0 do
time += dt
pos += vel*dt + 0.5*acc * dt * dt
vel += acc*dt
end

return time

end
func velocityVerlet(pos, acc, dt float64) (time, vel float64) {
time, vel = 0, 0

for pos > 0 {
time += dt
pos += vel*dt + .5*acc*dt*dt
vel += acc * dt
}

return
}
# xmm0 - pos
# xmm1 - acc
# xmm2 - dt
# RET xmm0 - time
# RET xmm1 - velocity
velocity_verlet:
pxor   xmm7, xmm7                  # Holds 0 for comparisons
pxor   xmm3, xmm3                  # Holds the velocity value
pxor   xmm4, xmm4                  # Holds the time value
comisd xmm0, xmm7                  # Check if pos is greater then 0.0
jbe    velocity_verlet_return
movsd  xmm5, half                  # xmm5 = 0.5 * dt * dt * acc
mulsd  xmm5, xmm2
mulsd  xmm5, xmm2
mulsd  xmm5, xmm1
velocity_verlet_loop:
movsd  xmm6, xmm3                  # Move velocity into register
mulsd  xmm6, xmm2                  # Calculate new position
addsd  xmm4, xmm2                  # Incrementing time
movsd  xmm3, xmm4                  # Updating velocity
mulsd  xmm3, xmm1
comisd xmm0, xmm7
ja     velocity_verlet_loop
velocity_verlet_return:
movsd  xmm0, xmm4                  # Saving time and velocity
movsd  xmm1, xmm3
ret
fun velocityVerlet(_pos: Double, acc: Double, dt: Double): VerletValues {
var pos = _pos
var time = 0.0
var vel = 0.0
while (pos > 0) {
time += dt
pos += vel * dt + 0.5 * acc * dt * dt
vel += acc * dt
}
return VerletValues(time, vel)
}
proc velocityVerlet(pos_in, acc, dt: float): (float, float) =
var
pos: float = pos_in
time: float = 0.0
vel: float = 0.0

while pos > 0.0:
time += dt
pos += vel * dt + 0.5 * acc * dt * dt
vel += acc * dt

return (time, vel)
(defun velocity-verlet (pos acc dt)
"Integrates Newton's equation for motion while pos > 0 using the velocity in calculations."
(loop
for time = 0 then (incf time dt)
for vel = 0 then (incf vel (* acc dt))
for p = pos then (incf p (+ (* vel dt) (* 0.5 acc dt dt)))
while (> p 0)
finally (return (list time vel))))

Even though this method is more widely used than the simple Verlet method mentioned above, it unforunately has an error term of , 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 to .

Video Explanation

Here is a video describing Verlet integration:

Example Code

Both of these methods work simply by iterating timestep-by-timestep and can be written straightforwardly in any language. For reference, here are snippets of code that use both the classic and velocity Verlet methods to find the time it takes for a ball to hit the ground after being dropped from a given height.

function verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0

while (pos > 0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos
end

return time
end

function stormer_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
temp_pos = pos
pos = pos * 2 - prev_pos + acc * dt * dt
prev_pos = temp_pos

# Because acceleration is constant, velocity is straightforward
vel += acc * dt
end

return time, vel
end

function velocity_verlet(pos::Float64, acc::Float64, dt::Float64)
prev_pos = pos
time = 0.0
vel = 0.0

while (pos > 0.0)
time += dt
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
end

return time, vel
end

function main()
time = verlet(5.0, -10.0, 0.01);
println("Time for Verlet integration is: $(time)\n") time, vel = stormer_verlet(5.0, -10.0, 0.01); println("Time for Stormer Verlet integration is:$(time)")
println("Velocity for Stormer Verlet integration is: $(vel)\n") time, vel = velocity_verlet(5.0, -10.0, 0.01); println("Time for velocity Verlet integration is:$(time)")
println("Velocity for velocity Verlet integration is: $(vel)\n") end main() #include <iomanip> #include <iostream> struct timestep { double time; double vel; }; double verlet(double pos, double acc, double dt) { double prev_pos = pos; double time = 0; while (pos > 0) { time += dt; double next_pos = pos * 2 - prev_pos + acc * dt * dt; prev_pos = pos; pos = next_pos; } return time; } timestep stormer_verlet(double pos, double acc, double dt) { double prev_pos = pos; double time = 0; double vel = 0; while (pos > 0) { time += dt; double next_pos = pos * 2 - prev_pos + acc * dt * dt; prev_pos = pos; pos = next_pos; // The acceleration is constant, so the velocity is // straightforward vel += acc * dt; } return timestep { time, vel }; } timestep velocity_verlet(double pos, double acc, double dt) { double time = 0; double vel = 0; while (pos > 0) { time += dt; pos += vel * dt + 0.5 * acc * dt * dt; vel += acc * dt; } return timestep { time, vel }; } int main() { std::cout << std::fixed << std::setprecision(8); // Note that depending on the simulation, you might want to have the // Verlet loop outside. // For example, if your acceleration chages as a function of time, // you might need to also change the acceleration to be read into // each of these functions. double time = verlet(5.0, -10, 0.01); std::cout << "Time for Verlet integration is: " \ << time << std::endl; timestep timestep_sv = stormer_verlet(5.0, -10, 0.01); std::cout << "Time for Stormer Verlet integration is: " \ << timestep_sv.time << std::endl; std::cout << "Velocity for Stormer Verlet integration is: " \ << timestep_sv.vel << std::endl; timestep timestep_vv = velocity_verlet(5.0, -10, 0.01); std::cout << "Time for velocity Verlet integration is: " \ << timestep_vv.time << std::endl; std::cout << "Velocity for velocity Verlet integration is: " \ << timestep_vv.vel << std::endl; return 0; } #include <stdio.h> void verlet(double *time, double pos, double acc, double dt) { double prev_pos, temp_pos; prev_pos = pos; *time = 0.0; while (pos > 0) { *time += dt; temp_pos = pos; pos = pos * 2 - prev_pos + acc * dt * dt; prev_pos = temp_pos; } } void stormer_verlet(double *time, double *vel, double pos, double acc, double dt) { double prev_pos, temp_pos; prev_pos = pos; *vel = 0.0; *time = 0.0; while (pos > 0) { *time += dt; temp_pos = pos; pos = pos * 2 - prev_pos + acc * dt * dt; prev_pos = temp_pos; *vel += acc * dt; } } void velocity_verlet(double *time, double *vel, double pos, double acc, double dt) { *vel = 0.0; *time = 0.0; while (pos > 0) { *time += dt; pos += (*vel) * dt + 0.5 * acc * dt * dt; *vel += acc * dt; } } int main() { double time, vel; verlet(&time, 5.0, -10, 0.01); printf("Time for Verlet integration is: %lf\n", time); stormer_verlet(&time, &vel, 5.0, -10, 0.01); printf("Time and velocity for Stormer Verlet integration is: %lf, %lf\n", time, vel); velocity_verlet(&time, &vel, 5.0, -10, 0.01); printf("Time and velocity for velocity Verlet integration is: %lf, %lf\n", time, vel); return 0; } public class VerletValues { public double time; public double vel; public VerletValues(double time, double vel) { this.time = time; this.vel = vel; } } public class Verlet { static double verlet(double pos, double acc, double dt) { // Note that we are using a temp variable for the previous position double prev_pos, temp_pos, time; prev_pos = pos; time = 0; while (pos > 0) { time += dt; temp_pos = pos; pos = pos*2 - prev_pos + acc * dt * dt; prev_pos = temp_pos; } return time; } static VerletValues stormer_verlet(double pos, double acc, double dt) { // Note that we are using a temp variable for the previous position double prev_pos, temp_pos, time, vel; prev_pos = pos; vel = 0; time = 0; while (pos > 0) { time += dt; temp_pos = pos; pos = pos*2 - prev_pos + acc * dt * dt; prev_pos = temp_pos; // The acceleration is constant, so the velocity is straightforward vel += acc*dt; } return new VerletValues(time, vel); } static VerletValues velocity_verlet(double pos, double acc, double dt) { // Note that we are using a temp variable for the previous position double time, vel; vel = 0; time = 0; while (pos > 0) { time += dt; pos += vel*dt + 0.5*acc * dt * dt; vel += acc*dt; } return new VerletValues(time, vel); } public static void main(String[] args) { double verletTime = verlet(5.0, -10, 0.01); System.out.println("Time for Verlet integration is: " + verletTime); VerletValues stormerVerlet = stormer_verlet(5.0, -10, 0.01); System.out.println("Time for Stormer Verlet integration is: " + stormerVerlet.time); System.out.println("Velocity for Stormer Verlet integration is: " + stormerVerlet.vel); VerletValues velocityVerlet = velocity_verlet(5.0, -10, 0.01); System.out.println("Time for velocity Verlet integration is: " + velocityVerlet.time); System.out.println("Velocity for velocity Verlet integration is: " + velocityVerlet.vel); } } def verlet(pos, acc, dt): prev_pos = pos time = 0 while pos > 0: time += dt next_pos = pos * 2 - prev_pos + acc * dt * dt prev_pos, pos = pos, next_pos return time def stormer_verlet(pos, acc, dt): prev_pos = pos time = 0 vel = 0 while pos > 0: time += dt next_pos = pos * 2 - prev_pos + acc * dt * dt prev_pos, pos = pos, next_pos vel += acc * dt return time, vel def velocity_verlet(pos, acc, dt): time = 0 vel = 0 while pos > 0: time += dt pos += vel * dt + 0.5 * acc * dt * dt vel += acc * dt return time, vel def main(): time = verlet(5, -10, 0.01) print("Verlet") print("Time: {:.10f}".format(time)) print() time, vel = stormer_verlet(5, -10, 0.01) print("Stormer-Verlet") print("Time: {:.10f}".format(time)) print("Velocity: {:.10f}".format(vel)) print() time, vel = velocity_verlet(5, -10, 0.01) print("Velocity Verlet") print("Time: {:.10f}".format(time)) print("Velocity: {:.10f}".format(vel)) print() if __name__ == '__main__': main() -- submitted by Jie type Time = Double type Position = Double type Speed = Double type Acceleration = Double type Particle = (Position, Speed, Acceleration, Time) type Model = Particle -> Acceleration type Method = Model -> Time -> Particle -> Particle -> Particle verlet :: Method verlet acc dt (xOld, _, _, _) (x, _, a, t) = (x', v', a', t + dt) where x' = 2 * x - xOld + a * dt ^ 2 v' = 0 a' = acc (x', v', a, t + dt) stormerVerlet :: Method stormerVerlet acc dt (xOld, _, _, _) (x, _, a, t) = (x', v', a', t + dt) where x' = 2 * x - xOld + a * dt ^ 2 v' = (x' - x) / dt a' = acc (x', v', a, t + dt) velocityVerlet :: Method velocityVerlet acc dt (xOld, _, aOld, _) (x, v, a, t) = (x', v', a', t + dt) where x' = 2 * x - xOld + a * dt ^ 2 v' = v + 0.5 * (aOld + a) * dt a' = acc (x', v', a, t + dt) trajectory :: Method -> Model -> Time -> Particle -> [Particle] trajectory method acc dt p0@(x, v, a, t0) = traj where traj = p0 : p1 : zipWith (method acc dt) traj (tail traj) p1 = (x', v', acc (x', v', a, t0 + dt), t0 + dt) x' = x + v * dt + 0.5 * a * dt ^ 2 v' = v + a * dt main :: IO () main = do let p0 = (5, 0, -10, 0) dt = 0.001 freefall _ = -10 aboveGround (x, _, _, _) = x > 0 integrate m = last$ takeWhile aboveGround $trajectory m freefall dt p0 print$ integrate verlet
print $integrate stormerVerlet print$ integrate velocityVerlet

Submitted by Jie % Submitted by P. Mekhail
% Parameters to change
n = 400;        % Number of steps
x0 = 5;         % Ball starting height(in metres)
v0 = 0;         % Ball starting velocity (+ive is up)
dt = 0.01;      % Time step (in seconds)
eff = 0.4;      % Ball efficency when bouncing
A = @(x) -10;   % Acceleration as a function of position
bounce = 1;     % Do you want the ball to bounce?

% Making position and time vectors
x = zeros(n,1);
t = 0:dt:n*dt-dt;

% Setting the initial conditions
x(1) = x0;
x(2) = x0 + v0*dt + 0.5*A(x0)*dt^2;

% Runnin Verlet Integration
for i = 2:n-1
xnew = 2*x(i)-x(i-1)+A(x(i))*dt^2;
if bounce
if xnew > 0
% If you haven't hit the ground keep going
x(i+1) = xnew;
else
% If you have calculated velocity and invert its sign
v = sqrt(eff)*(xnew-x(i-1))/(2*dt);
x(i+1) = x(i) - v*dt + 0.5*A(x(i))*dt^2;
end
else
x(i+1) = xnew;
end
end
plot(t,x)
title('Ball''s Trajectory')
xlabel('Time (s)'); ylabel('Height (m)');

Submitted by P. Mekhail function verlet(pos, acc, dt) {
let prevPos = pos;
let time = 0;
let tempPos;

while (pos > 0) {
time += dt;
tempPos = pos;
pos = pos * 2 - prevPos + acc * dt * dt;
prevPos = tempPos;
}

return time;
}

function stormerVerlet(pos, acc, dt) {
let prevPos = pos;
let time = 0;
let vel = 0;
let tempPos;

while (pos > 0) {
time += dt;
tempPos = pos;
pos = pos * 2 - prevPos + acc * dt * dt;
prevPos = tempPos;

vel += acc * dt;
}

return { time, vel };
}

function velocityVerlet(pos, acc, dt) {
let time = 0;
let vel = 0;

while (pos > 0) {
time += dt;
pos += vel * dt + 0.5 * acc * dt * dt;
vel += acc * dt;
}

return { time, vel };
}

const time = verlet(5, -10, 0.01);
console.log(Time for Verlet integration is: ${time}\n); const stormer = stormerVerlet(5, -10, 0.01); console.log(Time for Stormer Verlet integration is:${stormer.time});
console.log(Velocity for Stormer Verlet integration is: ${stormer.vel}\n); const velocity = velocityVerlet(5, -10, 0.01); console.log(Time for Velocity Verlet integration is:${velocity.time});
console.log(Velocity for Velocity Verlet integration is: ${velocity.vel}\n); fn verlet(mut pos: f64, acc: f64, dt: f64) -> f64 { let mut prev_pos = pos; let mut time = 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; } time } 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; } (time, vel) } fn velocity_verlet(mut pos: f64, acc: f64, dt: f64) -> (f64, f64) { let mut time = 0.0; let mut vel = 0.0; while pos > 0.0 { time += dt; pos += vel * dt + 0.5 * acc * dt * dt; vel += acc * dt; } (time, vel) } fn main() { 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 ); } func verlet(pos: Double, acc: Double, dt: Double) -> Double { var pos = pos var temp_pos, time: Double var prev_pos = pos time = 0.0 while (pos > 0) { time += dt temp_pos = pos pos = pos*2 - prev_pos + acc * dt * dt prev_pos = temp_pos } return time } func stormerVerlet(pos: Double, acc: Double, dt: Double) -> (time: Double, vel: Double) { var pos = pos var temp_pos, time, vel: Double var prev_pos = pos vel = 0 time = 0 while (pos > 0) { time += dt temp_pos = pos pos = pos*2 - prev_pos + acc * dt * dt prev_pos = temp_pos vel += acc*dt } return (time:time, vel:vel) } func velocityVerlet(pos: Double, acc: Double, dt: Double) -> (time: Double, vel: Double) { var pos = pos var time, vel : Double vel = 0 time = 0 while (pos > 0) { time += dt pos += vel*dt + 0.5*acc * dt * dt vel += acc*dt } return (time:time, vel:vel) } func main() { let verletTime = verlet(pos: 5.0, acc: -10.0, dt: 0.01) print("Time for Verlet integration is: \(verletTime)") let stormer = stormerVerlet(pos: 5.0, acc: -10.0, dt: 0.01); print("Time for Stormer Verlet integration is: \(stormer.time)") print("Velocity for Stormer Verlet integration is: \(stormer.vel)") let velVerlet = velocityVerlet(pos: 5.0, acc: -10, dt: 0.01) print("Time for velocity Verlet integration is: \(velVerlet.time)") print("Velocity for velocity Verlet integration is: \(velVerlet.vel)") } main() SUBROUTINE verlet(pos, acc, dt, time) IMPLICIT NONE REAL(8), INTENT(INOUT) :: pos, acc, dt, time REAL(8) :: prev_pos, next_pos prev_pos = pos time = 0d0 DO IF (pos > 0d0) THEN time = time + dt next_pos = pos * 2d0 - prev_pos + acc * dt ** 2 prev_pos = pos pos = next_pos ELSE EXIT END IF END DO END SUBROUTINE verlet SUBROUTINE stormer_verlet(pos, acc, dt, time, vel) IMPLICIT NONE REAL(8), INTENT(INOUT) :: pos, acc, dt, time, vel REAL(8) :: prev_pos, next_pos prev_pos = pos time = 0d0 vel = 0d0 DO IF (pos > 0d0) THEN time = time + dt next_pos = pos * 2 - prev_pos + acc * dt ** 2 prev_pos = pos pos = next_pos vel = vel + acc * dt ELSE EXIT END IF END DO END SUBROUTINE stormer_verlet SUBROUTINE velocity_verlet(pos, acc, dt, time, vel) IMPLICIT NONE REAL(8), INTENT(INOUT) :: pos, acc, dt, time, vel time = 0d0 vel = 0d0 DO IF (pos > 0d0) THEN time = time + dt pos = pos + vel * dt + 0.5d0 * acc * dt ** 2 vel = vel + acc * dt ELSE EXIT END IF END DO END SUBROUTINE velocity_verlet PROGRAM verlet_integration IMPLICIT NONE REAL(8) :: pos,acc, dt, time, vel INTERFACE SUBROUTINE verlet(pos, acc, dt, time) REAL(8), INTENT(INOUT) :: pos, acc, dt, time REAL(8) :: prev_pos, next_pos END SUBROUTINE END INTERFACE INTERFACE SUBROUTINE stormer_verlet(pos, acc, dt, time, vel) REAL(8), INTENT(INOUT) :: pos, acc, dt, time, vel REAL(8) :: prev_pos, next_pos END SUBROUTINE END INTERFACE INTERFACE SUBROUTINE velocity_verlet(pos, acc, dt, time, vel) REAL(8), INTENT(INOUT) :: pos, acc, dt, time, vel REAL(8) :: prev_pos, next_pos END SUBROUTINE END INTERFACE pos = 5d0 acc = -10d0 dt = 0.01d0 ! Verlet CALL verlet(pos, acc, dt, time) WRITE(*,*) 'Time for Verlet integration: ', time ! stormer Verlet pos = 5d0 CALL stormer_verlet(pos, acc, dt, time, vel) WRITE(*,*) 'Time for Stormer-Verlet integration: ', time ! Velocity Verlet pos = 5d0 CALL velocity_verlet(pos, acc, dt, time, vel) WRITE(*,*) 'Time for Velocity-Verlet integration: ', time END PROGRAM verlet_integration def verlet(pos, acc, dt) prev_pos = pos time = 0 while pos > 0 do time += dt temp_pos = pos pos = pos*2 - prev_pos + acc * dt * dt prev_pos = temp_pos end return time end def stormer_verlet(pos, acc, dt) prev_pos = pos vel = 0 time = 0 while pos > 0 do time += dt temp_pos = pos pos = pos*2 - prev_pos + acc * dt * dt prev_pos = temp_pos vel += acc*dt end return time end def velocity_verlet(pos, acc, dt) vel = 0 time = 0 while pos > 0 do time += dt pos += vel*dt + 0.5*acc * dt * dt vel += acc*dt end return time end p verlet(5.0, -10, 0.01) p stormer_verlet(5.0, -10, 0.01) p velocity_verlet(5.0, -10, 0.01) package main import "fmt" func verlet(pos, acc, dt float64) (time float64) { prevPos := pos time = 0 for pos > 0 { time += dt nextPos := pos*2 - prevPos + acc*dt*dt prevPos, pos = pos, nextPos } return } func stormerVerlet(pos, acc, dt float64) (time, vel float64) { prevPos := pos time, vel = 0, 0 for pos > 0 { time += dt vel += acc * dt nextPos := pos*2 - prevPos + acc*dt*dt prevPos, pos = pos, nextPos } return } func velocityVerlet(pos, acc, dt float64) (time, vel float64) { time, vel = 0, 0 for pos > 0 { time += dt pos += vel*dt + .5*acc*dt*dt vel += acc * dt } return } func main() { time := verlet(5., -10., .01) fmt.Println("Verlet") fmt.Println("Time:", time) fmt.Println() time, vel := stormerVerlet(5., -10., .01) fmt.Println("Stormer-Verlet") fmt.Println("Time:", time) fmt.Println("Velocity:", vel) fmt.Println() time, vel = velocityVerlet(5., -10., .01) fmt.Println("Velocity Verlet") fmt.Println("Time:", time) fmt.Println("Velocity:", vel) } .intel_syntax noprefix .section .rodata zero: .double 0.0 two: .double 2.0 half: .double 0.5 verlet_fmt: .string "Time for Verlet integration is: %lf\n" stormer_fmt: .string "Time and Velocity for Stormer Verlet Integration is: %lf, %lf\n" velocity_fmt: .string "Time and Velocity for Velocity Verlet Integration is: %lf, %lf\n" pos: .double 5.0 acc: .double -10.0 dt: .double 0.01 .section .text .global main .extern printf # xmm0 - pos # xmm1 - acc # xmm2 - dt # RET xmm0 - time verlet: pxor xmm7, xmm7 # Holds 0 for comparisons pxor xmm3, xmm3 # Holds time value comisd xmm0, xmm7 # Check if pos is greater then 0.0 jbe verlet_return movsd xmm6, xmm1 # xmm6 = acc * dt * dt mulsd xmm6, xmm2 mulsd xmm6, xmm2 movsd xmm5, xmm0 # Holds previous position verlet_loop: addsd xmm3, xmm2 # Adding dt to time movsd xmm4, xmm0 # Hold old value of posistion addsd xmm0, xmm0 # Calculating new position subsd xmm0, xmm5 addsd xmm0, xmm6 movsd xmm5, xmm4 comisd xmm0, xmm7 # Check if position is greater then 0.0 ja verlet_loop verlet_return: movsd xmm0, xmm3 # Saving time value ret # xmm0 - pos # xmm1 - acc # xmm2 - dt # RET xmm0 - time # RET xmm1 - velocity stormer_verlet: pxor xmm7, xmm7 # Holds 0 for comparisons pxor xmm3, xmm3 # Holds time value comisd xmm0, xmm7 # Check if pos is greater then 0.0 jbe stormer_verlet_return movsd xmm6, xmm1 # xmm6 = acc * dt * dt mulsd xmm6, xmm2 mulsd xmm6, xmm2 movsd xmm5, xmm0 # Holds previous position stormer_verlet_loop: addsd xmm3, xmm2 # Adding dt to time movsd xmm4, xmm0 # Hold old value of posistion addsd xmm0, xmm0 # Calculating new position subsd xmm0, xmm5 addsd xmm0, xmm6 movsd xmm5, xmm4 comisd xmm0, xmm7 # Check if position is greater then 0.0 ja stormer_verlet_loop stormer_verlet_return: movsd xmm0, xmm3 # Saving time and velocity mulsd xmm3, xmm1 movsd xmm1, xmm3 ret # xmm0 - pos # xmm1 - acc # xmm2 - dt # RET xmm0 - time # RET xmm1 - velocity velocity_verlet: pxor xmm7, xmm7 # Holds 0 for comparisons pxor xmm3, xmm3 # Holds the velocity value pxor xmm4, xmm4 # Holds the time value comisd xmm0, xmm7 # Check if pos is greater then 0.0 jbe velocity_verlet_return movsd xmm5, half # xmm5 = 0.5 * dt * dt * acc mulsd xmm5, xmm2 mulsd xmm5, xmm2 mulsd xmm5, xmm1 velocity_verlet_loop: movsd xmm6, xmm3 # Move velocity into register mulsd xmm6, xmm2 # Calculate new position addsd xmm6, xmm5 addsd xmm0, xmm6 addsd xmm4, xmm2 # Incrementing time movsd xmm3, xmm4 # Updating velocity mulsd xmm3, xmm1 comisd xmm0, xmm7 ja velocity_verlet_loop velocity_verlet_return: movsd xmm0, xmm4 # Saving time and velocity movsd xmm1, xmm3 ret main: push rbp movsd xmm0, pos # Calling verlet movsd xmm1, acc movsd xmm2, dt call verlet mov rdi, OFFSET verlet_fmt # Print output mov rax, 1 call printf movsd xmm0, pos # Calling stormer_verlet movsd xmm1, acc movsd xmm2, dt call stormer_verlet mov rdi, OFFSET stormer_fmt # Print output mov rax, 1 call printf movsd xmm0, pos # Calling velocity_verlet movsd xmm1, acc movsd xmm2, dt call velocity_verlet mov rdi, OFFSET velocity_fmt # Print output mov rax, 1 call printf pop rbp ret data class VerletValues(val time: Double, val vel: Double) fun verlet(_pos: Double, acc: Double, dt: Double): Double { var pos = _pos // Since function parameter are val and can't be modified var prevPos = pos var time = 0.0 while (pos > 0) { time += dt val nextPos = pos * 2 - prevPos + acc * dt * dt prevPos = pos pos = nextPos } return time } fun stormerVerlet(_pos: Double, acc: Double, dt: Double): VerletValues { var pos = _pos var prevPos = pos var time = 0.0 var vel = 0.0 while (pos > 0) { time += dt val nextPos = pos * 2 - prevPos + acc * dt * dt prevPos = pos pos = nextPos vel += acc * dt } return VerletValues(time, vel) } fun velocityVerlet(_pos: Double, acc: Double, dt: Double): VerletValues { var pos = _pos var time = 0.0 var vel = 0.0 while (pos > 0) { time += dt pos += vel * dt + 0.5 * acc * dt * dt vel += acc * dt } return VerletValues(time, vel) } fun main(args: Array<String>) { val verletTime = verlet(5.0, -10.0, 0.01) println("Time for Verlet integration is:$verletTime")

val stormerVerlet = stormerVerlet(5.0, -10.0, 0.01)
println("Time for Stormer Verlet integration is: $stormerVerlet.time") println("Velocity for Stormer Verlet integration is:$stormerVerlet.vel")

val velocityVerlet = velocityVerlet(5.0, -10.0, 0.01)
println("Time for Velocity Verlet integration is: $velocityVerlet.time") println("Velocity for Velocity Verlet integration is:$velocityVerlet.vel")
}
proc verlet(pos_in, acc, dt: float): float =
var
pos: float = pos_in
prevPos: float = pos
time: float = 0.0
tempPos: float

while pos > 0.0:
time += dt
tempPos = pos
pos = pos * 2 - prevPos + acc * dt * dt
prevPos = tempPos

return time

proc stormerVerlet(pos_in, acc, dt: float): (float, float) =
var
pos: float = pos_in
prevPos: float = pos
time: float = 0.0
vel: float = 0.0
tempPos: float

while pos > 0.0:
time += dt
tempPos = pos
pos = pos * 2 - prevPos + acc * dt * dt
prevPos = tempPos

vel += acc * dt

return (time, vel)

proc velocityVerlet(pos_in, acc, dt: float): (float, float) =
var
pos: float = pos_in
time: float = 0.0
vel: float = 0.0

while pos > 0.0:
time += dt
pos += vel * dt + 0.5 * acc * dt * dt
vel += acc * dt

return (time, vel)

let timeV = verlet(5.0, -10.0, 0.01)
echo "Time for Verlet integration is: ", timeV

let (timeSV, velSV) = stormerVerlet(5.0, -10.0, 0.01)
echo "Time for Stormer Verlet integration is: ", timeSV
echo "Velocity for Stormer Verlet integration is: ", velSV

let (timeVV, velVV) = velocityVerlet(5.0, -10.0, 0.01)
echo "Time for velocity Verlet integration is: ", timeVV
echo "Velocity for velocity Verlet integration is: ", velVV
;;;; Verlet integration implementation in Common Lisp

(defun verlet (pos acc dt)
"Integrates Newton's equation for motion while pos > 0 using Verlet integration."
(loop
with prev-pos = pos
for time = 0 then (incf time dt)
while (> pos 0)
;; The starting speed is assumed to be zero.
do (psetf
pos (+ (* pos 2) (- prev-pos) (* acc dt dt))
prev-pos pos)
finally (return time)))

(defun stormer-verlet (pos acc dt)
"Integrates Newton's equation for motion while pos > 0 using the Stormer-Verlet method."
(loop
with prev-pos = pos
for time = 0 then (incf time dt)
for vel = 0 then (incf vel (* acc dt))
while (> pos 0)
;; Variables are changed simultaneously by 'psetf', so there's no need for a temporary variable.
do (psetf
pos (+ (* pos 2) (- prev-pos) (* acc dt dt))
prev-pos pos)
finally (return (list time vel))))

(defun velocity-verlet (pos acc dt)
"Integrates Newton's equation for motion while pos > 0 using the velocity in calculations."
(loop
for time = 0 then (incf time dt)
for vel = 0 then (incf vel (* acc dt))
for p = pos then (incf p (+ (* vel dt) (* 0.5 acc dt dt)))
while (> p 0)
finally (return (list time vel))))

(format T "Time for Verlet integration: ~d~%" (verlet 5 -10 0.01))

(defvar stormer-verlet-result (stormer-verlet 5 -10 0.01))
(format T "Time for Stormer Verlet integration is: ~d~%" (first stormer-verlet-result))
(format T "Velocity for Stormer Verlet integration is: ~d~%" (second stormer-verlet-result))

(defvar velocity-verlet-result (velocity-verlet 5 -10 0.01))
(format T "Time for velocity Verlet integration is: ~d~%" (first velocity-verlet-result))
(format T "Velocity for velocity Verlet integration is: ~d~%" (second velocity-verlet-result))

Code Examples

The code examples are licensed under the MIT license (found in LICENSE.md).

Text

The text of this chapter was written by James Schloss and is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. Pull Requests

After initial licensing (#560), the following pull requests have modified the text or graphics of this chapter:

• none