Crash Course in Scientific Computing
XIII. Integrals
We will now start to apply algorithms, or precise recipes, to numerical problems. We will start with a simple one, that has a clear physical interpretation: integration of 1D functions of one variable, i.e., the computation of areas. We will also look at a precise application of such computations (which are of interest on their own, of course). Namely, we will look at how to generate random numbers on our own, that is, following any given probability distribution (a positive, normalized function). For instance, let us say we want to simulate the ground state of a particle in a box, with probability distribution over $[0,1]$ given by:
$$\psi^2(x)=2\sin^2(\pi x)$$
We remind that one can use unicode characters (by entering \psi+TAB) to call our density of probability ψ2, which in quantum mechanics is given by the modulus square of the probability amplitude (in 1D the wavefunction can always be taken real so we do not need worry about the modulus):
ψ2(x)=2*sin(pi*x)^2
That's our density of probability:
using LaTeXStrings default(; w=3) plot(ψ2,0,1, title=L"2\sin^2(x)", xlabel=L"x", ylabel=L"\psi^2(x)",dpi=150,legend=false)
where we decorated the plot with a lot of options, but the width which we gave as a default attribute, not to do it in the future (and since we like thick lines).
To draw random numbers that follow this probability distribution, we will need to compute quantities such as $F(y)\equiv\int_0^y\psi^2(x)\,dx$, or the so-called cumulative. Of course we should find that $F(1)=1$ (normalization condition), so we can check that as well. The simplest way to compute such integrals numerically is to bring to the computer the mathematical definition we have given of the integral, namely, in terms of Riemann sums, i.e., as a sum of parallelograms defined on a partition of the support of the function (here $[0,1]$) cut into $n$ intervals. Calling $f$ the function to integrate, this is defined as:
$$R_n=\sum_{i=1}^{n}f(x_i)\delta(i)$$
a, b = 0, 1; # boundaries lx=range(a, stop=b, length=11) # list of x intervals δ=diff(lx) # widths of intervals lxi=lx[1:end-1] # x_i points of parallelograms
We can now compute our integral over the entire interval to check the normalization:
sum([ψ2(lxi[i])*δ[i] for i=1:length(lxi)])
1.0
The cumulative is itself obtained as:
Ψ=[sum([ψ2(lxi[i])*δ[i] for i=1:j]) for j=1:length(lxi)]
This is the result for 100 intervals (length=101 above):
plot(lxi, Ψ, legend=false, title="Cumulative")
Now the generation of a random number following the initial distribution works as follows. We select randomly (uniformly) a number on the $y$-axis of the cumulative and find the corresponding $x$ such that $F(x)=y$. Those $x$ are $f$ distributed.
This can be simply achieved as follows (later we'll see a more sophisticated way but let's just use the most straightforward way to get to the point):
yr=rand() fx=lx[[findmin(abs.(Ψ.-rand()))[2] for i=1:10]]
Here are some positions where our ground state collapsed:
10-element Array{Float64,1}:
0.55
0.96
0.24
0.5
0.58
0.16
0.67
0.74
0.08
0.46
histogram(lx[[findmin(abs.(Ψ.-rand()))[2] for i=1:10^7]], bins=0:.01:1,norm=true)
Note that the binning cannot be less than the $\delta$ step, otherwise we will get "holes" (and also the wrong normalization) (?!)Too fine binning, resulting in holes and wrong normalization.
.
Putting our results together in a function of our own:
function myRiemann(f, n) a, b = 0, 1; # boundaries lx=range(a, stop=b, length=n) # list of x intervals δ=diff(lx); # widths of intervals lxi=lx[1:end-1]; # x_i points of parallelograms sum([f(lxi[i])*δ[i] for i=1:length(lxi)]) end
Then we can try:
myRiemann(g, 10^3)
1.717421971020861
The result is fairly close to the exact
$$e-1\approx1.718281828459045\tag{1}$$
However we'll find that our method is fairly limited:
scatter([abs(myRiemann(g, i)-(exp(1)-1)) for i=2:100:1002], ylims=(10^-4,1), yscale=:log10, xticks=(1:11, 2:100:1002), xlabel="Num steps", ylabel="error")
It is the purpose of numerical methods to learn how to use algorithms that are efficient, in the sense that they are accurate, fast and resource-effective.
It is easy to be brutal and terribly inefficient with a computer. In fact a fairly trivial enhancement of our method leads to considerable improvement. We have used Riemann integrals by taking the leftmost point. We could also take the rightmost point, or the central point, i.e., their average:
function myRiemann(f, n, method="mid") a, b = 0, 1; # boundaries lx=range(a, stop=b, length=n) # list of x intervals δ=diff(lx); # widths of intervals if method == "left" lxi=lx[1:end-1]; # x_i points of parallelograms elseif method == "right" lxi=lx[2:end]; else lxi=(lx[1:end-1]+lx[2:end])/2; end sum([f(lxi[i])*δ[i] for i=1:length(lxi)]) end
There is a huge gain in precision by evaluating the function elsewhere? Why? Clearly because we find a compensation of errors. Note that evaluating at the central point is equivalent to computing the area of a trapezoid, which also explain why the area is better estimated. For this reason, this method of integration is called the "trapezoidal rule".
Can it be further improved? Of course, following the same idea, one could approximate the function not by a straight line, but by a curve. Polynomials make easy curves and they are easy to integrate too, not even numerically, but symbolically. That brings forward the so-called "Newton-Cotes quadrature" that specify the weights $\omega_i$ and the locations $x_i$ where to sample the function to integrate:
$$\int _{a}^{b}f(x)\,dx\approx \sum _{i=0}^{n}w_{i}\,f(x_{i})$$
A great choice of interpolating polynomials is the so-called Lagrange-interpolating polynomials, which provide the unique polynomial of order $n$ to pass by the $n$ points $(x_i, y_i)$ (with $1\le i\le n$). Note that the order of the polynomial is the same as the number of points one wishes to pass through. These polynomials are obtained as a linear superposition of the Lagrange-basis polynomials:
$$\ell _{j}(x)\equiv\prod _{\begin{smallmatrix}1\leq k\leq n\\k\neq j\end{smallmatrix}}{\frac {x-x_{k}}{x_{j}-x_{k}}}={\frac {(x-x_{1})}{(x_{j}-x_{1})}}\cdots {\frac {(x-x_{j-1})}{(x_{j}-x_{j-1})}}{\frac {(x-x_{j+1})}{(x_{j}-x_{j+1})}}\cdots {\frac {(x-x_{n})}{(x_{j}-x_{n})}}$$
They are such that $l_j(x_k)=0$ for all $1\le k\le n$ such that $k\neq j$ where $l_j(x_k)=1$. The interpolating polynomial $P_f$ of the function $f$ sampled at the points $x_k$ is then obtained as a linear superposition of the Lagrange polynomials with the values of the function as the weighting coefficients:
$$P_f(x)=\sum_{j=1}^n f(x_j)\ell_j(x)$$
Given some points $(x_i, y_i)$, e.g.,
xi=[1, 2, 3, 5]; yi=[1, 3, 1, 2];
plot the Lagrange polynomials allowing to interpolate through all these points along with the interpolation polynomial itself. Your code should display something like this on the previous input:
function plotLagrange(xi, yi) scatter((xi, yi), label="Data points", legend=:topleft); for j=1:length(xi) plot!(x->prod([(x-xi[k])/(xi[j]-xi[k]) for k∈union(1:j-1,j+1:length(xi))]), minimum(xi)-1, maximum(xi)+1, label="l$(j)(x)", ls=:dash) end plot!(x->sum([yi[j]*prod([(x-xi[k])/(xi[j]-xi[k]) for k∈union(1:j-1,j+1:length(xi))]) for j∈1:length(xi)]), minimum(xi)-1, maximum(xi)+1, lw=3, lc=:red, label="Lagrange interpolation", xlabel="x", ylabel="f(x)") end
For instance, let us interpolate with two points (a line), at $x_1$ and $x_2=x_1+h$. The Lagrange basis in this case reads:
\begin{align} l_1(x)&={x-x_2\over x_1-x_2}\\ l_2(x)&={x-x_2\over x_2-x_1} \end{align}
and the interpolating Lagrange polynomial becomes:
$$P_f(x)=f_1l_1(x)+f_2l_2(x)$$
where $f_i\equiv f(x_i)$, which simplifies to (collecting $x$ and substituting $h=x_2-x_1$):
$$P_f(x)=x{f_2-f_1\over h}+(f_1+{x_1\over h}f_1-{x_1\over h}f_2)$$
So the integral can be found by performing it symbolically on the polynomials (which are easy to integrate exactly) and turn the problem into mere additions of weighted values of the functions at given points (quadrature):
\begin{align} \int_{x_1}^{x_2}f(x)\,dx&\approx\int_{x_1}^{x_2}P_f(x)\,dx\\ &={f_2-f_1\over h}\left.{x^2\over 2}\right|_{x_1}^{x_1+h}+(f_1+{x_1\over h}f_1-{x_1\over h}f_2)x\left.\vphantom{x^2\over2}\right|_{x_1}^{x_1+h}\\ &={f_1+f_2\over 2h} \end{align}
and we recover the trapezoidal rule.
With three points, at $x_1$, $x_2=x_1+h$ and $x_3=x_1+2h$, the Lagrange polynomials read:
\begin{align} l_1(x)&={x-x_2\over x_1-x_2}{x-x_3\over x_1-x_3}\\ l_2(x)&={x-x_1\over x_2-x_1}{x-x_3\over x_2-x_3}\\ l_3(x)&={x-x_1\over x_3-x_1}{x-x_2\over x_3-x_2} \end{align}
and the Lagrange interpolating polynomial:
$$P_f(x)=f_1l_1(x)+f_2l_2(x)+f_3l_3(x)$$
simplifies to:
$$ {f_1\over 2h^2}(x^2-x(x_2+x_3)+x_2x_3) -{f_2\over h^2}(x^2-x(x_1+x_3)+x_1x_3) +{f_3\over 2h^2}(x^2-x(x_1+x_2)+x_1x_2) $$
Now one can integrate each polynomial in turn. We work out the one that goes with the term $f_1$ (others follow similarly):
\begin{align} \int_{x_1}^{x_1+2h}&(x^2-x(x_2+x_3)+x_2x_3)\,dx={}\\ &{(x_1+2h)^3-x_1^3\over 3}-{(x_1+2h)^2-x_1^2\over 2}(x_2+x_3)+((x_1+2h)-x_1)x_2x_3\\ &={6x_1^2h+12x_1h^2+8h^3\over 3}-{4x_1h+4h^2\over 2}(2x_1+3h)+2h(x_1^2+3x_1h+2h^2)\\ &={2\over 3}h^3 \end{align}
So, incredibly, there is massive simplifications and one is left with a simple ${f_1\over2h^2}\times{2\over 3}h^3=f_1{h\over 3}$ expression meaning a $h/3$ weighting of the first sampled point. We let you compute the other coefficients yourself:
Show that simplifications of the integrated Lagrange polynomials result in the simple quadrature formula for a cubic (or three-points) integration of $f$:
$$\int_{x_1}^{x_3}P_f(x)\,dx={h\over 3}(f_1+4f_2+f_3)\,.\tag{2}$$
The following piece of code implements the 3-points ($n=3$) Newton Cotes method:
# Integrate f between a and b with a n-point method: function NewtonCotes(f, a, b, n=3) h=(b-a)/n if n==3 (h/3)*(f(a)+4f((a+b)/2)+f(b)) else println("Sorry, $(n)-point method is not available.") end end julia> NewtonCotes(exp,0,1) 1.7188611518765928
so quite accurate by using 3 points only (the exact result is Eq. (1))! To increase accuracy one can do two things, increase $n$, or break the functions into smaller intervals and iterate the method. The latter is natural in spirit as extending the Riemann idea to Lagrange-interpolated functions. This is called a "composite Newton-Cotes method" and can be implemented either through computer-assisted iteration:
julia> δx=.1 0.1 julia> sum([NewtonCotes(exp,i,i+δx) for i=0:δx:1-δx]) 1.7182818881038564
or one can provide directly the expression that sums the Newton Cotes formula on each subinterval, e.g.,
$${\displaystyle {\begin{aligned}\int _{a}^{b}f(x)\,dx&\approx {\frac {h}{3}}\sum _{j=1}^{n/2}{\bigg [}f(x_{2j-2})+4f(x_{2j-1})+f(x_{2j}){\bigg ]}\\{}&={\frac {h}{3}}{\bigg [}f(x_{0})+2\sum _{j=1}^{n/2-1}f(x_{2j})+4\sum _{j=1}^{n/2}f(x_{2j-1})+f(x_{n}){\bigg ]},\tag{3}\end{aligned}}}$$
where $x_j=a+jh$ for $0\le j\le n$ with $h=(b-a)/n$.
Prove Eq. (3).
This is the error made with this 3-point method:
lϵ=abs.([sum([NewtonCotes(exp,i,i+δx) for i=0:δx:1-δx])-(exp(1)-1) for δx in lδx]) 6-element Array{Float64,1}: 5.964481131215393e-8 5.9676708019651414e-12 2.220446049250313e-15 1.7186252421197423e-13 4.369393735714766e-12 2.1261659099991448e-11 plot(lϵ, yscale=:log10, label="3-point Newton Cotes") scatter!(lϵ, yscale=:log10, label=false, xlabel=L"h=10^{-k}", ylabel="error", dpi=300)
The idea can be extended to higher-orders of interpolations, that is, with more points in the interval. Possibly the most famous formula is the Simpson 3/8 rule which is a healthy practice in awkward (but easy) algebra to bring to its conclusion to find that:
Show that the 4-point interpolation provides the (so-called Simpson's 3/8) rule:
$$\int_{x_1}^{x_4} P_f(x)\,dx={\displaystyle {\frac {3h}{8}}(f_{0}+3f_{1}+3f_{2}+f_{3})}$$
where $f_i\equiv f(x_1+ih)$ with step $h\equiv(x_4-x_1)/3$ for $1\le i\le 3$.
This is better than the 3-point method but still of order $h^5$. To go to order $h^7$ one needs Boole's rule which reads:
$$\int_{x_1}^{x^5}P_f(x)\,dx={\displaystyle {\frac {2h}{45}}(7f_{0}+32f_{1}+12f_{2}+32f_{3}+7f_{4})}$$
with the same definitions as above (now $h\equiv(x_5-x_1)/4$) and that you could similarly demonstrate, if interested. You can also find the cases elsewhere (e.g., here or here). Here is an implementation of cases up to 7-point integration:
function NewtonCotes(f, a, b, n=4) h=(b-a)/(n-1) if n==2 (h/2)*(f(a)+f(b)) elseif n==3 (h/3)*(f(a)+4f((a+b)/2)+f(b)) elseif n==4 (3h/8)*(f(a)+3f(a+h)+3f(a+2h)+f(b)) elseif n==5 (2h/45)*(7f(a)+32f(a+h)+12f(a+2h)+32f(a+3h)+7f(a+4h)) elseif n==6 (5h/288)*(19f(a)+75f(a+h)+50f(a+2h)+50f(a+3h)+75f(a+4h)+19f(a+5h)) elseif n==7 (h/140)*(41f(a)+216f(a+h)+27f(a+2h)+272f(a+3h)+27f(a+4h)+216f(a+5h)+41f(a+6h)) else println("Sorry, $(n)-point method is not available.") end end
which it is now easy to benchmark by compositing the function over steps of smaller sizes, e.g.:
julia> lδx=[2.0^(-k) for k=1:16] 16-element Array{Float64,1}: 0.5 0.25 0.125 0.0625 0.03125 0.015625 0.0078125 0.00390625 0.001953125 0.0009765625 0.00048828125 0.000244140625 0.0001220703125 6.103515625e-5 3.0517578125e-5 1.52587890625e-5
which gives us for $\int_0^1 e^x\,dx$:
p=plot(); for n=2:7 lϵ=filter(!iszero, abs.([sum([NewtonCotes(exp,i,i+δx,n) for i=0:δx:1-δx])-(exp(1)-1) for δx in lδx])); plot!(lϵ, yscale=:log10, label="$(n)-point") scatter!(lϵ, label=false, c=:red, m=2) end plot!(xticks = (1:length(lδx), ["\$2^{-$(k)}\$" for k∈1:length(lδx)]), dpi=200, xformatter = :latex, xlabel = "Step in composite Newton Cotes", ylabel="Error")
Error in the various $n$-point Newton Cotes composite integrations on $\int_0^1e^x\,dx$. The machine precision is $\approx 10^{-16}$. Points not shown mean that the error is zero.
The higher-order method quickly reach machine precision and thus cannot be more accurate: this is a limitation of the computer, not of the method itself. When points are not shown, this is because the result is exactly the same as found by exp(1)-1 from the computer, and one cannot hope for a better result than this one (with machine precision). So we have basically achieved (numerical) perfection!
Of course, other people already wrote such functions, that are available through packages. For numerical integration of 1D functions, one can use QuadGK which is based on the so-called Gauss–Kronrod quadrature formula, which basically figures out how to best choose the points $x_i$ and how to weight them:
using QuadGK @time quadgk(g, 0, 1)
Unlike our naive implementation, the result is pretty accurate:
0.000034 seconds (9 allocations: 304 bytes) (1.718281828459045, 0.0)
That's the best we can achieve (with mid-points and a hundred-million intervals!)
@time myRiemann(g, 10^8)
9.762264 seconds (500.00 M allocations: 8.941 GiB, 11.01% gc time)
1.7182818284590453
It differs from the exact result only in the last digit, but took about 10s & 9GiB of memory.