There is a more compact notation to $x_1 + x_2 + \cdots + x_n$, this using the *summation notation* or capital sigma. We have:
```math
\Sigma_{i = 1}^n x_i = x_1 + x_2 + \cdots + x_n
```
The notation includes three pieces of information:
- The $\Sigma$ is an indication of a sum
- The ${i=1}$ and $n$ sub- and superscripts indicate the range to sum over.
- The term $x_i$ is a general term describing the $i$th entry, where it is understood that $i$ is just some arbitrary indexing value.
With this notation, a Riemann sum can be written as $\Sigma_{i=1}^n f(c_i)(x_i-x_{i-1})$.
### Other sums
The choice of the $c_i$ will give different answers for the
approximation, though for an integrable function these differences
will vanish in the limit. Some common choices are:
* Using the right hand endpoint of the interval $[x_{i-1}, x_i]$ giving the right-Riemann sum, ``R_n``.
* The choice $c_i = x_{i-1}$ gives the left-Riemann sum, ``L_n``.
* The choice $c_i = (x_i + x_{i-1})/2$ is the midpoint rule, ``M_n``.
* If the function is continuous on the closed subinterval $[x_{i-1}, x_i]$, then it will take on its minimum and maximum values. By the extreme value theorem, we could take $c_i$ to correspond to either the maximum or the minimum. These choices give the "upper Riemann-sums" and "lower Riemann-sums".
The choice of partition can also give different answers. A common
choice is to break the interval into $n+1$ equal-sized pieces. With
$\Delta = (b-a)/n$, these pieces become the arithmetic sequence $a =
a + 0 \cdot \Delta < a + 1 \cdot \Delta < a + 2 \cdot \Delta < \cdots
a + n \cdots < \Delta = b$ with $x_i = a + i (b-a)/n$. (The
`range(a, b, length=n+1)` command will compute these.) An alternate choice
made below for one problem is to use a geometric progression:
```math
a = a(1+\alpha)^0 < a(1+\alpha)^1 < a (1+\alpha)^2 < \cdots < a (1+\alpha)^n = b.
```
The general statement allows for any partition such that the largest gap
goes to ``0``.
----
Riemann sums weren't named after Riemann because he was the first to
approximate areas using rectangles. Indeed, others had been using even
more efficient ways to compute areas for centuries prior to
Riemann's work. Rather, Riemann put the definition of the area under
the curve on a firm theoretical footing with the following
theorem which gives a concrete notion of what functions
are integrable:
> **Riemann Integral**: A function $f$ is Riemann integrable over the
> interval $[a,b]$ and its integral will have value $V$ provided for every
> $\epsilon > 0$ there exists a $\delta > 0$ such that for any
> When the integral exists, it is written $V = \int_a^b f(x) dx$.
!!! info "History note"
The expression $V = \int_a^b f(x) dx$ is known as the *definite integral* of $f$ over $[a,b]$. Much earlier than Riemann, Cauchy had defined the definite integral in terms of a sum of rectangular products beginning with $S=(x_1 - x_0) f(x_0) + (x_2 - x_1) f(x_1) + \cdots + (x_n - x_{n-1}) f(x_{n-1})$ (the left Riemann sum). He showed the limit was well defined for any continuous function. Riemann's formulation relaxes the choice of partition and the choice of the $c_i$ so that integrability can be better understood.
### Some immediate consequences
The following formulas are consequences when $f(x)$ is
integrable. These mostly follow through a judicious rearranging of the
approximating sums.
The area is $0$ when there is no width to the interval to integrate over:
> $\int_a^a f(x) dx = 0.$
Even our definition of a partition doesn't really apply, as we assume
$a < b$, but clearly if $a=x_0=x_n=b$ then our only"approximating" sum
could be $f(a)(b-a) = 0$.
The area under a constant function is found from the area of
rectangle, a special case being $c=0$ yielding $0$ area:
The scaling operation shifts ``a`` to ``ca`` and ``b`` to ``cb`` so the limits of integration make sense. However, the area stretches by ``c`` in the ``x`` direction, so must contract by ``c`` in the ``y`` direction to stay in balance. Hence the factor of ``1/c``.
Combining two operations above, the operation ``g(x) = \frac{1}{h}f(\frac{x-c}{h})`` will leave the area between ``a`` and ``b`` under ``g`` the same as the area under ``g`` between ``(a-c)/h`` and ``(b-c)/h``.
----
The area between $a$ and $b$ can be broken up into the sum of the area
Cauchy showed this using a *geometric series* for the partition, not the arithmetic series $x_i = a + i (b-a)/n$. The series defined by $1 + \alpha = (b/a)^{1/n}$, then $x_i = a \cdot (1 + \alpha)^i$. Here the bases $x_{i+1} - x_i$ simplify to $x_i \cdot \alpha$ and $f(x_i) = (a\cdot(1+\alpha)^i)^k = a^k (1+\alpha)^{ik}$, or $f(x_i)(x_{i+1}-x_i) = a^{k+1}\alpha[(1+\alpha)^{k+1}]^i$,
so, using $u=(1+\alpha)^{k+1}=(b/a)^{(k+1)/n}$, $f(x_i) \cdot(x_{i+1} - x_i) = a^{k+1}\alpha u^i$. This gives
Again, Cauchy showed this using a geometric series. The expression $f(x_i) \cdot(x_{i+1} - x_i)$ becomes just $\alpha$. So the approximating sum becomes:
Certainly other integrals could be computed with various tricks, but
we won't pursue this. There is another way to evaluate integrals using
the forthcoming Fundamental Theorem of Calculus.
### Some other consequences
* The definition is defined in terms of any partition with its norm
bounded by $\delta$. If you know a function $f$ is Riemann
integrable, then it is enough to consider just a regular partition
$x_i = a + i \cdot (b-a)/n$ when forming the sums, as was done above. It is just that
showing a limit for just this particular type of partition would not be
sufficient to prove Riemann integrability.
* The choice of $c_i$ is arbitrary to allow for maximum
flexibility. The Darboux integrals use the maximum and minimum over
the subinterval. It is sufficient to prove integrability to show
that the limit exists with just these choices.
* Most importantly,
> A continuous function on $[a,b]$ is Riemann integrable on $[a,b]$.
The main idea behind this is that the difference between the maximum
and minimum values over a partition gets small. That is if
$[x_{i-1}, x_i]$ is like $1/n$ is length, then the difference between
the maximum of $f$ over this interval, $M$, and the minimum, $m$ over
this interval will go to zero as $n$ gets big. That $m$ and $M$ exists
is due to the extreme value theorem, that this difference goes to $0$
is a consequence of continuity. What is needed is that this value goes
to $0$ at the same rate -- no matter what interval is being discussed --
is a consequence of a notion of uniform continuity, a concept
discussed in advanced calculus, but which holds for continuous
functions on closed intervals. Armed with this, the Riemann sum for a
general partition can be bounded by this difference times $b-a$, which
will go to zero. So the upper and lower Riemann sums will converge to
the same value.
* A "jump", or discontinuity of the first kind, is a value $c$ in $[a,b]$ where $\lim_{x \rightarrow c+} f(x)$ and $\lim_{x \rightarrow c-}f(x)$ both exist, but are not equal. It is true that a function that is not continuous on $I=[a,b]$, but only has discontinuities of the first kind on $I$ will be Riemann integrable on $I$.
For example, the function $f(x) = 1$ for $x$ in $[0,1]$ and $0$
otherwise will be integrable, as it is continuous at all but two
points, $0$ and $1$, where it jumps.
* Some functions can have infinitely many points of discontinuity and still be integrable. The example of $f(x) = 1/q$ when $x=p/q$ is rational, and $0$ otherwise is often used as an example.
## Numeric integration
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 compute $\int_0^1 x^2 dx$.
Let us fix some values:
```julia;
a, b = 0, 1
f(x) = x^2
```
Then for a given $n$ we have some steps to do: create the partition,
find the $c_i$, multiply the pieces and add up. Here is one way to do all this:
```julia;
n = 5
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.
```
Now to multiply the values. We want to sum the product `f(cs[i]) * deltas[i]`, here is one way to do so:
```julia;
sum(f(cs[i]) * deltas[i] for i in 1:length(deltas))
```
Our 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:
```julia; hold=true
n = 50_000
xs = a:(b-a)/n:b
deltas = diff(xs)
cs = xs[1:end-1]
sum(f(cs[i]) * deltas[i] for i in 1:length(deltas))
```
This 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 the Riemann sum for
us with an extra argument to specifying one of four methods for
computing $c_i$:
```julia; eval=false
function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
if method == "right"
meth = f -> (lr -> begin l,r = lr; f(r) * (r-l) end)
elseif method == "left"
meth = f -> (lr -> begin l,r = lr; f(l) * (r-l) end)
elseif method == "trapezoid"
meth = f -> (lr -> begin l,r = lr; (1/2) * (f(l) + f(r)) * (r-l) end)
elseif method == "simpsons"
meth = f -> (lr -> begin l,r=lr; (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l) end)
(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. 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`):
```julia;
𝒇(x) = exp(x)
riemann(𝒇, 0, 5, 10) # S_10
```
Or with more intervals in the partition
```julia;
riemann(𝒇, 0, 5, 50_000)
```
(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.)
## "Negative" area
So far, we have had the assumption that $f(x) \geq 0$, as that allows
us to define the concept of area. We can define the signed area
between $f(x)$ and the $x$ axis through the definite integral:
```math
A = \int_a^b f(x) dx.
```
The right hand side is defined whenever the Riemann limit exists and
in that case we call $f(x)$ Riemann integrable. (The definition does
not suppose $f$ is non-negative.)
Suppose $f(a) = f(b) = 0$ for $a < b$ and for all $a < x < b$ we have $f(x) < 0$. Then we can see easily from the geometry (or from the Riemann sum approximation) that
If we think of the area below the $x$ axis as "signed" area carrying a minus sign, then the total area can be seen again as a sum, only this time some of the summands may be negative.
##### Example
Consider a function $g(x)$ defined through its piecewise linear graph:
```julia; echo=false
g(x) = abs(x) > 2 ? 1.0 : abs(x) - 1.0
plot(g, -3,3)
plot!(zero)
```
* Compute $\int_{-3}^{-1} g(x) dx$. The area comprised of a square of area $1$ and a triangle with area $1/2$, so should be $3/2$.
* Compute $\int_{-3}^{0} g(x) dx$. In addition to the above, there is a triangle with area $1/2$, but since the function is negative, this area is added in as $-1/2$. In total then we have $1 + 1/2 - 1/2 = 1$ for the answer.
* Compute $\int_{-3}^{1} g(x) dx$:
We could add the signed area over $[0,1]$ to the above, but instead see a square of area $1$, a triangle with area $1/2$ and a triangle with signed area $-1$. The total is then $1/2$.
* Compute $\int_{-3}^{3} g(x) dx$:
We could add the area, but let's use a symmetry trick. This is clearly
twice our second answer, or $2$. (This is because $g(x)$ is an even
function, as we can tell from the graph.)
##### Example
Suppose $f(x)$ is an odd function, then $f(x) = - f(-x)$ for any
$x$. So the signed area between $[-a,0]$ is related to the signed area
between $[0,a]$ but of different sign. This gives $\int_{-a}^a f(x) dx
= 0$ for odd functions.
An immediate consequence would be $\int_{-\pi}^\pi \sin(x) = 0$, as would $\int_{-a}^a x^k dx$ for any *odd* integer $k > 0$.
##### Example
Numerically estimate the definite integral $\int_0^e x\log(x) dx$. (We
redefine the function to be $0$ at $0$, so it is continuous.)
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:
```julia;
𝒉(x) = x > 0 ? x * log(x) : 0.0
```
This is actually inefficient, as the check for the size of `x` will
slow things down a bit. Since we will call this function 50,000 times, we
would like to avoid this, if we can. In this case just using the right
sum will work:
```julia;
h(x) = x * log(x)
riemann(h, 0, 2, 50_000, method="right")
```
(The default is `"right"`, so no method specified would also work.)
##### Example
Let $j(x) = \sqrt{1 - x^2}$. The area under the curve between $-1$ and
$1$ is $\pi/2$. Using a Riemann sum with 4 equal subintervals and the
midpoint, estimate $\pi$. How close are you?
The partition is $-1 < -1/2 < 0 < 1/2 < 1$. The midpoints are $-3/4,
-1/4, 1/4, 3/4$. We thus have that $\pi/2$ is approximately:
```julia; hold=true
xs = range(-1, 1, length=5)
deltas = diff(xs)
cs = [-3/4, -1/4, 1/4, 3/4]
j(x) = sqrt(1 - x^2)
a = sum(j(c)*delta for (c,delta) in zip(cs, deltas))
a, pi/2 # π ≈ 2a
```
(For variety, we used an alternate way to sum over two vectors.)
So $\pi$ is about `2a`.
##### Example
We have the well-known triangle
[inequality](http://en.wikipedia.org/wiki/Triangle_inequality) which
says for an individual sum:
$\lvert a + b \rvert \leq \lvert a \rvert +\lvert b \rvert$.
This can be used to give bounds on the size of an integral. For example, suppose you know that $f(x)$ is continuous on $[a,b]$ and takes its maximum value of $M$ and minimum value of $m$. Letting $K$ be the larger of $\lvert M\rvert$ and $\lvert m \rvert$, gives this bound when $a < b$:
The formula for Simpson's rule was the *composite* formula. If just a single rectangle is approximated over $[a,b]$ by a parabola interpolating the points $x_1=a$, $x_2=(a+b)/2$, and $x_3=b$, the formula is:
```math
\frac{b-a}{6}(f(x_1) + 4f(x_2) + f(x_3)).
```
This formula will actually be exact for any 3rd degree polynomial. In
fact an entire family of similar approximations using $n$ points can
be made exact for any polynomial of degree $n-1$ or lower. But with
non-evenly spaced points, even better results can be found.
The formulas for an approximation to the integral $\int_{-1}^1 f(x)
*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$, provided $f$ can
be well approximated by a polynomial over $[-1,1]$. (Which is the case
for the "nice" functions we encounter.) Some examples are given in the questions.
### The quadgk function
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, we have:
```julia; hold=true
f(x) = x * log(x)
quadgk(f, 0, 2)
```
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 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:
```julia; hold=true
f(x) = x^5 - x + 1
quadgk(f, -2, 2)
```
The error term is $0$, answer is $4$ up to the last unit of precision
(1 ulp), so any error is only in floating point approximations.
For the numeric computation of definite integrals, the `quadgk`
function should be used over the Riemann sums or even Simpson's rule.
Here are some sample integrals computed with `quadgk`:
```math
\int_0^\pi \sin(x) dx
```
```julia;
quadgk(sin, 0, pi)
```
(Again, the actual answer is off only in the last digit, the error estimate is an upper bound.)
```math
\int_0^2 x^x dx
```
```julia;
quadgk(x -> x^x, 0, 2)
```
```math
\int_0^5 e^x dx
```
```julia;
quadgk(exp, 0, 5)
```
When composing the answer with other functions it may be desirable to
drop the error in the answer. Two styles can be used for this. The
first is to just name the two returned values:
```julia; hold=true
A, err = quadgk(cos, 0, pi/4)
A
```
The second is to ask for just the first component of the returned value:
```julia; hold=true
A = quadgk(tan, 0, pi/4)[1] # or first(quadgk(tan, 0, pi/4))
```
----
To visualize the choice of nodes by the algorithm, we have for ``f(x)=\sin(x)`` over ``[0,\pi]`` relatively few nodes used to get a high-precision estimate:
```julia; echo=false
function FnWrapper(f)
xs=Any[]
ys=Any[]
x -> begin
fx = f(x)
push!(xs, x)
push!(ys, fx)
fx
end
end
nothing
```
```julia; hold=true; echo=false
let
a, b= 0, pi
f(x) = sin(x)
F = FnWrapper(f)
ans,err = quadgk(F, a, b)
plot(f, a, b, legend=false, title="Error ≈ $(round(err,sigdigits=2))")
scatter!(F.xs, F.ys)
end
```
For a more oscillatory function, more nodes are chosen:
```julia; hold=true; echo=false
let
a, b= 0, pi
f(x) = exp(-x)*sinpi(x)
F = FnWrapper(f)
ans,err = quadgk(F, a, b)
plot(f, a, b, legend=false, title="Error ≈ $(round(err,sigdigits=2))")
scatter!(F.xs, F.ys)
end
```
##### Example
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$
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
```julia;
k(x) = exp(cos(x))
A,err = quadgk(k, -pi, pi)
```
So
```julia;
C = 1/A
k₁(x) = C * exp(cos(x))
```
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)$:
```julia;
K(x) = quadgk(k₁, -pi, x)[1]
```
(The trailing `[1]` is so only the answer - and not the error - is returned.)
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 `[find_zero(x -> K(x) - p, (-pi, pi)) for p in [0.25, 0.5, 0.75]]`, but that is a bit less performant than using the `solve` interface for this task:
```julia; hold=true
Z = ZeroProblem((x,p) -> K(x) - p, (-pi, pi))
solve.(Z, (1/4, 1/2, 3/4))
```
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$.
## Questions
###### Question
Using geometry, compute the definite integral:
```math
\int_{-5}^5 \sqrt{5^2 - x^2} dx.
```
```julia; hold=true; echo=false
f(x) = sqrt(5^2 - x^2)
val, _ = quadgk(f, -5,5)
numericq(val)
```
###### Question
Using geometry, compute the definite integral:
```math
\int_{-2}^2 (2 - \lvert x\rvert) dx
```
```julia; hold=true; echo=false
f(x) = 2- abs(x)
a,b = -2, 2
val, _ = quadgk(f, a,b)
numericq(val)
```
###### Question
Using geometry, compute the definite integral:
```math
\int_0^3 3 dx + \int_3^9 (3 + 3(x-3)) dx
```
```julia; hold=true; echo=false
f(x) = x <= 3 ? 3.0 : 3 + 3*(x-3)
a,b = 0, 9
val, _ = quadgk(f, a, b)
numericq(val)
```
###### Question
Using geometry, compute the definite integral:
```math
\int_0^5 \lfloor x \rfloor dx
```
(The notation $\lfloor x \rfloor$ is the integer such that $\lfloor x \rfloor \leq x < \lfloor x \rfloor + 1$.)
```julia; hold=true; echo=false
f(x) = floor(x)
a, b = 0, 5
val, _ = quadgk(f, a, b)
numericq(val)
```
###### Question
Using geometry, compute the definite integral between $-3$ and $3$ of this graph comprised of lines and circular arcs:
```julia; hold=true; echo=false
function f(x)
if x < -1
abs(x+1)
elseif -1 <= x <= 1
sqrt(1 - x^2)
else
abs(x-1)
end
end
plot(f, -3, 3, aspect_ratio=:equal)
```
The value is:
```julia; hold=true; echo=false
val = (1/2 * 2 * 2) * 2 + pi*1^2/2
numericq(val)
```
###### Question
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$.
```julia; hold=true; echo=false
f(x) = sin(x)
xs = -1:1/2:1
deltas = diff(xs)
val = sum(map(f, xs[1:end-1]) .* deltas)
numericq(val)
```
###### Question
Without doing any *real* work, find this integral:
```math
\int_{-\pi/4}^{\pi/4} \tan(x) dx.
```
```julia; hold=true; echo=false
val = 0
numericq(val)
```
###### Question
Without doing any *real* work, find this integral:
```math
\int_3^5 (1 - \lvert x-4 \rvert) dx
```
```julia; hold=true; echo=false
val = 1
numericq(val)
```
###### Question
Suppose you know that for the integrable function
$\int_a^b f(u)du =1$ and $\int_a^c f(u)du = p$. If $a < c < b$ what is
