This section uses these add-on packages:
using CalculusWithJulia
using Plots; plotly()
using QuadGK
using RootsPrecompiling packages... 1301.9 ms ✓ Roots 528.7 ms ✓ Roots → RootsForwardDiffExt 899.8 ms ✓ CalculusWithJulia 3 dependencies successfully precompiled in 3 seconds. 64 already precompiled. Precompiling packages... 1657.4 ms ✓ CalculusWithJulia → CalculusWithJuliaPlotsExt 1 dependency successfully precompiled in 2 seconds. 213 already precompiled.
The fundamental theorem of calculus gives an easy to compute answer to the value of a definite integral when a computable (elementary) anti-derivative can be found. This is not always the case. See Liousville’s theorem to read more. If there is no computable anti-derivative the definite integral can be approximated numerically, as discussed in this section, where we begin with a Riemann sum approach, but end with a much more efficient Gauss-quadrature approach we will utilize in subsequent sections.
The Riemann sum approach gives a method to approximate the value of a definite integral. We just compute an approximating sum for a large value of \(n\), so large that the limiting value and the approximating sum are close.
To see the mechanics, let’s again return to Archimedes’ problem and approximate \(\int_0^1 x^2 dx\).
Let us fix some values, \(a\), \(b,\) and \(f\) are part of the question, \(n\) is related to the approximation.
a, b = 0, 1
f(x) = x^2
n = 55Then for a given \(n\) we have some steps to do: create the partition, find the \(c_i\), multiply the pieces, and add them up. Here is one way to do all this:
xs = a:(b-a)/n:b # also range(a, b, length=n)
deltas = diff(xs) # forms x2-x1, x3-x2, ..., xn-xn-1
cs = xs[1:end-1] # finds left-hand end points. xs[2:end] would be right-hand ones.0.0:0.2:0.8We want to sum the products \(f(c_i)\Delta_i\). Here is one way to do so using zip to iterate over the paired off values in cs and deltas.
sum(f(ci)*Δi for (ci, Δi) in zip(cs, deltas))0.24000000000000002Our answer is not so close to the value of \(1/3\), but what did we expect—we only used \(n=5\) intervals. Trying again with \(50,000\) gives us:
n = 50_000
xs = a:(b-a)/n:b
deltas = diff(xs)
cs = xs[1:end-1]
sum(f(ci)*Δi for (ci, Δi) in zip(cs, deltas))0.3333233333999998This value is about \(10^{-5}\) off from the actual answer of \(1/3\).
We should expect that larger values of \(n\) will produce better approximate values, as long as numeric issues don’t get involved.
Before continuing, we define a function to compute approximating sums for us with an extra argument to specifying one of four common methods for estimating \(\int_{x_{i-1}}^{x_i}f(x)dx\). Leaving explanations for later, Figure 39.1 shows the different approximations.
function riemann(f, xs; method="right")
Ms = (left = (f,a,b) -> f(a),
right = (f,a,b) -> f(b),
trapezoid = (f,a,b) -> (f(a) + f(b))/2,
simpsons = (f,a,b) -> (c = a/2 + b/2; (1/6) * (f(a) + 4*f(c) + f(b)))
)
M = Ms[Symbol(method)}
xs′ = zip(xs[1:end-1], xs[2:end])
sum(M(f, a, b) * (b-a) for (a,b) ∈ xs′)
end
riemann(f, a, b, n; method="right") =
riemann(f, range(a,b,n+1); method)(This function is defined in CalculusWithJulia and need not be copied over if that package is loaded.)
With this, we can easily find an approximate answer for a definite integral. We wrote the function to use the familiar template action(function, arguments...), so we pass in a function and arguments to describe the problem (a, b, and n and, optionally, the method):
f(x) = exp(x)
riemann(f, 0, 5, 10)187.324835773627Or with more intervals in the partition
riemann(f, 0, 5, 50_000)147.42052988337647(The answer is \(e^5 - e^0 = 147.4131591025766\dots\), which shows that even \(50,000\) partitions is not enough to guarantee many digits of accuracy.)
Numerically estimate the definite integral \(\int_0^2 x\log(x) dx\).
This particular integrand is continuous on \((0,2]\) but we can redefine it to be \(0\) at \(0\) to make it continuous on \([0,2]\), hence the integral above is well defined. Numerically though, we have to be a bit careful with the Riemann sum, as the left Riemann sum will have an issue at \(0=x_0\)—0*log(0) returns NaN which will poison any subsequent arithmetic operations, so the value returned will be NaN and not an approximate answer. We could define our function with a check, instead we avoid this value by using the right Riemann sum:
h(x) = x * log(x)
riemann(h, 0, 2, 50_000; method="right")0.38632208884775826(The default is "right", so no method specified would also work.)
The Riemann sum above is actually extremely inefficient in that it can take a large number of rectangles to produce an “accurate” approximation for a definite integral, even for nice functions. To see how much so, we can derive an estimate for the error in approximating the value using an arithmetic progression as the partition. Let’s assume that our function \(f(x)\) is increasing, so that the right sum gives an upper estimate and the left sum a lower estimate, so the error in any Riemann sum estimate will be smaller than the distance between these two values:
\[ \begin{align*} \text{error} &\leq \text{upper sum} - \text{lower sum}\\ &= \left(f(x_1) \cdot (x_{1} - x_0) + f(x_2) \cdot (x_{2} - x_1) + \cdots \right.\\ &\quad + \left. f(x_{n-1})(x_{n-1} - x_{n-2}) + f(x_n) \cdot (x_n - x_{n-1})\right)\\ &\quad - \left(f(x_0) \cdot (x_{1} - x_0) + f(x_1) \cdot (x_{2} - x_1) + \cdots \right.\\ &\quad + \left. f(x_{n-1})(x_n - x_{n-1}) \right)\\ &= \left(f(x_1) \cdot \Delta + f(x_2) \cdot \Delta + \cdots + f(x_{n-1})\Delta + f(x_n) \cdot \Delta y\right)\\ &\quad - \left(f(x_0) \cdot \Delta + f(x_1) \cdot \Delta + \cdots + f(x_{n-1})\Delta\right) \\ &= \left(\left[f(x_1) + f(x_2) + \cdots + f(x_n)\right] - \left[f(x_0) + \cdots + f(x_{n-1})\right]\right) \cdot \Delta \\ &= \left(f(b) - f(a)\right) \cdot \frac{b-a}{n}. \end{align*} \]
We see the error goes to \(0\) at a rate of \(1/n\) with the constant depending on \(b-a\) and the function \(f\). In general, a similar bound holds when \(f\) is not monotonic.
There are other ways to approximate the integral that use fewer points in the partition. Riemann sums approximate the definite integral over each “piece” of the partition—\(\int_{x_{i-1}}^{x_i} f(x) dx\)—using a rectangle. Other geometric shapes are possible.
The trapezoid rule uses a trapezoid formed to approximate this area, namely the one formed by \((x_{i-1}, 0)\), \((x_{i-1}, f(x_{i-1}))\), \((x_i, f(x_i))\), and \((x_i, 0)\) with area
\[ \frac{1}{2} \left(f(x_{i-1}) + f(x_i) \right) \cdot (x_i - x_{i-1}). \]
If we use an equally spaced partition (\(\Delta=(b-a)/n\)) and add all the \(n\) terms, we get single contributions from the endpoints and double from the others giving
\[ A \approx (\frac{f(x_0)}{2} + \frac{f(x_n)}{2})\Delta + \sum_{i=1}^{n-1} f(x_i) \Delta/ \]
In a later section, we will see that the error in using trapezoids to estimate the area is bounded, for some constant \(K\):
\[ \text{error} \leq \frac{K (b-a)^3}{12n^2}. \]
The \(n^2\) means roughly that the error in the estimate using a Riemann sum with \(n\) terms is similar to the error in the estimate using the trapezoid rule with \(\sqrt{n}\) terms.1
Consider the integral
\[ \int_0^2 x e^{-x} dx = 1 - 3 e^{-2} = 0.59399415\cdots \]
For comparison sake, we define the exact answer as a constant:
A = 1 - 3 * exp(-2)0.5939941502901619The error of a Riemann sum with \(n=10^4\) is then:
a, b = 0, 2
f(x) = x * exp(-x)
riemann(f, a, b, 10^4; method="right") - A2.7063272193594834e-5and this is comparable to the error of the trapezoid method with \(n=10^2\):
riemann(f, a, b, 10^2; method="trapezoid") - A-3.7843872858767114e-5Simpson’s rule is one, where instead of approximating the area with rectangles that go through some \(c_i\) in \([x_{i-1}, x_i]\) instead the function is approximated by the quadratic polynomial going through \(x_{i-1}\), \((x_i + x_{i-1})/2\), and \(x_i\) and the exact area under that polynomial is used in the approximation. The explicit formula for a single partition is2
\[ \int_{x_{i-1}}^{x_i} f(x) dx \approx \frac{x_i - x_{i-1}}{6}\left(f(x_{i-1}) + 4f(\frac{x_{i-1} + x_i}{2}) + f(x_i)\right) \]
The error in this approximation can be shown to be
\[ \text{error} \leq \frac{(b-a)^5}{180n^4} \text{max}_{\xi \text{ in } [a,b]} \lvert f^{(4)}(\xi) \rvert. \]
That is, the error is like \(1/n^4\) with constants depending on the length of the interval, \((b-a)^5\), and the maximum value of the fourth derivative over \([a,b]\). This is significant, the error in \(10\) steps of Simpson’s rule is on the scale of the error of \(10,000\) steps of the Riemann sum for well-behaved functions.
The Wikipedia article mentions that Kepler used a similar formula \(100\) years prior to Simpson, or about \(200\) years before Riemann published his work. Again, the value in Riemann’s work is not the computation of the answer, but the framework it provides in determining if a function is Riemann integrable or not.
Continuing the previous example, the accuracy of Simpson’s rule with \(10\) steps is comparable to that of a Riemann sum with \(10^4\) steps:
riemann(f, a, b, 10; method="simpsons") - A-1.5884464277249322e-6There are function types where the above approximations are actually exact:
This pattern could be extended by taking more intermediate points. In fact an entire family of similar approximations using \(n\) points can be made exact for any polynomial of degree \(n-1\) or lower. However, by choosing points judiciously—not necessarily evenly spaced out and not necessarily including the end points—\(n\) points can be exact for polynomials of degree higher than \(n\). (Simpson’ rule actually being exact for cubic polynomials is something that hints at this.)
The formulas for an approximation to the integral \(\int_{-1}^1 f(x) dx\) discussed so far can be written as:
\[ \begin{align*} S &= f(x_1) \Delta_1 + f(x_2) \Delta_2 + \cdots + f(x_n) \Delta_n\\ &= w_1 f(x_1) + w_2 f(x_2) + \cdots + w_n f(x_n)\\ &= \sum_{i=1}^n w_i f(x_i). \end{align*} \]
The \(w\)s are “weights” and the \(x\)s are nodes. Restricting to the interval \([-1,1]\) presents no loss in generality.
A Gaussian quadrature rule is a set of weights and nodes for \(i=1, \dots n\) for which the sum is exact for any \(f\) which is a polynomial of degree \(2n-1\) or less. Such choices then also approximate well the integrals of functions which are not polynomials of degree \(2n-1\) or less, provided \(f\) can be well approximated by a polynomial over \([-1,1]\). (Which is the case for the “nice” functions we encounter, though not for highly oscillatory functions.) More details are discussed in the section on orthogonal polynomials and some examples are given in the questions.
In Julia a modification of the Gauss quadrature rule is implemented in the quadgk function (from the QuadGK package) to give numeric approximations to integrals. The quadgk function also has the familiar interface action(function, arguments...). Unlike our riemann function, there is no n specified, as the number of steps is adaptively determined. (There is more partitioning occurring where the function is changing rapidly.) Instead, the algorithm outputs an estimate on the possible error along with the answer. Instead of \(n\), some trickier problems require a specification of an error threshold.
To use the function to integrate f over an interval [a,b] we have:
f(x) = x * log(x)
quadgk(f, 0, 2)(0.38629436103070175, 4.856575112145755e-9)As mentioned, there are two values returned: an approximate answer, and an error estimate. In this example we see that the value of \(0.3862943610307017\) is accurate to within \(10^{-9}\). (The actual answer is \(-1 + 2\cdot \log(2)\) and the error is only \(10^{-11}\). The reported error is an estimated upper bound, and may be conservative, as with this problem.) Our previous answer using \(50,000\) right-Riemann sums was \(0.38632208884775737\) and is only accurate to \(10^{-5}\). By contrast, this method uses just \(256\) function evaluations in the above problem.
The method should be exact for polynomial functions:
f(x) = x^5 - x + 1
quadgk(f, -2, 2)(3.9999999999999973, 1.3322676295501878e-15)The error term is \(0\), the answer is \(4\) up to the last unit of precision (1 ulp), so any error is only in floating point approximations.
For the numeric approximation of a definite integral, the quadgk function should be preferred over the other methods previously discussed.
Here are some sample integrals computed with quadgk:
\[ \int_0^\pi \sin(x) dx \]
quadgk(sin, 0, pi)(2.0000000000000004, 1.7901236049056024e-12)(Again, the actual answer is off only in the last digit, the error estimate is an upper bound.)
\[ \int_0^5 e^x dx \]
quadgk(exp, 0, 5)(147.41315910257657, 2.6594506152832764e-8)\[ \int_0^2 x^x dx \]
u(x) = x^x
quadgk(u, 0, 2)(2.8338767448900546, 1.948175154531384e-8)The function \(x^x\) is not continuous at \(0\), but can be defined to be so. In this case, the numeric definition of 0^0 matches the limit, so no discussion of redefining the function is necessary, as was done earlier with the function \(x\cdot \log(x)\).
In fact, the specified endpoints to quadgk are never evaluated, so such concerns are not needed. (Which can be exploited when integrals involving functions with vertical asymptotes are discussed.) This is why the first example—which integrated x*log(x)—did not return NaN but rather an estimate for the integral.
When composing the answer with other functions it may be desirable to drop the error in the answer, we discuss three styles that can be used for this. The first is to just name the two returned values:
A, err = quadgk(cos, 0, pi/4)
A0.7071067811865475The second is to ask for just the first component of the returned value:
A = first(quadgk(tan, 0, pi/4))0.3465735902799726Finally, direct indexing can be applied, as with
quadgk(tan, 0, pi/4)[1]0.3465735902799726Though we try to avoid this style in favor of being more explicit when that is convenient.
In probability theory, a univariate density is a function, \(f(x)\) such that \(f(x) \geq 0\) and \(\int_a^b f(x) dx = 1\), where \(a\) and \(b\) are the range of the distribution.
The Von Mises distribution, takes the form
\[ k(x) = C \cdot \exp(\cos(x)), \quad -\pi \leq x \leq \pi. \]
Compute \(C\) (numerically).
The fact that \(1 = \int_{-\pi}^\pi C \cdot \exp(\cos(x)) dx = C \int_{-\pi}^\pi \exp(\cos(x)) dx\) implies that \(C\) is the reciprocal of the definite integral:
k(x) = exp(cos(x))
A, err = quadgk(k, -pi, pi)(7.954926521012919, 3.9023298370466364e-8)So
C = 1/A
k₁(x) = C * exp(cos(x))k₁ (generic function with 1 method)The cumulative distribution function for \(k(x)\) is \(K(x) = \int_{-\pi}^x k(u) du\), \(-\pi \leq x \leq \pi\). We just showed that \(K(\pi) = 1\) and it is trivial that \(K(-\pi) = 0\). The quantiles of the distribution are the values \(q_1\), \(q_2\), and \(q_3\) for which \(K(q_i) = i/4\). Can we find these?
First we define a function, that computes \(K(x)\). We only need the first of the two answers given by quadgk.
K(x) = first(quadgk(k₁, -pi, x))K (generic function with 1 method)The question asks us to solve \(K(x) = 0.25\), \(K(x) = 0.5\) and \(K(x) = 0.75\). The Roots package can be used for such work, in particular find_zero. We will use a bracketing method, as clearly \(K(x)\) is increasing, as \(k(u)\) is positive, so we can just bracket our answer with \(-\pi\) and \(\pi\). (We solve \(K(x) - p = 0\), so \(K(\pi) - p > 0\) and \(K(-\pi)-p < 0\).). We could do this with a comprehension, but for variety use broadcasting with solve below.
Z = ZeroProblem((x,p) -> K(x) - p, (-pi, pi))
solve.(Z, (1/4, 1/2, 3/4))(-0.8097673745015915, 0.0, 0.8097673745016285)The middle one is clearly \(0\). This distribution is symmetric about \(0\), so half the area is to the right of \(0\) and half to the left, so clearly when \(p=0.5\), \(x\) is \(0\). The other two show that the area to the left of \(-0.809767\) is equal to the area to the right of \(0.809767\) and equal to \(0.25\).
quadgkTo visualize the choice of nodes by the algorithm, In Figure 39.2 the nodes chosen are shown for \(f(x)=\sin(x)\) over \([0,\pi]\). Relatively few nodes used to get a high-precision estimate.
quadgk for \(f(x) = \sin(x)\) over \([0, \pi]\)
For a more oscillatory function, more nodes are chosen, as seen in Figure 39.3.
quadgk for \(f(x) = e^{x} \sin(\pi x)\) over \([0, \pi]\). There are more nodes chosen then for the function \(\sin(x)\), as seen in Figure 39.2
In both Figure 39.2 and Figure 39.3 it can be verified that no node is one of the endpoints.
The QuadGK.gauss(n) function returns a pair of \(n\) quadrature points and weights to integrate a function over the interval \((-1,1)\), with an option to use a different interval \((a,b)\). For a given \(n\), these values exactly integrate any polynomial of degree \(2n-1\) or less. In this example, these \(5\) points produce an answer accurate already to the \(5\)th decimal point.
xs, ws = QuadGK.gauss(5)([-0.906179845938664, -0.5384693101056831, 0.0, 0.5384693101056831, 0.906179845938664], [0.2369268850561892, 0.4786286704993663, 0.5688888888888887, 0.4786286704993663, 0.2369268850561892])f(x) = exp(cos(x))
sum(w * f(x) for (x, w) in zip(xs, ws))4.68315881456892The pattern to integrate can be expressed in other ways, but using the zip function to iterate over the xs and ws as pairs of values is pretty direct.
For the function \(f(x) = \sin(\pi x)\), estimate the integral for \(-1\) to \(1\) using a left-Riemann sum with the partition \(-1 < -1/2 < 0 < 1/2 < 1\).
For the right Riemann sum approximating \(\int_0^{10} e^x dx\) with \(n=100\) subintervals, what would be a good estimate for the error?
Use quadgk to find the following definite integral:
\[ \int_1^4 x^x dx . \]
Use quadgk to find the following definite integral:
\[ \int_0^3 e^{-x^2} dx . \]
Use quadgk to find the following definite integral:
\[ \int_0^{9/10} \tan(u \frac{\pi}{2}) du. \]
Use quadgk to find the following definite integral:
\[ \int_{-1/2}^{1/2} \frac{1}{\sqrt{1 - x^2}} dx \]
Let \(A=1.98\) and \(B=1.135\) and
\[ f(x) = \frac{1 - e^{-Ax}}{B\sqrt{\pi}x} e^{-x^2}. \]
Find \(\int_0^1 f(x) dx\)
A bound for the complementary error function ( positive function) is
\[ \text{erfc}(x) \leq \frac{1}{2}e^{-2x^2} + \frac{1}{2}e^{-x^2} \leq e^{-x^2} \quad x \geq 0. \]
Let \(f(x)\) be the first bound, \(g(x)\) the second. Assuming this is true, confirm numerically using quadgk that
\[ \int_0^3 f(x) dx \leq \int_0^3 g(x) dx \]
The value of \(\int_0^3 f(x) dx\) is
The value of \(\int_0^3 g(x) dx\) is
The function in the interactive graph of Figure 39.4 is
\[ f(x) = \frac{1}{\sqrt{ x^4 + 10x^2 - 60x + 100}}. \]
When \(n=5\) what is the area of the Riemann sum?
When \(n=50\) what is the area of the Riemann sum?
Using quadgk what is the area under the curve?
Gauss nodes for approximating the integral \(\int_{-1}^1 f(x) dx\) for \(n=4\) are:
ns = [-0.861136, -0.339981, 0.339981, 0.861136]4-element Vector{Float64}:
-0.861136
-0.339981
0.339981
0.861136The corresponding weights are
wts = [0.347855, 0.652145, 0.652145, 0.347855]4-element Vector{Float64}:
0.347855
0.652145
0.652145
0.347855Use these to estimate the integral \(\int_{-1}^1 \cos(\pi/2 \cdot x)dx\) with \(w_1f(x_1) + w_2 f(x_2) + w_3 f(x_3) + w_4 f(x_4)\).
The actual answer is \(4/\pi\). How far off is the approximation based on 4 points?
Using the Gauss nodes and weights from the previous question, estimate the integral of \(f(x) = e^x\) over \([-1, 1]\). The value is:
There are functions where the trapezoid method has much faster convergence, even exponential. (cf. this article and these notes for some background.)↩︎
The riemann function sums this expression, but this approach is inefficient computationally. Alternative formulations would be suggested.↩︎