Filon quadrature

<!--

-->
 UoN Applied & Industrial Maths | 2025 Oct 17

# Filon Quadrature

### Dave Smith

Quadrature examples at 
 <a href="https://www.desmos.com/calculator/hdvosguhh4">www.desmos.com/calculator/hdvosguhh4</a>
 
 
 ---

## Quadrature
            **Aim**: evaluate an integral numerically

Important because analytical evaluation of integrals is usually impossible.

eg. if
 $\displaystyle f(x) = \cos\left(\frac{x^2}{50}+\frac x{10}\right) + \frac{2x}5 + 2$
 then
 $\displaystyle \int_a^b f(x) dx$
 is difficult analytically.
 
 ---

## Approximation
            Approximate the function with something easy to integrate analytically.

### Constant approximation
            Define
            $\displaystyle f_{\mathrm{rect}}(x) = f_{\mathrm{rect}} := f(c)$
            for some $c\in\[a,b\]$.

Then
            $\displaystyle f_{\mathrm{rect}}(x) \approx f(x)$.

Therefore,
            $\displaystyle \int_a^b f(x) dx \approx (b-a)f_{\mathrm{rect}}$.

But how do we choose $c$?
            No clear answer, though IVT can be used to show $\exists$ $c\in[a,b]$ such that error is zero.

### Error
            $\displaystyle \mathcal E_{\mathrm{rect}} \leq (b-a) \left(\max_{\xi\in[a,b]}\left\lvert f(\xi)\right\rvert - \min_{\xi\in[a,b]}\left\lvert f(\xi)\right\rvert \right)$.

By MVT,

$\displaystyle \mathcal E_{\mathrm{rect}} \leq (b-a) \max_{\xi\in[a,b]}\left\lvert f'(\xi)\right\rvert$.
 
 ---

## Approximation
            Approximate the function with something easy to integrate analytically.

### Piecewise constant approximation
            Let $m=(a+b)/2$ be the midpoint of the interval.

Define
 $\displaystyle f_{\mathrm{rect}2}(x) := \begin{cases} f(c_1) & \text{if }a \leq x \leq m \\\\ f(c_2) & \text{if }m \leq x \leq b.\end{cases}$
 
 $\exists$ $c_1\in\[a,m\]$, $c_2\in\[m,b\]$.
 Then
 $\displaystyle f_{\mathrm{rect}2}(x) \approx f(x)$.

Therefore,

$\displaystyle \int_a^b f(x) dx \approx \frac{b-a}2(f(c_1)+f(c_2))$.

But how do we choose $c_j$?
            Still no clear answer.

### Error
 $\displaystyle \mathcal E_{\mathrm{rect}2} \leq (m-a) \left(\max_{\xi\in[a,m]}\left\lvert f(\xi)\right\rvert - \min_{\xi\in[a,m]}\left\lvert f(\xi)\right\rvert \right) + (b-m) \left(\max_{\xi\in[m,b]}\left\lvert f(\xi)\right\rvert - \min_{\xi\in[m,b]}\left\lvert f(\xi)\right\rvert \right) \leq \mathcal E_{\mathrm{rect}}$.
 
 ---

## Approximation
            Approximate the function with something easy to integrate analytically.

### Trapezium rule
            Approximate by a straight line interpolation between $(a,f(a))$ and $(b,f(b))$.

Then
            $\displaystyle f_{\mathrm{tpz}}(x) \approx f(x)$.

Therefore,

$\displaystyle \int_a^b f(x) dx \approx (b-a)\frac{f(b)+f(a)}2$.

Note: This is the same integral as the $f_{\mathrm{rect}2}$ approximation with $c_1=a$ and $c_2=b$.
 
 ---

## Approximation
            Approximate the function with something easy to integrate analytically.

### Trapezium rule
            Approximate by a straight line interpolation between $(a,f(a))$ and $(b,f(b))$.

### Multiple trapezia
            Approximate by a piecewise straight line interpolation between $(a,f(a))$, $(m,f(m))$, and $(b,f(b))$.
            Then

$$\begin{align*}
            \int_a^b f(x) dx \&\approx \frac{b-a}2 \times \frac{f(a)+2f(m)+f(b)}2 \\\\
            \&= \frac h2 \left(f(a)+2f(m)+f(b)\right),
            \end{align*}$$
            where
            $h = \frac{b-a}2$.

With $N$ equal subdivisions
 $$
 \int_a^b f(x) dx \approx \frac h2 \left(f(a)+2\sum\_{j=1}\^{N-1}f(m_j)+f(b)\right).
 $$
 
 ---

## Approximation
            Approximate the function with something easy to integrate analytically.

### Multiple trapezia
            With $N$ equal subdivisions
            $$
            \int_a^b f(x) dx \approx \frac h2 \left(f(a)+2\sum\_{j=1}\^{N-1}f(m_j)+f(b)\right).
            $$

### Error
 $\displaystyle \mathcal E_{\mathrm{rect}2} \leq \frac{(b-a)h^2}{12} \max\_{\xi\in[a,b]}\left\lvert f\^{\prime\prime}(\xi) \right\rvert$.
 
 ---

## Approximation
            Approximate the function with something easy to integrate analytically.

### Multiple rectangles
            Piecewise constant approximation: polynomial order $0$.

Error $\mathcal E\_{\mathrm{rect}}\lesssim (b-a)h\left\lVert f' \right\rVert_\infty$ for $h$ the subinterval width.

### Multiple trapezia
            Piecewise affine approximation: polynomial order $1$.

Error $\mathcal E\_{\mathrm{tpz}}\lesssim (b-a)h^2\left\lVert f'' \right\rVert_\infty$.

### Next step?
            Piecewise quadratic approximation: polynomial order $2$.

Called *Simpson's rule*.

Hypothesis: error $\mathcal E\_{\mathrm{Smp}}\lesssim (b-a)h^3\left\lVert f''' \right\rVert_\infty$.
 
 ---

## Simpson's rule
            Approximate $f$ by a quadratic, exactly at the endpoints and the midpoint.

### Quadratic approximation
            Let
            $f_{\mathrm{Smp}}(x) = \rho x^2 + \sigma x + \tau \approx f(x)$,
            with
            $$
            \begin{pmatrix} a^2 \& a \& 1 \\\\ m^2 \& m \& 1 \\\\ b^2 \& b \& 1 \end{pmatrix} \begin{pmatrix} \rho \\\\ \sigma \\\\ \tau \end{pmatrix} = \begin{pmatrix} f(a) \\\\ f(m) \\\\ f(b) \end{pmatrix}
            $$

Then
            $$\begin{align*}
            \int_a^b f(x) dx \&\approx \frac \rho3(a^3-b^3) + \frac \sigma2(a^2-b^2) + \tau(a-b) \\\\
            \&= \frac{b-a}3(f(a) + 4f(m) + f(b)).
            \end{align*}$$

### Error
 $\mathcal E_{\mathrm{Smp}} \lesssim (b-a)h^4\left\lVert f\^{(4)} \right\rVert_\infty$ where $h$ is the interval width $(b-a)$.
 
 ---

## Simpson's rule
            Approximate $f$ by a quadratic, exactly at the endpoints and the midpoint.

### Error
            $\mathcal E_{\mathrm{Smp}} \lesssim (b-a)h^4\left\lVert f\^{(4)} \right\rVert_\infty$
            but we expected only
            $\mathcal E_{\mathrm{Smp}} \lesssim (b-a)h^3\left\lVert f\^{\prime\prime\prime} \right\rVert_\infty$.
            Why is the error so good?
            It looks like we approximated using a cubic interpolant instead of quadratic.

Consider a cubic interpolant of these three points.
            Clearly there is one degree of freedom.
            Can write the cubic as
            $$
            f_{\mathrm{cubic}}(x) = f_{\mathrm{Smp}}(x) + C\begin{pmatrix}\text{cubic which evaluates}\\\\\text{to } 0 \text{ at } x=a,m,b.\end{pmatrix}
            $$
            for $C$ tuned to give the "best" approximation.
            So
            $$
            f_{\mathrm{cubic}}(x) = f_{\mathrm{Smp}}(x) + C(x-a)(x-m)(x-b)
            $$

Then
            $\displaystyle \int_a^b f_{\mathrm{cubic}}(x) dx = \int_a^b f_{\mathrm{Smp}}(x) dx$,
            because the cubic correction is odd on the interval.

Good to use even polynomial approximants with symmetric sample point distribution; get extra degree of error for free.

---

## Simpson's rule
            Approximate $f$ by a quadratic, exactly at the endpoints and the midpoint.

But now make these quadratic approximations on each of the $N$ subintervals $[a,m_2],[m_2,m_4],\ldots,[m_{2N-2},b]$.
            The midpoint of each subinterval is then $m_{2j-1}$.

### Composite quadratic approximation
 $\displaystyle \int_a^b f(x) dx \approx \int_a^b f_{\mathrm{cSmp}}(x) dx =$
 $$
 \frac{h}3\left(f(a) + 4\sum\_{j=1}\^{N}f(m_\{2j-1}) + 2\sum\_{j=1}\^{N-1}f(m_\{2j}) + f(b)\right)
 $$
 $\displaystyle =h\left(\frac23S_{\mathrm{even}} + \frac43S_{\mathrm{odd}}\right)$.

### Error in composite Simpson's rule
 $\mathcal E_{\mathrm{cSmp}} \lesssim (b-a)h^4\left\lVert f\^{(4)} \right\rVert_\infty$, where $h$ is the subinterval width.
 
 ---

## When is Simpson quadrature expensive?
 If $\left\lVert f\^{(4)} \right\rVert_\infty$ is large, we need $h$ small so that 
 $\mathcal E_{\mathrm{cSmp}} \lesssim (b-a)h^4\left\lVert f\^{(4)} \right\rVert_\infty$
 is small.

$\left\lVert f\^{(4)} \right\rVert_\infty$ is large in highly oscillatory functions, such as
            $$
            \left(\cos\left(\frac{x^2}{50}+\frac x{10}\right) + \frac{2x}5 + 2\right)\left(\frac7{10}+\frac3{10}\sin\left(kx\right)\right),
            $$
            with $k\gg b-a$.

It is not just Simpson quadrature; they are all bad unless $h$ is very small (ie many subintervals).

We need a new approach.
 
 ---

## Integrals with trigonometric oscillation
            Suppose $f$ is well approximated by its quadratic interpolant, $k\gg b-a$.
            Aim to evaluate numerically

$$
            \int_a^b f(x) \left(\frac7{10}+\frac3{10}\sin\left(kx\right)\right) dx
            $$
            $$
            = \frac7{10}\int_a^b f(x) dx + \frac3{10}\int_a^b f(x) \sin\left(kx\right) dx;
            $$
            first integral we can use Simpson's rule.

So we need only evaluate

$$
 \int_a^b f(x) \sin\left(kx\right) dx.
 $$
 
 ---

## Filon quadrature
            Suppose $f$ is well approximated by its quadratic interpolant, $k\gg b-a$.
            Aim to evaluate numerically

$$
            \int_a^b f(x) \sin\left(kx\right) dx.
            $$

Approximate $f$ by a piecewise quadratic over $N$ subintervals width $h$.
            Then integral can be explicitly calculated.

Doing so, and after some calculation,

$\displaystyle \int_a^b f(x) \sin\left(kx\right) dx$
 $$
 \approx h\Big( \alpha [f(a)\cos(ka) - f(b)\cos(kb)]
 $$
 $\displaystyle + \beta S_{\mathrm{even}} + \gamma S_{\mathrm{odd}} \Big)$,

where

$$\begin{align*}
 \scriptstyle S_{\mathrm{even}} \&\scriptstyle = \scriptstyle\sum\_{j=1}\^{N-1}f(m_\{2j})\sin km_\{2j} + \frac{f(a)\sin ka + f(b) \sin kb}2, \\\\
 \scriptstyle S_{\mathrm{odd}} \&\scriptstyle = \scriptstyle\sum\_{j=1}\^{N}f(m_\{2j-1})\sin km_\{2j-1}.
 \end{align*}$$
 
 ---

## Filon quadrature
            Approximate $f$ by a piecewise quadratic over $N$ subintervals width $h$.
            Then integral can be explicitly calculated.

where

$$\begin{align*}
 \scriptstyle \alpha \&\scriptstyle = \scriptstyle \frac1{kh} + \frac{\cos kh \sin kh}{(kh)^2} - \frac{2\sin^2kh}{(kh)^3}, \\\\
 \scriptstyle \beta \&\scriptstyle = \scriptstyle \frac{2(1+\cos^2kh)}{(kh)^2} - \frac{4\sin kh\cos kh}{(kh)^3}, \\\\
 \scriptstyle \gamma \&\scriptstyle = \scriptstyle \frac{-4\cos kh}{(kh)^2} + \frac{4\sin kh}{(kh)^3}, \\\\
 \scriptstyle S_{\mathrm{even}} &\scriptstyle = \scriptstyle \sum\_{j=1}\^{N-1}f(m_\{2j})\sin km_\{2j} + \frac{f(a)\sin ka + f(b) \sin kb}2, \\\\
 \scriptstyle S_{\mathrm{odd}} &\scriptstyle = \scriptstyle \sum\_{j=1}\^{N}f(m_\{2j-1})\sin km_\{2j-1}.
 \end{align*}$$
 
 ---

## Filon quadrature

where

Here $\alpha,\beta,\gamma$ are given as *exact* functions of $kh$.

Terms $(kh)^{-1},(kh)^{-2},(kh)^{-3}$ appear because we integrate $\sin(kx)$ thrice (by parts twice).
 
 ---

## Recovering Simpson's rule

$\scriptstyle \int_a^b f(x) \sin\left(kx\right) dx \approx h\big( \alpha [f(a)\cos(ka) - f(b)\cos(kb)] + \beta S_{\mathrm{even}} + \gamma S_{\mathrm{odd}} \big)$,

where

Here $\alpha,\beta,\gamma$ are given as *exact* functions of $kh$.

More subintervals: $h\to0^+$ so $kh\to0^+$.
            Small $kh$ approximation of coefficients:
            $$\begin{align*}
            \scriptstyle \alpha \&\scriptstyle = \scriptstyle (kh)^3\frac2{45} - \mathcal O \left((kh)^5\right), \\\\
            \scriptstyle \beta \&\scriptstyle = \scriptstyle \frac23 + (kh)^2\frac2{15} - \mathcal O \left((kh)^4\right), \\\\
            \scriptstyle \gamma \&\scriptstyle = \scriptstyle \frac43 - (kh)^2\frac2{15} - \mathcal O \left((kh)^4\right).
            \end{align*}$$

So
            $\scriptstyle \int_a^b f(x) \sin\left(kx\right) dx \approx h\big( \frac23 S_{\mathrm{even}} + \frac43 S_{\mathrm{odd}} \big)$
            for $h$ small.

But we should not need so many subintervals.
 
 ---

# Thanks!

### These slides

## ![qr code for URL of slides](q/filon-quadrature.svg)
 
 ---

## Further topics

* Error bounds for Filon quadrature
            * Adaptive Simpson quadrature
            * Adaptive Filon quadrature
            * Filon quadrature on a semiinfinite interval
            * Filon-Curtis-Clenshaw quadrature