Error Propagation of Verlet Algorithm

Error Propagation of Verlet Algorithm Binbin Li Civil and Environmental Engineering, UC Berkeley

1. Introduction Verlet Algorithm is a popular numerical integration method, especially in molecular dynamics (MD) simulation. The discussion of errors in Verlet Algorithm appears in almost all books about MD simulation, but many of their descriptions are conflicted, and none of them provides strictly mathematical, thus it confuses most of readers, like me. Here, I try to talk about the error analysis of Verlet Algorithm using mathematical derivations, but it is not guaranteed to be errorfree. 2. Verlet Algorithm Two Verlet Algorithms will be briefly discussed: the Original Verlet and the Velocity Verlet, and the updating procedures are as follows: Orignal Verlet: 𝑟(𝑡 + Δ𝑡) = 2𝑟(𝑡) − 𝑟(𝑡 − Δ𝑡) + 𝑎(𝑡)Δ𝑡 2 (1) 𝑣(𝑡) =

𝑟(𝑡+Δ𝑡)−𝑟(𝑡−Δ𝑡)

(2)

2Δ𝑡

Velocity Verlet: 1

𝑟(𝑡 + Δ𝑡) = 𝑟(𝑡) + 𝑣(𝑡)Δ𝑡 + 2 𝑎(𝑡)Δ𝑡 2 𝑣(𝑡 + Δ𝑡) = 𝑣(𝑡) +

Δ𝑡 2

[𝑎(𝑡) + 𝑎(𝑡 + Δ𝑡)]

(3) (4)

where 𝑟(𝑡), 𝑣(𝑡) and 𝑎(𝑡) are displacement, velocity and acceleration at 𝑡, respectively, and Δ𝑡 is the integration time step. In fact, these two Verlet Algorithms are mathematically equivalent, i.e. we could derive one from the other1. 3. Error Analysis The errors considered here include the truncation error (local error) and global error. The basic mathematical tool used is Taylor expansion. 3.1 Original Verlet First, let’s consider the Original Verlet. Taking the Taylor expansion of 𝑟(𝑡 + Δ𝑡) and 𝑟(𝑡 − Δ𝑡) at 𝑟(𝑡) gives 1

1

1

1

𝑟(𝑡 + Δ𝑡) = 𝑟(𝑡) + 𝑣(𝑡)Δ𝑡 + 2 𝑎(𝑡)Δ𝑡 2 + 6 𝑗(𝑡)Δ𝑡 3 + 𝑂(Δ𝑡 4 ) 𝑟(𝑡 − Δ𝑡) = 𝑟(𝑡) − 𝑣(𝑡)Δ𝑡 + 2 𝑎(𝑡)Δ𝑡 2 − 6 𝑗(𝑡)Δ𝑡 3 + 𝑂(Δ𝑡 4 )

1

http://www.physics.udel.edu/~bnikolic/teaching/phys660/numerical_ode/node5.html

(5) (6)

where 𝑗(𝑡) is the “jerk” at 𝑡 (i.e. the third derivative of displacement). Summarizing Eqn. (5) and (6) and doing some elementary operations yields 𝑟(𝑡 + Δ𝑡) = 2𝑟(𝑡) − 𝑟(𝑡 − Δ𝑡) + 𝑎(𝑡)Δ𝑡 2 + 𝑂(Δ𝑡 4 ) (7) Subtracting Eqn. (6) by (5) and doing some elementary operations gives 𝑣(𝑡) =

𝑟(𝑡+Δ𝑡)−𝑟(𝑡−Δ𝑡) 2Δ𝑡

+ 𝑂(Δ𝑡 2 )

(8)

By comparing Eqn. (1) and (7), (2) and (8), we can conclude that the truncation error of displacement and velocity are 𝑂(Δ𝑡 4 ) and (Δ𝑡 2 ) , respectively. Based on these truncation errors, the global error could be derived as follows2: First, we have 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + Δ𝑡)] = 𝑂(Δ𝑡 4 ) and 𝑟(𝑡 + 2Δ𝑡) = 2𝑟(𝑡 + Δ𝑡) − 𝑟(𝑡) + 𝑎(𝑡 + Δ𝑡)Δ𝑡 2 + 𝑂(Δ𝑡 4 ) Then, 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 2Δ𝑡)] = 2 × 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + Δ𝑡)] + 𝑂(Δ𝑡 4 ) = 3𝑂(Δ𝑡 4 ) 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 3Δ𝑡)] = 2 × 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 2Δ𝑡)] + 𝑂(Δ𝑡 4 ) = 6𝑂(Δ𝑡 4 ) 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 4Δ𝑡)] = 2 × 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 2Δ𝑡)] + 𝑂(Δ𝑡 4 ) = 10𝑂(Δ𝑡 4 ) and so on. By induction, we get 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 𝑛Δ𝑡)] =

𝑛(𝑛+1) 2

𝑂(Δ𝑡 4 )

Considering the global error in displacement between 𝑟(𝑡) and 𝑟(𝑡 + 𝑇) with 𝑇 = 𝑛Δ𝑡, it is clear that 1 𝑇

𝑇

𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 𝑇)] = 2 Δ𝑡 (Δ𝑡 + 1) 𝑂(Δ𝑡 4 ) = 𝑂(Δ𝑡 2 )

(9)

Using the same way, we can derive the global error for velocity: 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + Δ𝑡)] =

𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡+2Δ𝑡)]−𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡)]

𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 2Δ𝑡)] =

+ 𝑂(Δ𝑡 2 ) = 𝑂(Δ𝑡 3 ) + 𝑂(Δ𝑡 2 ) = 𝑂(Δ𝑡 2 )

2Δ𝑡 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡+3Δ𝑡)]−𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡+Δ𝑡)] 2Δ𝑡

+ 𝑂(Δ𝑡 2 ) = 𝑂(Δ𝑡 2 )

… 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 𝑛Δ𝑡)] = 𝑂(Δ𝑡 2 ) Therefore, 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 𝑇)] = 𝑂(Δ𝑡 2 ) (10) We see that the error of velocity does not increase with time. Conclusively, Original Verlet is a second-order integrator both for displacement and velocity. 3.2 Velocity Verlet The error analysis of Velocity Verlet is a little more complex than Original Verlet. Let’s first show the truncation error for velocity. Taking Taylor expansion of 𝑎(𝑡 +

Δ𝑡 2

) and 𝑎(𝑡 + Δ𝑡) at 𝑎(𝑡) gives

𝑎 (𝑡 +

Δ𝑡 2

1

) = 𝑎(𝑡) + 2 𝑗(𝑡)Δ𝑡 + 𝑂(Δ𝑡 2 )

𝑎(𝑡 + Δ𝑡) = 𝑎(𝑡) + 𝑗(𝑡)Δ𝑡 + 𝑂(Δ𝑡 2 ) Eqn. (11) minus half of Eqn. (12) gives 2

http://www.saylor.org/site/wp-content/uploads/2011/06/MA221-6.1.pdf

(11) (12)

𝑎 (𝑡 +

Δ𝑡

)= 2

𝑎(𝑡)+𝑎(𝑡+Δ𝑡) 2

+ 𝑂(Δ𝑡 2 )

Then, taking Taylor expansion of 𝑣(𝑡)and 𝑣(𝑡 + Δ𝑡) at 𝑣(𝑡 + 𝑣(𝑡) = 𝑣 (𝑡 +

Δ𝑡

) − 𝑎 (𝑡 + 2

𝑣(𝑡 + Δ𝑡) = 𝑣 (𝑡 +

Δ𝑡 Δ𝑡

) 2

Δ𝑡

2

) + 𝑎 (𝑡 + 2

Δ𝑡

1

+ 2 𝑗 (𝑡 +

Δ𝑡 Δ𝑡

) 2

2

2

) gives

Δ𝑡

1

(13)

Δ𝑡 2

) ( 2 ) + 𝑂(Δ𝑡 3 ) 2

+ 2 𝑗 (𝑡 +

Δ𝑡 2

Δ𝑡 2

) ( 2 ) + 𝑂(Δ𝑡 3 )

(14) (15)

Eqn. (15) minus Eqn. (14) gives 𝑣(𝑡 + Δ𝑡) = 𝑣(𝑡) + 𝑎 (𝑡 +

Δ𝑡 2

) Δ𝑡 + 𝑂(Δ𝑡 3 )

(16)

Substituting Eqn. (13) into Eqn. (16) turns out 𝑣(𝑡 + Δ𝑡) = 𝑣(𝑡) + [ = 𝑣(𝑡) +

𝑎(𝑡)+𝑎(𝑡+Δ𝑡)

Δ𝑡 2

2

+ 𝑂(Δ𝑡 2 )] Δ𝑡 + 𝑂(Δ𝑡 3 )

[𝑎(𝑡) + 𝑎(𝑡 + Δ𝑡)] + 𝑂(Δ𝑡 3 )

(17) 3 ).

By comparing Eqn. (4) and Eqn. (17), it easy to see the truncation error of velocity is 𝑂(Δ𝑡 For the truncation error of displacement, just taking Taylor expansion of 𝑟(𝑡 + Δ𝑡) at 𝑟(𝑡) gives the answer 1

𝑟(𝑡 + Δ𝑡) = 𝑟(𝑡) + 𝑣(𝑡)Δ𝑡 + 2 𝑎(𝑡)Δ𝑡 2 + 𝑂(Δ𝑡 3 )

(18)

Thus, the truncation error of displacement is also 𝑂(Δ𝑡 3 ). Now, let’s consider the global error of Velocity Verlet. Starting from 𝑣(𝑡 + Δ𝑡) = 𝑂(Δ𝑡 3 ) and (𝑡 + 2Δ𝑡) = 𝑣(𝑡 + Δ𝑡) +

Δ𝑡 2

[𝑎(𝑡 + Δ𝑡) + 𝑎(𝑡 + 2Δ𝑡)] + 𝑂(Δ𝑡 3 ) , we have

𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 2Δ𝑡)] = 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + Δ𝑡)] + 𝑂(Δ𝑡 3 ) = 2𝑂(Δ𝑡 3 ) 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 3Δ𝑡)] = 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 2Δ𝑡)] + 𝑂(Δ𝑡 3 ) = 3𝑂(Δ𝑡 3 ) …… 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 𝑛Δ𝑡)] = 𝑛𝑂(Δ𝑡 3 ) For the global error in velocity between 𝑣(𝑡) and 𝑣(𝑡 + 𝑇) with 𝑇 = 𝑛Δ𝑡, it is clear that 𝑇

𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 𝑇)] = Δ𝑡 𝑂(Δ𝑡 3 ) = 𝑂(Δ𝑡 2 )

(19)

About the global error in displacement, 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + Δ𝑡)] = 𝑂(Δ𝑡 3 ), 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + Δ𝑡)] = 1

𝑂(Δ𝑡 3 ) and 𝑟(𝑡 + 2Δ𝑡) = 𝑟(𝑡 + Δ𝑡) + 𝑣(𝑡 + Δ𝑡)Δ𝑡 + 2 𝑎(𝑡 + Δ𝑡)Δ𝑡 2 + 𝑂(Δ𝑡 3 ), then 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 2Δ𝑡)] = 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + Δ𝑡)] + 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + Δ𝑡)]Δ𝑡 + 𝑂(Δ𝑡 3 ) = 2𝑂(Δ𝑡 3 ) 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 3Δ𝑡)] = 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 2Δ𝑡)] + 𝑒𝑟𝑟𝑜𝑟[𝑣(𝑡 + 2Δ𝑡)]Δ𝑡 + 𝑂(Δ𝑡 3 ) = 3𝑂(Δ𝑡 3 ) …… 𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 𝑛Δ𝑡)] = 𝑛𝑂(Δ𝑡 3 ) Therefore, the global error in displacement between 𝑟(𝑡) and 𝑟(𝑡 + 𝑇) with 𝑇 = 𝑛Δ𝑡 is 𝑇

𝑒𝑟𝑟𝑜𝑟[𝑟(𝑡 + 𝑇)] = Δ𝑡 𝑂(Δ𝑡 3 ) = 𝑂(Δ𝑡 2 )

(20)

We see that the global error in Velocity Verlet are both 𝑂(Δ𝑡 2 ) for displacement and velocity, and it is the same as Original Verlet. Conclusion This paper is concluded by the follow table:

Order Truncation error Global error

Table 1 Error in Verlet Algorithm Original Verlet Velocity Verlet Displacement Velocity Displacement Velocity 𝑂(Δ𝑡 4 ) 𝑂(Δ𝑡 2 ) 𝑂(Δ𝑡 3 ) 𝑂(Δ𝑡 3 ) 2) 2) 2) 𝑂(Δ𝑡 𝑂(Δ𝑡 𝑂(Δ𝑡 𝑂(Δ𝑡 2 )