Numerical Error¶
Arithmetic with floating-point numbers is not exact, causing accumulation of numerical error.
In modern computers, the floating point is presented as $$x=(-1)^s\,m\,2^E,$$ where:
- $s$ is the sign bit,
- $m$ is the significand,
- $E$ is the binary exponent.
The largest and the smallest floating point number depends on the type. Most often we will
use python float, which needs 8bytes=64bits and can store numbers between 2.22507e-308 to
1.79769e+308. It is composed of roughly: 11-bits exponent, 52-bits mantissa, 1-bit sign.
More precisely:
- 1 bit stores the sign,
- 11 bits store a biased exponent ($-1022\le E\le 1023$).
- 52 bits store the fractional part of the significand,
- the leading binary digit is implicit.
The overflow error occurs if we want to store $x > 1.79769\; 10^{308}$ and underflow when
$x < 5\; 10^{−324}$. These numbers are sufficiently large/small that usually do not cause problems.
If they do, we should probably work with logarithms of the numbers (log representation).
The roundoff error is the hardest problem to avoid, which occurs when : $1+\epsilon == 1$.
For 64-bit float it occurs around (only!) $2\, 10^{−16}$.
import numpy as np
info = np.finfo(float)
print("machine epsilon:", info.eps, 2**-52) # 52 bit mantissa
print("largest value:", info.max, (2-2**-52)*2**(2**10-1) ) # 11-bit exponent, but 10 bit positive and 10 negative
print("smallest normal:", info.tiny, 2**-(2**10-2)) # just smallest exponent possible, but mantisa is one
print("smallest subnormal:", np.nextafter(0.0, 1.0), 2**-52 * 2**-(2**10-2)) # smallest mantisa & smallest exponent
machine epsilon: 2.220446049250313e-16 2.220446049250313e-16 largest value: 1.7976931348623157e+308 1.7976931348623157e+308 smallest normal: 2.2250738585072014e-308 2.2250738585072014e-308 smallest subnormal: 5e-324 5e-324
np.finfo(float).tiny is the smallest positive normal binary64 number, not the smallest representable positive number. Below tiny, IEEE-754 uses subnormal numbers, which provide gradual underflow down to approximately $5\; 10^{-324}$. Subnormal numbers have progressively fewer significant bits, so their relative precision deteriorates as zero is approached.
for i in range(1,10):
x=1+i*1e-16-1
print('1 + '+str(i)+'e-16 -1 = ', x)
1 + 1e-16 -1 = 0.0 1 + 2e-16 -1 = 2.220446049250313e-16 1 + 3e-16 -1 = 2.220446049250313e-16 1 + 4e-16 -1 = 4.440892098500626e-16 1 + 5e-16 -1 = 4.440892098500626e-16 1 + 6e-16 -1 = 6.661338147750939e-16 1 + 7e-16 -1 = 6.661338147750939e-16 1 + 8e-16 -1 = 8.881784197001252e-16 1 + 9e-16 -1 = 8.881784197001252e-16
Many numerical methods rely on iterations, which cause so called error accumulation, and just within a few steps the numerical error makes result useless.
Typical example is the three term linear recurrence relation, $$x_{n+1}= a_n x_n + b_n x_{n-1},$$ which can very quickly become unstable. Here $a_n$ and $b_n$ are some real or complex numbers. A three-term recurrence relation is related to the finite-difference form of a general second-order differential equation. It is commonly used to obtain special functions, such as Bessel functions, associated Legendre functions, and regular Coulomb wave functions.
For constant coefficient, $a_n=a$ and $b_n=b$, the exact solution can be found by the zeros of the characteristic polynomial. More specifically, we search for the solution in the form $$x_n = C r^n, $$ which gives $$r^{n+1} = a r^n + b r^{n-1}$$ which is solved by characteristic polynomial of the form $$r^2 - a r - b=0.$$ For an $n$-term recurrence relation, the characteristic polynomial has degree $n-1$. For a three-term recurrence relation, we therefore have two independent solutions $$r_{1,2} = \frac{1}{2}a \pm \frac{1}{2}\sqrt{a^2+ 4 b}. $$ The general solution is then of the form $$ x_n = C_1 r_1^n + C_2 r_2^n. $$ If $a^2+4b=0$, and the characteristic polynomial has only one solution, the solution of the three-term recurrence is instead: $$ x_n = C_1 r_1^n + C_2 n r_1^n.$$
Even when coefficients are not constant, the $n$-term homogeneous linear recurrence relation generally has an $n-1$-dimensional space of solutions. They can not be found analytically or exactly, but we know that they exist. However, only one solution out of $n-1$ is easy to follow, namely, the one that grows the fastest in the direction of recurrence.
For a three-term recurrence relation, we have two solutions: one is typically growing faster than the other. If $|r_1/r_2|>1$ we could follow $r_1^n$ with upward recursion, but not $r_2$. Because even if we start with $r_2$ and tiny admixture of $r_1$, we will end up with admixture of $r_1$ of the form $x \approx r_2^n+\epsilon r_1^n$, which can be very different from $r_2^n$ even for small $\epsilon$.
But if we iterate in the opposite direction the other solution appears larger, and we could follow $r_2^n$. Mathematically, if we start at large $N$ and iterate down to $N-n$ we get $x \approx r_2^k + \epsilon (r_2/r_1)^N r_1^k = r_2^k(1+\epsilon (r_2/r_1)^{N-k})$. Thus, the relative contamination is $\epsilon\left(\frac{r_2}{r_1}\right)^{N-k}$.
Hence, for a three-term recurrence relation, this instability can be easily fixed by so-called Miller’s algorithm. If iteration in one direction is unstable, it is likely stable in the opposite direction.
Examples of three-term recurrence relations include:
Spherical Bessel Functions $$ j_{l+1}(x) = \frac{2l+1}{x} j_l(x) - j_{l-1}(x)$$ Ordinary Bessel Functions $$ J_{l+1}(x) = \frac{2l}{x} J_l(x) - J_{l-1}(x)$$ Modified Spherical Bessel Functions $$ j_{l+1}(x) = -\frac{2l+1}{x} j_l(x) + j_{l-1}(x)$$ Legendre Polynomials $$(n+1) P_{n+1}(x) = (2n+1) x P_n(x)-n P_{n-1}(x)$$ Associated Legendre polynomials $$(l-m+1) P^m_{l+1}(x)=(2l+1) x P^m_l(x) - (l+m)P^m_{l-1}$$ $$2 m x P^m_l(x) = -\sqrt{1-x^2}\left[ P^{m+1}_l(x)+(l+m)(l-m+1)P^{m-1}_l(x)\right]$$ in terms of them the spherical harmonics are given $$Y_{l,m}(\theta,\phi)=\sqrt{\frac{(2l+1)(l-m)!}{4\pi(l+m)!}}P^m_l(\cos\theta)e^{i m\phi},$$ which are used to solve the 3D partial differential equation of the form $$\nabla^2 \psi(\theta,\phi) + \lambda\; \psi(\theta,\phi)=0.$$ Confluent Hypergeometric Series $$(b-n) M_{n-1} + (2n-b+x)M_n - n M_{n+1}=0$$ with solution $M_n=M(n,b; x)$.
Spherical Bessel functions¶
The spherical Bessel functions satisfy the three-term recurrence relation:
$$ j_{l+1}(x)=\frac{2l+1}{x} j_l(x)-j_{l-1}(x)$$
with the initial condition $$j_0(x)=\frac{\sin x }{x}$$ and $$j_1(x)=\frac{\sin x}{x^2}-\frac{\cos x}{x}$$
This recurrence relation should be sufficient to compute $j_l(x)$ for any $l$ and any $x$. However, as we will see below, the numeric iteration becomes extremely unstable for $x \lesssim l$.
Any three-term linear recurrence relation must have exactly two solutions. In case of Bessel functions $j_l(x)$, we also have Neumann functions $n_l(x)$, which satisfy exactly the same recurrence relation, but have different initial conditions. If $n_l(x)$ is much larger than $j_l(x)$ the iteration will very quickly start to follow $n_l(x)$ even though we started with $j_l(x)$. As is well known, $n_l(x)$ diverges at small $x$, while $j_l(x)\propto x^l$ becomes very small. Obviously small $x$ and large $l$ are challenging for this algorithm.
To refresh our memory, spherical Bessel functions appear in wave equations with spherical geometry. Examples include the radial Schrödinger equation for a particle in a central potential, electromagnetic scattering by spherical objects, acoustic waves, and heat conduction in spherical geometry. Problems with cylindrical symmetry instead lead to ordinary Bessel functions $J_m(x)$ and $Y_m(x)$.
The spherical bessel functions are the solution of the following differential equation:
$$\left[-\frac{1}{2}\frac{d^2}{dr^2}+\frac{l(l+1)}{2 r^2}\right]\left(r j_l(k r)\right) = \frac{k^2}{2} \left(r j_l(k r)\right)$$
We will first iterate the recurrence relation for $j_l(x)$ starting at $l=0$. This is called upward recurrence. If this process proves to be unstable, we will demonstrate that downward recurrence becomes stable and serves as the solution to this problem.
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import spherical_jn
xs = np.linspace(0.1, 20, 200)
vals = [spherical_jn(0, xs),spherical_jn(1, xs),spherical_jn(2, xs)] # j_2(x) for many x
plt.plot(xs, vals[0], label='$j_0$')
plt.plot(xs, vals[1], label='$j_1$')
plt.plot(xs, vals[2], label='$j_2$')
plt.legend(loc='best')
plt.show()
Upward recurrence¶
We will evaluate bessel upward recurrence using the formula
\begin{equation} j_{l+1}(x) = \frac{2l+1}{x} j_l - j_{l-1} \end{equation}
import numpy as np
from scipy import special
def bessel_upward(l,x):
"returns array of j_i from i=0 to i=l, including l"
res = np.zeros(l+1)
if abs(x)<1e-30: # first take care of the special case, which is numerically hard.
res[0]=1.
return res
j0 = np.sin(x)/x
res[0]=j0
if l==0: return res
j1 = j0/x - np.cos(x)/x
res[1] = j1
for i in range(1,l):
j2 = (2*i+1)/x*j1 - j0 # (j2,j1,j0)==(j_{i+1},j_i,j_{i-1})
res[i+1]=j2 # store j_{i+1}
j0,j1 = j1,j2 # prepare for the next step (j_{i-1},j_i) <- (j_i,j_{i+1})
return res
help(special.spherical_jn)
from scipy.special import spherical_jn
l=10
x=0.1
dat0 = bessel_upward(l,x)
dat1 = spherical_jn(range(l+1),x)
diff = dat0-dat1
print("%4s %16s %16s %16s %16s" % ('l','exact','upward','abs-error','rel-error'))
for i in range(len(dat0)):
print("%4d %16.13f %16.13f %16.13f %16.8g" % (i,dat1[i],dat0[i],diff[i], diff[i]/dat1[i]))
l exact upward abs-error rel-error 0 0.9983341664683 0.9983341664683 0.0000000000000 0 1 0.0333000119026 0.0333000119026 0.0000000000000 3.7507521e-15 2 0.0006661906084 0.0006661906084 0.0000000000000 6.4359747e-12 3 0.0000095185197 0.0000095185199 0.0000000000002 2.2510964e-08 4 0.0000001057720 0.0000001057870 0.0000000000150 0.00014176421 5 0.0000000009616 0.0000000023109 0.0000000013493 1.4031448 6 0.0000000000074 0.0000001484162 0.0000001484088 20061.914 7 0.0000000000000 0.0000192917991 0.0000192917990 3.9116462e+08 8 0.0000000000000 0.0028936214469 0.0028936214469 9.9738775e+12 9 0.0000000000000 0.4918963541799 0.4918963541799 3.2213554e+17 10 0.0000000000000 93.4574136727252 93.4574136727252 1.2852544e+22
What happens with upward recurrence?
x=0.1
j9 = spherical_jn(9,x)
j10 = spherical_jn(10,x)
j11 = spherical_jn(11,x)
j11a = (2*10+1)/x * j10 - j9 # upward
j9a = (2*10+1)/x * j10 - j11 # downward
print('part1 upward or downward=', (2*10+1)/x*j10)
print('part2 upward=', -j9)
print('part2 downward=', -j11)
print('j9=', j9a, 'j9_exact=', j9, 'abs-err=', j9a-j9, 'rel-err=', (j9a-j9)/j9)
print('j11=', j11a, 'j11_exact=',j11, 'abs-err=', j11a-j11, 'rel-err=', (j11a-j11)/j11)
part1 upward or downward= 1.5270173093098784e-18 part2 upward= -1.5269856934948229e-18 part2 downward= -3.16158150515107e-23 j9= 1.526985693494827e-18 j9_exact= 1.5269856934948229e-18 abs-err= 4.044452883213195e-33 rel-err= 2.648651457864433e-15 j11= 3.161581505555956e-23 j11_exact= 3.16158150515107e-23 abs-err= 4.048855109557025e-33 rel-err= 1.280642331364966e-10
Downward recurrence : Miller's algorithm¶
Now we will use recurrence:
\begin{eqnarray} j_{l-1} = (2l+1)/x j_l - j_{l+1} \end{eqnarray}
Because the recurrence relation is linear and homogeneous, we are allowed to multiply all accumulated values $j_l$ with an arbitrary constant, and they will still represent solution of the same recurrence relation. We can use this to rescale all values $j_l$. We know $j_0=\frac{\sin(x)}{x}$, hence we can use this to rescale all $j_l$ in the series.
We can start with $j_{l_{max}+1}=0$ and $j_{l_{max}}=1$ as we know that the values of $j_l$ fall off quickly, and we are allowed to make substantial error at large enough $l$.
def bessel_downward(l,x):
"downward recursion"
if abs(x)<1e-20: # again take care of the special case where the algorithm would be unstable.
res = np.zeros(l+1)
res[0]=1
return res
extra = int(np.sqrt(10*l)) # this is a reasonable extra number of steps from NRC
lstart = l + extra
true_j0 = np.sin(x)/x
values = np.zeros(lstart + 2)
values[lstart + 1] = 0.0
values[lstart] = 1.0
for i in range(lstart, 0, -1):
values[i - 1] = (2*i + 1) / x * values[i] - values[i + 1]
# The overall normalization is arbitrary, so rescaling is allowed.
if abs(values[i - 1]) > 1e100:
values[i - 1:lstart + 2] *= 1e-100
values *= true_j0 / values[0]
return values[:l + 1]
l=10
x=0.1
dat0 = bessel_downward(l,x)
dat1 = spherical_jn(range(l+1),x)
diff=dat0-dat1
print("%4s %16s %16s %16s %16s" % ('l','exact','downward','abs-error','rel-error'))
for i in range(len(dat0)):
print("%4d %16.13f %16.13f %16.13f %16.8g" % (i,dat1[i],dat0[i],diff[i], diff[i]/dat1[i]))
l exact downward abs-error rel-error 0 0.9983341664683 0.9983341664683 0.0000000000000 0 1 0.0333000119026 0.0333000119026 -0.0000000000000 -6.2512535e-16 2 0.0006661906084 0.0006661906084 -0.0000000000000 -9.7647925e-16 3 0.0000095185197 0.0000095185197 -0.0000000000000 -1.4238062e-15 4 0.0000001057720 0.0000001057720 0.0000000000000 6.2563287e-16 5 0.0000000009616 0.0000000009616 -0.0000000000000 -2.1504626e-16 6 0.0000000000074 0.0000000000074 0.0000000000000 8.7358062e-16 7 0.0000000000000 0.0000000000000 -0.0000000000000 -1.5355307e-15 8 0.0000000000000 0.0000000000000 -0.0000000000000 -1.1895996e-15 9 0.0000000000000 0.0000000000000 -0.0000000000000 -1.8918939e-15 10 0.0000000000000 0.0000000000000 -0.0000000000000 -4.3453538e-15
def bessel_j(l,x):
"combines upward and downward recursion"
if l<=x : return bessel_upward(l,x)
# Upward recurrence is generally stable when l is not substantially larger than x.
# For l > x, use downward recurrence for the higher orders.
lcritical = int(x)
if lcritical<=0 : return bessel_downward(l,x) # for very small x we only need downward recursion
_ju_ = bessel_upward(lcritical-1,x) # upward works
_jd_ = bessel_downward(l,x)
return np.hstack( (_ju_, _jd_[lcritical:]) )
l=20
x=10
dat0 = bessel_upward(l,x)
dat1 = bessel_downward(l,x)
dat2 = bessel_j(l,x)
date = special.spherical_jn(range(l+1),x)
diff = date-dat2
#print('difference=', date-dat2)
#print('updard-diff=', date-dat0)
#print('down-diff=', date-dat1)
print("%4s %16s %16s %16s %16s" % ('l','exact','combined','abs-error','rel-error'))
for i in range(len(dat0)):
print("%4d %16.13f %16.13f %16.13f %16.8g" % (i, date[i], dat2[i], diff[i], diff[i]/date[i]))
l exact combined abs-error rel-error 0 -0.0544021110889 -0.0544021110889 0.0000000000000 -0 1 0.0784669417988 0.0784669417988 0.0000000000000 0 2 0.0779421936286 0.0779421936286 0.0000000000000 0 3 -0.0394958449845 -0.0394958449845 0.0000000000000 -0 4 -0.1055892851177 -0.1055892851177 0.0000000000000 -0 5 -0.0555345116215 -0.0555345116215 0.0000000000000 -2.4989484e-16 6 0.0445013223341 0.0445013223341 0.0000000000000 3.1185113e-16 7 0.1133862306558 0.1133862306558 0.0000000000000 0 8 0.1255780236496 0.1255780236496 -0.0000000000000 -2.2102255e-16 9 0.1000964095485 0.1000964095485 -0.0000000000000 -5.5457685e-16 10 0.0646051544926 0.0646051544926 -0.0000000000000 -4.7258045e-15 11 0.0355744148859 0.0355744148859 -0.0000000000000 -5.0713762e-15 12 0.0172159997450 0.0172159997450 -0.0000000000000 -4.8365897e-15 13 0.0074655844766 0.0074655844766 -0.0000000000000 -4.9957984e-15 14 0.0029410783418 0.0029410783418 -0.0000000000000 -5.1609745e-15 15 0.0010635427146 0.0010635427146 -0.0000000000000 -5.3010107e-15 16 0.0003559040735 0.0003559040735 -0.0000000000000 -5.6357153e-15 17 0.0001109407280 0.0001109407280 -0.0000000000000 -5.3750431e-15 18 0.0000323884744 0.0000323884744 -0.0000000000000 -5.4396774e-15 19 0.0000088966273 0.0000088966273 -0.0000000000000 -5.7124993e-15 20 0.0000023083720 0.0000023083720 -0.0000000000000 -5.8710327e-15
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
def scaled_error(numerical, exact, floor=1e-15):
return np.abs(numerical - exact) / np.maximum(np.abs(exact), floor)
l=15
x = np.linspace(1e-6,50,100)
dat0 = np.array([bessel_upward(l,t) for t in x])
dat1 = np.array([bessel_downward(l,t) for t in x])
dat2 = np.array([bessel_j(l,t) for t in x])
date = np.array([spherical_jn(range(l+1),t) for t in x])
for i in range(5,l,3):
plt.semilogy(x, scaled_error(dat2[:, i], date[:, i]), label=fr"$l={i}$")
plt.title('combination of up and down recursion')
plt.ylim([1e-20,1e-11])
plt.legend(loc='best')
plt.show()
for i in range(5,l,3):
plt.semilogy(x, scaled_error(dat1[:,i],date[:,i]), label=fr"$l={i}$")
plt.title('downward recursion')
plt.legend(loc='best')
plt.show()
for i in range(5,l,3):
plt.semilogy(x, scaled_error(dat0[:,i],date[:,i]), label=fr"$l={i}$")
plt.title('upward recursion')
plt.legend(loc='best')
plt.show()
Propagation of Floating-Point Underflow/Overflow¶
It's a different problem related to finite exponents of a floating point numbers. $\approx 2 10^{308}$
The common trick is to use logarithms throughout the algorithm.
Example: compute $n!$ for numbers $>200$. $100!\approx 9.3\times 10^{157}$ and $200!\approx 7.9 \times 10^{374}$.
Stable calculation requires one to use $log(n!)$ and only evaluate result when combined with other factors to obtain finite value. Typically $n!$ cancels with something else in the expression. To evaluate it, we do
\begin{equation} \log(n!) = \log(1\times 2 \times 3 \times \cdots \times n) = \sum_{k=2}^{n} \log(k) \end{equation}
import numpy as np
def log_factorial(n):
return np.sum(np.log(np.arange(2,n+1)))
def factorial(n):
fct=1.0 # Force binary64 arithmetic to demonstrate overflow.
for i in range(2,n+1):
fct *= i
return fct
np.log(factorial(171)+0.0) , log_factorial(171)
(np.float64(inf), np.float64(711.7147258022899))
- The problem can occur when computing the ratios:
$$ R = \frac{a_1 a_2 \cdots a_N}{b_1 b_2 \cdots b_M} $$ with an obvious solution $$ \log(R) = \sum_i \log(a_i)-\sum_j \log(b_j)$$
- It can also occur with matrices
$$ R = (A^1 A^2 \cdots A^N)({ B^1 B^2 \cdots B^M})^{-1}$$
A possible trick here is to rewrite the matrices in terms of a common exponent:
$${A^i}_{mn} = e^{a_i} c^i_{mn}$$ $${B^i}_{mn} = e^{b_i} d^i_{mn}$$ and largest $c_{mn}$ are of the order of unity, and consequently some $c_{mn}$ are very small.
We have $$R = e^{(\sum_i a_i-\sum_j b_j)} (c^1 c^2 \cdots c^N)(d^1 d^2 \cdots d^M)^{-1}$$
- Many times it occurs when dividing by the fermi or bose function:
$$ R = g(x)/f(x) $$ where $g(x)\approx e^{-x}$ and $f(x) = \frac{1}{1+e^x}$
One can reformulate the problem in terms of $\log(f(x))$, which has nice property $\log(f(x))=\log(f(-x))-x$.
Long-time integration: Runge–Kutta versus velocity Verlet¶
Roundoff is not the only error that accumulates in an iterative calculation. When a differential equation is replaced by finite time steps, each step also introduces a discretization error. A method can be very accurate over a short interval and nevertheless develop a systematic error during a sufficiently long simulation.
We will first implement both integrators for generic differential equations. We will then apply them to a nonlinear pendulum, whose conserved energy gives us a sensitive test of accumulated numerical error.
Generic fixed-step fourth-order Runge–Kutta¶
Consider a general first-order system
\begin{equation} \dot{\mathbf y}=f(t,\mathbf y). \end{equation}
The classical RK4 update is
\begin{align} k_1 &= f(t_n,\mathbf y_n),\\ k_2 &= f(t_n+\tfrac{1}{2}\Delta t,\mathbf y_n+\tfrac{1}{2}\Delta t\,k_1),\\ k_3 &= f(t_n+\tfrac{1}{2}\Delta t,\mathbf y_n+\tfrac{1}{2}\Delta t\,k_2),\\ k_4 &= f(t_n+\Delta t,\mathbf y_n+\Delta t\,k_3),\\ \mathbf y_{n+1} &= \mathbf y_n+\frac{\Delta t}{6}(k_1+2k_2+2k_3+k_4). \end{align}
RK4 has global error of order $O(\!\Delta t^4)$. The implementation below works for a scalar equation or a system of equations and allows nonuniform time points.
def rk4(derivative, y0, times):
"""Solve dy/dt = derivative(t, y) using fixed-step classical RK4."""
times = np.asarray(times, dtype=float)
y0 = np.asarray(y0, dtype=float)
solution = np.empty((len(times),) + y0.shape)
solution[0] = y0
for n in range(len(times)-1):
t = times[n]
dt = times[n+1] - t
y = solution[n]
k1 = derivative(t, y)
k2 = derivative(t + 0.5*dt, y + 0.5*dt*k1)
k3 = derivative(t + 0.5*dt, y + 0.5*dt*k2)
k4 = derivative(t + dt, y + dt*k3)
solution[n+1] = y + dt*(k1 + 2*k2 + 2*k3 + k4)/6
return solution
Generic velocity Verlet¶
Velocity Verlet applies to a second-order equation of the form
\begin{equation} \ddot{\mathbf q}=\mathbf a(t,\mathbf q), \end{equation}
where the acceleration does not depend on velocity. Its update is
\begin{align} \mathbf v_{n+1/2} &= \mathbf v_n+\frac{\Delta t}{2}\mathbf a(t_n,\mathbf q_n),\\ \mathbf q_{n+1} &= \mathbf q_n+\Delta t\,\mathbf v_{n+1/2},\\ \mathbf v_{n+1} &= \mathbf v_{n+1/2} +\frac{\Delta t}{2}\mathbf a(t_{n+1},\mathbf q_{n+1}). \end{align}
The method has global error of order $O(\!\Delta t^2)$. For an autonomous conservative system, it is symplectic and time reversible. It is therefore especially useful for long-time Hamiltonian dynamics.
This standard form does not directly apply when acceleration depends on velocity, as it does for friction or the magnetic part of the Lorentz force.
def velocity_verlet(acceleration, q0, v0, times):
"""Solve d²q/dt² = acceleration(t, q) using velocity Verlet."""
times = np.asarray(times, dtype=float)
q0 = np.asarray(q0, dtype=float)
v0 = np.asarray(v0, dtype=float)
position = np.empty((len(times),) + q0.shape)
velocity = np.empty((len(times),) + v0.shape)
position[0] = q0
velocity[0] = v0
for n in range(len(times)-1):
t = times[n]
dt = times[n+1] - t
a = acceleration(t, position[n])
v_half = velocity[n] + 0.5*dt*a
position[n+1] = position[n] + dt*v_half
a_new = acceleration(times[n+1], position[n+1])
velocity[n+1] = v_half + 0.5*dt*a_new
return position, velocity
Application: the nonlinear pendulum¶
For a pendulum of length $L$ in a uniform gravitational field,
\begin{equation} \ddot{\theta}=-\frac{g}{L}\sin\theta. \end{equation}
Taking $g/L=1$, the conserved energy per unit moment of inertia is
\begin{equation} E=\frac{1}{2}\dot{\theta}^{,2}+1-\cos\theta. \end{equation}
For RK4 we rewrite the equation as the first-order system
\begin{equation} \frac{d}{dt} \begin{pmatrix}\theta\\\omega\end{pmatrix} = \begin{pmatrix}\omega\\-\sin\theta\end{pmatrix}. \end{equation}
Velocity Verlet can use the second-order equation directly, with $a(\theta)=-\sin\theta$.
def pendulum_derivative(t, state):
theta, angular_velocity = state
return np.array([angular_velocity, -np.sin(theta)])
def pendulum_acceleration(t, theta):
return -np.sin(theta)
def pendulum_energy(theta, angular_velocity):
return 0.5*angular_velocity**2 + 1-np.cos(theta)
Short-time accuracy¶
We use the same step size for RK4 and velocity Verlet. A high-accuracy adaptive DOP853 calculation from SciPy supplies a reference solution. Because RK4 is fourth order while velocity Verlet is second order, RK4 should be substantially more accurate over this short interval.
from scipy.integrate import solve_ivp
theta0, angular_velocity0 = 2.0, 0.0
dt = 0.3
nsteps_short = int(20/dt)
t_short = dt*np.arange(nsteps_short+1)
# first solve with RK4
state_rk4 = rk4(pendulum_derivative, [theta0, angular_velocity0], t_short)
theta_rk4, angular_velocity_rk4 = state_rk4.T
# solve with verlet
theta_verlet, angular_velocity_verlet = velocity_verlet(pendulum_acceleration, theta0, angular_velocity0, t_short)
# also use scipy best method
reference = solve_ivp(
pendulum_derivative,
(t_short[0], t_short[-1]),
[theta0, angular_velocity0],
method='DOP853', t_eval=t_short,
rtol=1e-12, atol=1e-14
)
theta_reference = reference.y[0]
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(t_short, theta_reference, 'k-', label='reference')
axes[0].plot(t_short, theta_rk4, '--', label='RK4')
axes[0].plot(t_short, theta_verlet, ':', label='velocity Verlet')
axes[0].set_xlabel('time')
axes[0].set_ylabel(r'$\theta(t)$')
axes[0].set_title('Short-time trajectory')
axes[0].legend()
axes[1].semilogy(
t_short, np.maximum(np.abs(theta_rk4-theta_reference), 1e-16),
label='RK4'
)
axes[1].semilogy(
t_short, np.maximum(np.abs(theta_verlet-theta_reference), 1e-16),
label='velocity Verlet'
)
axes[1].set_xlabel('time')
axes[1].set_ylabel(r'$|\theta_{\rm numerical}-\theta_{\rm reference}|$')
axes[1].set_title('Angle error')
axes[1].legend()
plt.tight_layout()
plt.show()
Long-time energy behavior¶
We now integrate the pendulum for many oscillations. To expose the accumulated behavior clearly, we deliberately retain the rather large step $\Delta t=0.3$. Both algorithms use exactly the same time points.
RK4 initially gives a much more accurate trajectory, but it is not symplectic, and its small energy error develops a systematic drift. Velocity Verlet has a larger short-time trajectory error, yet its energy error remains bounded and oscillatory.
The conclusion is not that velocity Verlet is universally more accurate than RK4. Rather, high local accuracy and faithful long-time conservation are different numerical properties.
nsteps_long = 100_000
t_long = dt*np.arange(nsteps_long+1)
# solve with RK4
state_rk4 = rk4(pendulum_derivative, [theta0, angular_velocity0], t_long)
theta_rk4, angular_velocity_rk4 = state_rk4.T
# solve with Verlet
theta_verlet, angular_velocity_verlet = velocity_verlet(pendulum_acceleration, theta0, angular_velocity0, t_long)
# get energies for both
energy_rk4 = pendulum_energy(theta_rk4, angular_velocity_rk4)
energy_verlet = pendulum_energy(theta_verlet, angular_velocity_verlet)
energy_initial = pendulum_energy(theta0, angular_velocity0)
plt.figure(figsize=(10, 5))
plt.plot(t_long, (energy_rk4-energy_initial)/energy_initial, label='RK4')
plt.plot(t_long, (energy_verlet-energy_initial)/energy_initial, label='velocity Verlet')
plt.axhline(0.0, color='black', linewidth=0.7)
plt.xlabel('time')
plt.ylabel(r'$[E(t)-E(0)]/E(0)$')
plt.title('Accumulation of pendulum energy error over long times')
plt.legend()
plt.tight_layout()
plt.show()
At the same step size, RK4 is much more accurate over short times. The distinction is in the long-time structure of the error:
- RK4 has high local accuracy, but it does not preserve Hamiltonian phase-space geometry. A small systematic energy drift can eventually accumulate.
- Velocity Verlet has lower local accuracy, but its symplectic structure keeps the energy error bounded over long integrations of this regular conservative system.
Velocity Verlet does not conserve the exact pendulum energy at every step. It nearly conserves a modified Hamiltonian, which is why the measured energy oscillates. This behavior assumes a stable time step; no numerical method is stable for an arbitrarily large $\Delta t$.
In actual floating-point arithmetic, roundoff can also accumulate during an extremely long calculation. Symplecticity controls the systematic discretization error; it does not eliminate roundoff.
SciPy provides adaptive Runge–Kutta solvers through scipy.integrate.solve_ivp. We implemented fixed-step RK4 here because using identical time steps makes the comparison with velocity Verlet direct.
Homework¶
Implement the same forward and backward recursion for normal Bessel functions $J_n(x)$. Hence instead of using $j_0(x)=\frac{\sin(x)}{x}$ we will use scipy library for $J_0(x)$ and $J_1(x)$: special.jn(n, x).
You can check your results by comparing to the values obtained by special.jn(n, x).
Note that the recurence relation for ordinary Bessel function is slightly different:
$$ J_{l+1}(x) = \frac{2l}{x} J_l(x) - J_{l-1}(x)$$
hence $j_l(x) \propto J_{l+1/2}(x)$
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import jn
x = np.linspace(0.01,10,50)
plt.plot(x, jn(0,x), label=r'$J_0(x)$')
plt.plot(x, np.sin(x)/x, label=r'$j_0(x)$')
plt.legend(loc='best')
print( jn(0,x) )
[ 0.999975 0.98859675 0.95684376 0.90570153 0.83675353 0.75212657 0.65441767 0.54660502 0.43194604 0.31386663 0.19584514 0.08129554 -0.02654622 -0.12472852 -0.21067941 -0.28228758 -0.33796627 -0.37669788 -0.398058 -0.4022181 -0.3899271 -0.3624727 -0.32162411 -0.26955854 -0.2087743 -0.14199384 -0.07206054 -0.00183295 0.06591944 0.12861516 0.18395084 0.22997919 0.26517228 0.2884678 0.29929707 0.29759385 0.2837841 0.25875736 0.22382098 0.18063955 0.13116179 0.07753805 0.02203174 -0.03307184 -0.08555575 -0.13335715 -0.17464548 -0.20788998 -0.23191422 -0.24593576]
Old Homework:¶
We want to compute the series of integrals, defined by $$K_n(z,\alpha) = \int_0^1 dx \frac{x^n}{z+\alpha x}$$ when $n=0,1,...n_{max}=10$.
These occur, for example, when we evaluate matrix elements of the following correlation function $$<P_n(x)|G(x)|P_m(x)>,$$ where $P_n(x)$ is Legendre polynomial, and $G(x)=\frac{1}{z-x}$ is the Green's function.
We can derive recurrence relation by noting $$K_{n+1}(z,\alpha) = \int_0^1 dx \frac{x^n(x+z/\alpha) - x^n z/\alpha}{z+\alpha x}= \frac{1}{\alpha} \int_0^1 dx x^n -\frac{z}{\alpha}\int_0^1 dx \frac{x^n}{z+\alpha x}$$ which gives $$K_{n+1} = \frac{1}{\alpha(n+1)}-\frac{z}{\alpha} K_n$$
We can also calculate $K_0$: $$K_0=\frac{1}{\alpha} \log(1+\alpha/z)$$
Starting from $K_0$ you can compute $K_n$ up to $n_{max}$ using recurrence. This works quite well for $|\alpha/z|\gtrsim 1$.
Choosing $z$ and $\alpha$ so that $|\alpha/z|\ll 1$ (for example $\alpha/z=10^{-4}$) verify that upward recurrence is unstable.
Implement downward recurrence for $\alpha/z<1/2$. Since recurrence relation is not homogeneous, we can not start with arbitrary value and later normalize result. We thus need to start with very accurate value of $K_{n_{max}}$. Use power series expansion for $K_{n_{max}}$ in powers of $(\alpha/z)^k$, and evaluate as many terms as needed to achieve desired accuracy (say $10^{-12}$).
We can derive the power series in $\alpha/z$ by Taylor expansion: $$K_n(z,\alpha) = \frac{1}{z}\int_0^1 dx \frac{x^n}{1+\alpha/z x}= \frac{1}{z} \sum_{k=0}^{\infty}(-\alpha/z)^k \int_0^1 dx x^{n+k}= \frac{1}{z} \sum_{k=0}^{\infty}\frac{(-\alpha/z)^k}{n+k+1}$$ For $n$ large and $|\alpha/z|<1/2$ this series converges rapidly.
The first order recurrence relation can be solved analytically by the algorithm described in wikipedia:Recurrence relation
The solution is:
$$K_n = \left(-\frac{z}{\alpha}\right)^n\left[K_0 + \frac{1}{\alpha}\sum_{k=1}^{n} \frac{\left(-\frac{\alpha}{z}\right)^k}{k}\right]$$
By noting that Taylor expansion of $\log(1+x)=-\sum_{k=1}^{\infty} \frac{(-x)^k}{k}$, we see that $$K_0=-\frac{1}{\alpha}\sum_{k=1}^\infty \frac{\left(-\frac{\alpha}{z}\right)^k}{k}$$ and inserting this expression into previous formula for $K_n$ gives Taylor expansion that we derived for small $\alpha/z$ above.
You can use this expression to check your implementation of recurrence relation.