A velocity Verlet routine integrates a simple harmonic oscillator with , ,

PHYSICS · Multi-part · Fermi University Physics · Mathematical Methods for Physics · Computational Physics · Research software practice, verification and validation

Scientific Computing for Physicists, Ch Is the equation solved right? Testing, debugging and code verification

Problem

A velocity Verlet routine integrates a simple harmonic oscillator with

\omega = 2\,\mathrm{rad\,s^{-1}}

x(0) = 1\,\mathrm{m}

v(0) = 0\,\mathrm{m\,s^{-1}}

, whose exact solution is

x(t) = \cos(\omega t)

. ((a)) A unit test runs the routine with

h = 0.001\,\mathrm{s}

for

1000

steps to

t = 1\,\mathrm{s}

and compares

x(1)

with the exact value. State the exact value to six significant figures, and say whether the test passes a tolerance of

1\times10^{-6}\,\mathrm{m}

, given that the measured error at this step is

3.0\times10^{-7}\,\mathrm{m}

. ((b)) A property test runs the routine forward for

N = 2000

steps and then backward for

2000

steps with the same

h

, and checks that the state returns to its start. Each of the

2000

forward-backward step pairs contributes a rounding error of order

\epsilon_{\rm mach} = 2^{-52}

; take the property-test tolerance as

10N\epsilon_{\rm mach}/2

, relative to the amplitude

1\,\mathrm{m}

. Compute the tolerance, and say whether a measured return-difference of

9.7\times10^{-17}

(the larger of the position and velocity differences) passes it. ((c)) An order-of-accuracy test on the same problem gives maximum errors

e(0.02) = 1.21\times10^{-4}\,\mathrm{m}

e(0.01) = 3.03\times10^{-5}\,\mathrm{m}

and

e(0.005) = 7.58\times10^{-6}\,\mathrm{m}

. Compute the observed order

\hat{p}

from each consecutive pair using

\hat{p} = \ln[e(h)/e(h/2)]/\ln 2

, and say whether the code is verified against velocity Verlet's formal order

p_{\rm formal} = 2

Browse the problem library