- §3.11(i) Minimax Polynomial Approximations
- §3.11(ii) Chebyshev-Series Expansions
- §3.11(iii) Minimax Rational Approximations
- §3.11(iv) Padé Approximations
- §3.11(v) Least Squares Approximations
- §3.11(vi) Splines

Let $f(x)$ be continuous on a closed interval $[a,b]$. Then there exists a unique $n$th degree polynomial ${p}_{n}(x)$, called the
*minimax* (or *best uniform*) polynomial approximation to $f(x)$
on $[a,b]$, that
minimizes ${\mathrm{max}}_{a\le x\le b}\left|{\u03f5}_{n}(x)\right|$, where
${\u03f5}_{n}(x)=f(x)-{p}_{n}(x)$.

A sufficient condition for ${p}_{n}(x)$ to be the minimax polynomial is that $\left|{\u03f5}_{n}(x)\right|$ attains its maximum at $n+2$ distinct points in $[a,b]$ and ${\u03f5}_{n}(x)$ changes sign at these consecutive maxima.

If we have a sufficiently close approximation

3.11.1 | $${p}_{n}(x)={a}_{n}{x}^{n}+{a}_{n-1}{x}^{n-1}+\mathrm{\cdots}+{a}_{0}$$ | ||

to $f(x)$, then the coefficients ${a}_{k}$ can be computed iteratively. Assume that ${f}^{\prime}(x)$ is continuous on $[a,b]$ and let ${x}_{0}=a$, ${x}_{n+1}=b$, and ${x}_{1},{x}_{2},\mathrm{\dots},{x}_{n}$ be the zeros of ${\u03f5}_{n}^{\prime}(x)$ in $(a,b)$ arranged so that

3.11.2 | $$ | ||

Also, let

3.11.3 | $${m}_{j}={(-1)}^{j}{\u03f5}_{n}({x}_{j}),$$ | ||

$j=0,1,\mathrm{\dots},n+1$. | |||

(Thus the ${m}_{j}$ are approximations to $m$, where $\pm m$ is the maximum value of $\left|{\u03f5}_{n}(x)\right|$ on $[a,b]$.)

Then (in general) a better approximation to ${p}_{n}(x)$ is given by

3.11.4 | $$\sum _{k=0}^{n}({a}_{k}+\mathit{\delta}{a}_{k}){x}^{k},$$ | ||

where

3.11.5 | $$\sum _{k=0}^{n}{x}_{j}^{k}\mathit{\delta}{a}_{k}={(-1)}^{j}({m}_{j}-m),$$ | ||

$j=0,1,\mathrm{\dots},n+1$. | |||

This is a set of $n+2$ equations for the $n+2$ unknowns $\mathit{\delta}{a}_{0},\mathit{\delta}{a}_{1},\mathrm{\dots},\mathit{\delta}{a}_{n}$ and $m$.

The iterative process converges locally and quadratically (§3.8(i)).

A method for obtaining a sufficiently accurate first approximation is described in the next subsection.

The Chebyshev polynomials ${T}_{n}$ are given by

3.11.6 | $${T}_{n}\left(x\right)=\mathrm{cos}\left(n\mathrm{arccos}x\right),$$ | ||

$-1\le x\le 1$. | |||

They satisfy the recurrence relation

3.11.7 | $${T}_{n+1}\left(x\right)-2x{T}_{n}\left(x\right)+{T}_{n-1}\left(x\right)=0,$$ | ||

$n=1,2,\mathrm{\dots}$, | |||

with initial values ${T}_{0}\left(x\right)=1$, ${T}_{1}\left(x\right)=x$. They enjoy an orthogonal property with respect to integrals:

3.11.8 | $${\int}_{-1}^{1}\frac{{T}_{j}\left(x\right){T}_{k}\left(x\right)}{\sqrt{1-{x}^{2}}}dx=\{\begin{array}{cc}\pi ,\hfill & j=k=0,\hfill \\ \frac{1}{2}\pi ,\hfill & j=k\ne 0,\hfill \\ 0,\hfill & j\ne k,\hfill \end{array}$$ | ||

as well as an orthogonal property with respect to sums, as follows. When $n>0$ and $0\le j\le n$, $0\le k\le n$,

3.11.9 | $$\underset{\mathrm{\ell}=0\phantom{\prime \prime}}{\overset{n\phantom{\prime \prime}}{{\sum}^{\prime \prime}}}{T}_{j}\left({x}_{\mathrm{\ell}}\right){T}_{k}\left({x}_{\mathrm{\ell}}\right)=\{\begin{array}{cc}n,\hfill & j=k=0\text{or}n,\hfill \\ \frac{1}{2}n,\hfill & j=k\ne 0\text{or}n,\hfill \\ 0,\hfill & j\ne k,\hfill \end{array}$$ | ||

where ${x}_{\mathrm{\ell}}=\mathrm{cos}\left(\pi \mathrm{\ell}/n\right)$ and the double prime means that the first and last terms are to be halved.

For these and further properties of Chebyshev polynomials, see Chapter 18, Gil et al. (2007a, Chapter 3), and Mason and Handscomb (2003).

If $f$ is continuously differentiable on $[-1,1]$, then with

3.11.10 | $${c}_{n}=\frac{2}{\pi}{\int}_{0}^{\pi}f(\mathrm{cos}\theta )\mathrm{cos}\left(n\theta \right)d\theta ,$$ | ||

$n=0,1,2,\mathrm{\dots}$, | |||

the expansion

3.11.11 | $$f(x)=\underset{n=0\phantom{\prime}}{\overset{\mathrm{\infty}\phantom{\prime}}{{\sum}^{\prime}}}{c}_{n}{T}_{n}\left(x\right),$$ | ||

$-1\le x\le 1$, | |||

converges uniformly. Here the single prime on the summation symbol means that the first term is to be halved. In fact, (3.11.11) is the Fourier-series expansion of $f(\mathrm{cos}\theta )$; compare (3.11.6) and §1.8(i).

Furthermore, if $f\in {C}^{\mathrm{\infty}}[-1,1]$, then the convergence of (3.11.11) is usually very rapid; compare (1.8.7) with $k$ arbitrary.

For general intervals $[a,b]$ we rescale:

3.11.12 | $$f(x)=\underset{n=0\phantom{\prime}}{\overset{\mathrm{\infty}\phantom{\prime}}{{\sum}^{\prime}}}{d}_{n}{T}_{n}\left(\frac{2x-a-b}{b-a}\right).$$ | ||

Because the series (3.11.12) converges rapidly we obtain a very good first approximation to the minimax polynomial ${p}_{n}(x)$ for $[a,b]$ if we truncate (3.11.12) at its $(n+1)$th term. This is because in the notation of §3.11(i)

3.11.13 | $${\u03f5}_{n}(x)={d}_{n+1}{T}_{n+1}\left(\frac{2x-a-b}{b-a}\right),$$ | ||

approximately, and the right-hand side enjoys exactly those properties concerning its maxima and minima that are required for the minimax approximation; compare Figure 18.4.3.

More precisely, it is known that for the interval $[a,b]$, the ratio of the maximum value of the remainder

3.11.14 | $$\left|\sum _{k=n+1}^{\mathrm{\infty}}{d}_{k}{T}_{k}\left(\frac{2x-a-b}{b-a}\right)\right|$$ | ||

to the maximum error of the minimax polynomial ${p}_{n}(x)$ is bounded by $1+{L}_{n}$,
where ${L}_{n}$ is the $n$th *Lebesgue constant* for Fourier series; see
§1.8(i). Since ${L}_{0}=1$, ${L}_{n}$ is a monotonically increasing
function of $n$, and (for example) ${L}_{1000}=4.07\mathrm{\dots}$, this means that in
practice the gain in replacing a truncated Chebyshev-series expansion by the
corresponding minimax polynomial approximation is hardly worthwhile.
Moreover, the set of minimax approximations
${p}_{0}(x),{p}_{1}(x),{p}_{2}(x),\mathrm{\dots},{p}_{n}(x)$ requires the calculation and storage of
$\frac{1}{2}(n+1)(n+2)$ coefficients, whereas the corresponding set of
Chebyshev-series approximations requires only $n+1$ coefficients.

The ${c}_{n}$ in (3.11.11) can be calculated from (3.11.10), but in general it is more efficient to make use of the orthogonal property (3.11.9). Also, in cases where $f(x)$ satisfies a linear ordinary differential equation with polynomial coefficients, the expansion (3.11.11) can be substituted in the differential equation to yield a recurrence relation satisfied by the ${c}_{n}$.

For the expansion (3.11.11), numerical values of the Chebyshev polynomials ${T}_{n}\left(x\right)$ can be generated by application of the recurrence relation (3.11.7). A more efficient procedure is as follows. Let ${c}_{n}{T}_{n}\left(x\right)$ be the last term retained in the truncated series. Beginning with ${u}_{n+1}=0$, ${u}_{n}={c}_{n}$, we apply

3.11.15 | $${u}_{k}=2x{u}_{k+1}-{u}_{k+2}+{c}_{k},$$ | ||

$k=n-1,n-2,\mathrm{\dots},0$. | |||

Then the sum of the truncated expansion equals $\frac{1}{2}({u}_{0}-{u}_{2})$. For error analysis and modifications of Clenshaw’s algorithm, see Oliver (1977).

If $x$ is replaced by a complex variable $z$ and $f(z)$ is analytic, then the expansion (3.11.11) converges within an ellipse. However, in general (3.11.11) affords no advantage in $\u2102$ for numerical purposes compared with the Maclaurin expansion of $f(z)$.

For further details on Chebyshev-series expansions in the complex plane, see Mason and Handscomb (2003, §5.10).

Let $f$ be continuous on a closed interval $[a,b]$ and $w$ be a continuous
nonvanishing function on $[a,b]$: $w$ is called a *weight function*. Then
the *minimax* (or *best uniform*) rational approximation

3.11.16 | $${R}_{k,\mathrm{\ell}}(x)=\frac{{p}_{0}+{p}_{1}x+\mathrm{\cdots}+{p}_{k}{x}^{k}}{1+{q}_{1}x+\mathrm{\cdots}+{q}_{\mathrm{\ell}}{x}^{\mathrm{\ell}}}$$ | ||

of *type* $[k,\mathrm{\ell}]$ to $f$ on $[a,b]$ minimizes the maximum value of
$\left|{\u03f5}_{k,\mathrm{\ell}}(x)\right|$ on $[a,b]$, where

3.11.17 | $${\u03f5}_{k,\mathrm{\ell}}(x)=\frac{{R}_{k,\mathrm{\ell}}(x)-f(x)}{w(x)}.$$ | ||

The theory of polynomial minimax approximation given in §3.11(i) can be extended to the case when ${p}_{n}(x)$ is replaced by a rational function ${R}_{k,\mathrm{\ell}}(x)$. There exists a unique solution of this minimax problem and there are at least $k+\mathrm{\ell}+2$ values ${x}_{j}$, $$, such that ${m}_{j}=m$, where

3.11.18 | $${m}_{j}={(-1)}^{j}{\u03f5}_{k,\mathrm{\ell}}({x}_{j}),$$ | ||

$j=0,1,\mathrm{\dots},k+\mathrm{\ell}+1$, | |||

and $\pm m$ is the maximum of $\left|{\u03f5}_{k,\mathrm{\ell}}(x)\right|$ on $[a,b]$.

A collection of minimax rational approximations to elementary and special functions can be found in Hart et al. (1968).

A widely implemented and used algorithm for calculating the coefficients
${p}_{j}$ and ${q}_{j}$ in (3.11.16) is *Remez’s second algorithm*.
See Remez (1957), Werner et al. (1967), and
Johnson and Blair (1973).

With $w(x)=1$ and 14-digit computation, we obtain the following rational approximation of type $[3,3]$ to the Bessel function ${J}_{0}\left(x\right)$ (§10.2(ii)) on the interval $0\le x\le {j}_{0,1}$, where ${j}_{0,1}$ is the first positive zero of ${J}_{0}\left(x\right)$:

3.11.19 | $${R}_{3,3}(x)=\frac{{p}_{0}+{p}_{1}x+{p}_{2}{x}^{2}+{p}_{3}{x}^{3}}{1+{q}_{1}x+{q}_{2}{x}^{2}+{q}_{3}{x}^{3}},$$ | ||

with coefficients given in Table 3.11.1.

$j$ | ${p}_{j}$ | ${q}_{j}$ |
---|---|---|

0 | 0.99999 99891 7854 | |

1 | $-$0.34038 93820 9347 | $-$0.34039 05233 8838 |

2 | $-$0.18915 48376 3222 | 0.06086 50162 9812 |

3 | 0.06658 31942 0166 | $-$0.01864 47680 9090 |

The error curve is shown in Figure 3.11.1.

Let

3.11.20 | $$f(z)={c}_{0}+{c}_{1}z+{c}_{2}{z}^{2}+\mathrm{\cdots}$$ | ||

be a formal power series. The rational function

3.11.21 | $$\frac{{N}_{p,q}(z)}{{D}_{p,q}(z)}=\frac{{a}_{0}+{a}_{1}z+\mathrm{\cdots}+{a}_{p}{z}^{p}}{{b}_{0}+{b}_{1}z+\mathrm{\cdots}+{b}_{q}{z}^{q}}$$ | ||

is called a *Padé approximant at zero* of $f$ if

3.11.22 | $${N}_{p,q}(z)-f(z){D}_{p,q}(z)=O\left({z}^{p+q+1}\right),$$ | ||

$z\to 0$. | |||

It is denoted by ${[p/q]}_{f}\left(z\right)$. Thus if ${b}_{0}\ne 0$, then the Maclaurin expansion of (3.11.21) agrees with (3.11.20) up to, and including, the term in ${z}^{p+q}$.

The requirement (3.11.22) implies

3.11.23 | ${a}_{0}$ | $={c}_{0}{b}_{0},$ | ||

${a}_{1}$ | $={c}_{1}{b}_{0}+{c}_{0}{b}_{1},$ | |||

$\mathrm{\vdots}$ | ||||

${a}_{p}$ | $={c}_{p}{b}_{0}+{c}_{p-1}{b}_{1}+\mathrm{\cdots}+{c}_{p-q}{b}_{q},$ | |||

$0$ | $={c}_{p+1}{b}_{0}+{c}_{p}{b}_{1}+\mathrm{\cdots}+{c}_{p-q+1}{b}_{q},$ | |||

$\mathrm{\vdots}$ | ||||

$0$ | $={c}_{p+q}{b}_{0}+{c}_{p+q-1}{b}_{1}+\mathrm{\cdots}+{c}_{p}{b}_{q},$ | |||

where ${c}_{j}=0$ if $$. With ${b}_{0}=1$, the last $q$ equations give ${b}_{1},\mathrm{\dots},{b}_{q}$ as the solution of a system of linear equations. The first $p+1$ equations then yield ${a}_{0},\mathrm{\dots},{a}_{p}$.

The array of Padé approximants

3.11.24 | $$\begin{array}{cccc}{[0/0]}_{f}& {[0/1]}_{f}& {[0/2]}_{f}& \mathrm{\cdots}\\ {[1/0]}_{f}& {[1/1]}_{f}& {[1/2]}_{f}& \mathrm{\cdots}\\ {[2/0]}_{f}& {[2/1]}_{f}& {[2/2]}_{f}& \mathrm{\cdots}\\ \mathrm{\vdots}& \mathrm{\vdots}& \mathrm{\vdots}& \mathrm{\ddots}\end{array}$$ | ||

is called a *Padé table*. Approximants with the same denominator degree
are located in the same column of the table.

For convergence results for Padé approximants, and the connection with continued fractions and Gaussian quadrature, see Baker and Graves-Morris (1996, §4.7).

The Padé approximants can be computed by *Wynn’s cross rule*. Any five
approximants arranged in the Padé table as

satisfy

3.11.25 | $${(N-C)}^{-1}+{(S-C)}^{-1}={(W-C)}^{-1}+{(E-C)}^{-1}.$$ | ||

Starting with the first column ${[n/0]}_{f}$, $n=0,1,2,\mathrm{\dots}$, and initializing the preceding column by ${[n/-1]}_{f}=\mathrm{\infty}$, $n=1,2,\mathrm{\dots}$, we can compute the lower triangular part of the table via (3.11.25). Similarly, the upper triangular part follows from the first row ${[0/n]}_{f}$, $n=0,1,2,\mathrm{\dots}$, by initializing ${[-1/n]}_{f}=0$, $n=1,2,\mathrm{\dots}$.

For the recursive computation of ${[n+k/k]}_{f}$ by Wynn’s epsilon algorithm, see (3.9.11) and the subsequent text.

Numerical inversion of the Laplace transform (§1.14(iii))

3.11.26 | $$F(s)=\mathcal{L}f\left(s\right)={\int}_{0}^{\mathrm{\infty}}{\mathrm{e}}^{-st}f(t)dt$$ | ||

requires $f={\mathcal{L}}^{-1}F$ to be obtained from numerical values of $F$. A general procedure is to approximate $F$ by a rational function $R$ (vanishing at infinity) and then approximate $f$ by $r={\mathcal{L}}^{-1}R$. When $F$ has an explicit power-series expansion a possible choice of $R$ is a Padé approximation to $F$. See Luke (1969b, §16.4) for several examples involving special functions.

Suppose a function $f(x)$ is approximated by the polynomial

3.11.27 | $${p}_{n}(x)={a}_{n}{x}^{n}+{a}_{n-1}{x}^{n-1}+\mathrm{\cdots}+{a}_{0}$$ | ||

that minimizes

3.11.28 | $$S=\sum _{j=1}^{J}{\left(f({x}_{j})-{p}_{n}({x}_{j})\right)}^{2}.$$ | ||

Here ${x}_{j}$, $j=1,2,\mathrm{\dots},J$, is a given set of distinct real points and
$J\ge n+1$. From the equations $\partial S/\partial {a}_{k}=0$, $k=0,1,\mathrm{\dots},n$, we
derive the *normal equations*

3.11.29 | $$\left[\begin{array}{cccc}{X}_{0}& {X}_{1}& \mathrm{\cdots}& {X}_{n}\\ {X}_{1}& {X}_{2}& \mathrm{\cdots}& {X}_{n+1}\\ \mathrm{\vdots}& \mathrm{\vdots}& \mathrm{\ddots}& \mathrm{\vdots}\\ {X}_{n}& {X}_{n+1}& \mathrm{\cdots}& {X}_{2n}\end{array}\right]\left[\begin{array}{c}{a}_{0}\\ {a}_{1}\\ \mathrm{\vdots}\\ {a}_{n}\end{array}\right]=\left[\begin{array}{c}{F}_{0}\\ {F}_{1}\\ \mathrm{\vdots}\\ {F}_{n}\end{array}\right],$$ | ||

where

3.11.30 | ${X}_{k}$ | $={\displaystyle \sum _{j=1}^{J}}{x}_{j}^{k},$ | ||

${F}_{k}$ | $={\displaystyle \sum _{j=1}^{J}}f({x}_{j}){x}_{j}^{k}.$ | |||

(3.11.29) is a system of $n+1$ linear equations for the coefficients ${a}_{0},{a}_{1},\mathrm{\dots},{a}_{n}$. The matrix is symmetric and positive definite, but the system is ill-conditioned when $n$ is large because the lower rows of the matrix are approximately proportional to one another. If $J=n+1$, then ${p}_{n}(x)$ is the Lagrange interpolation polynomial for the set ${x}_{1},{x}_{2},\mathrm{\dots},{x}_{J}$ (§3.3(i)).

More generally, let $f(x)$ be approximated by a linear combination

3.11.31 | $${\mathrm{\Phi}}_{n}(x)={a}_{n}{\varphi}_{n}(x)+{a}_{n-1}{\varphi}_{n-1}(x)+\mathrm{\cdots}+{a}_{0}{\varphi}_{0}(x)$$ | ||

of given functions ${\varphi}_{k}(x)$, $k=0,1,\mathrm{\dots},n$, that minimizes

3.11.32 | $$\sum _{j=1}^{J}w({x}_{j}){\left(f({x}_{j})-{\mathrm{\Phi}}_{n}({x}_{j})\right)}^{2},$$ | ||

$w(x)$ being a given positive *weight function*, and again $J\ge n+1$.
Then (3.11.29) is replaced by

3.11.33 | $$\left[\begin{array}{cccc}{X}_{00}& {X}_{01}& \mathrm{\cdots}& {X}_{0n}\\ {X}_{10}& {X}_{11}& \mathrm{\cdots}& {X}_{1n}\\ \mathrm{\vdots}& \mathrm{\vdots}& \mathrm{\ddots}& \mathrm{\vdots}\\ {X}_{n0}& {X}_{n1}& \mathrm{\cdots}& {X}_{nn}\end{array}\right]\left[\begin{array}{c}{a}_{0}\\ {a}_{1}\\ \mathrm{\vdots}\\ {a}_{n}\end{array}\right]=\left[\begin{array}{c}{F}_{0}\\ {F}_{1}\\ \mathrm{\vdots}\\ {F}_{n}\end{array}\right],$$ | ||

with

3.11.34 | $${X}_{k\mathrm{\ell}}=\sum _{j=1}^{J}w({x}_{j}){\varphi}_{k}({x}_{j}){\varphi}_{\mathrm{\ell}}({x}_{j}),$$ | ||

and

3.11.35 | $${F}_{k}=\sum _{j=1}^{J}w({x}_{j})f({x}_{j}){\varphi}_{k}({x}_{j}).$$ | ||

Since ${X}_{k\mathrm{\ell}}={X}_{\mathrm{\ell}k}$, the matrix is again symmetric.

If the functions ${\varphi}_{k}(x)$ are linearly independent on the set ${x}_{1},{x}_{2},\mathrm{\dots},{x}_{J}$, that is, the only solution of the system of equations

3.11.36 | $$\sum _{k=0}^{n}{c}_{k}{\varphi}_{k}({x}_{j})=0,$$ | ||

$j=1,2,\mathrm{\dots},J$, | |||

is ${c}_{0}={c}_{1}=\mathrm{\cdots}={c}_{n}=0$, then the approximation ${\mathrm{\Phi}}_{n}(x)$ is determined uniquely.

Now suppose that ${X}_{k\mathrm{\ell}}=0$ when $k\ne \mathrm{\ell}$, that is, the functions
${\varphi}_{k}(x)$ *are orthogonal with respect to weighted summation on the
discrete set* ${x}_{1},{x}_{2},\mathrm{\dots},{x}_{J}$. Then the system (3.11.33)
is diagonal and hence well-conditioned.

A set of functions ${\varphi}_{0}(x),{\varphi}_{1}(x),\mathrm{\dots},{\varphi}_{n}(x)$ that is linearly
independent on the set ${x}_{1},{x}_{2},\mathrm{\dots},{x}_{J}$ (compare (3.11.36))
can always be orthogonalized in the sense given in the preceding paragraph by
the *Gram–Schmidt* procedure; see Gautschi (1997a).

We take $n$ complex exponentials ${\varphi}_{k}(x)={\mathrm{e}}^{\mathrm{i}kx}$, $k=0,1,\mathrm{\dots},n-1$, and approximate $f(x)$ by the linear combination (3.11.31). The functions ${\varphi}_{k}(x)$ are orthogonal on the set ${x}_{0},{x}_{1},\mathrm{\dots},{x}_{n-1}$, ${x}_{j}=2\pi j/n$, with respect to the weight function $w(x)=1$, in the sense that

3.11.37 | $$\sum _{j=0}^{n-1}{\varphi}_{k}({x}_{j})\overline{{\varphi}_{\mathrm{\ell}}({x}_{j})}=n{\delta}_{k,\mathrm{\ell}},$$ | ||

$k,\mathrm{\ell}=0,1,\mathrm{\dots},n-1$, | |||

${\delta}_{k,\mathrm{\ell}}$ being Kronecker’s symbol and the bar denoting complex conjugate. In consequence we can solve the system

3.11.38 | $${f}_{j}=\sum _{k=0}^{n-1}{a}_{k}{\varphi}_{k}({x}_{j}),$$ | ||

$j=0,1,\mathrm{\dots},n-1$, | |||

and obtain

3.11.39 | $${a}_{k}=\frac{1}{n}\sum _{j=0}^{n-1}{f}_{j}\overline{{\varphi}_{k}({x}_{j})},$$ | ||

$k=0,1,\mathrm{\dots},n-1$. | |||

With this choice of ${a}_{k}$ and ${f}_{j}=f({x}_{j})$, the corresponding sum (3.11.32) vanishes.

The pair of vectors $\{\mathbf{f},\mathbf{a}\}$

3.11.40 | $\mathbf{f}$ | $={[{f}_{0},{f}_{1},\mathrm{\dots},{f}_{n-1}]}^{\mathrm{T}},$ | ||

$\mathbf{a}$ | $={[{a}_{0},{a}_{1},\mathrm{\dots},{a}_{n-1}]}^{\mathrm{T}},$ | |||

is called a *discrete Fourier transform pair*.

The direct computation of the discrete Fourier transform (3.11.38), that is, of

3.11.41 | ${f}_{j}$ | $={\displaystyle \sum _{k=0}^{n-1}}{a}_{k}{\omega}_{n}^{jk},$ | ||

${\omega}_{n}$ | $={\mathrm{e}}^{2\pi \mathrm{i}/n}$, | |||

$j=0,1,\mathrm{\dots},n-1$, | ||||

requires approximately ${n}^{2}$ multiplications. The method of the *fast
Fourier transform* (FFT) exploits the structure of the matrix $\mathbf{\Omega}$
with elements ${\omega}_{n}^{jk}$, $j,k=0,1,\mathrm{\dots},n-1$. If $n={2}^{m}$, then
$\mathbf{\Omega}$ can be factored into $m$ matrices, the rows of which contain
only a few nonzero entries and the nonzero entries are equal apart from
signs. In consequence of this structure the number of operations can be reduced
to $nm=n{\mathrm{log}}_{2}n$ operations.

Splines are defined piecewise and usually by low-degree polynomials. Given
$n+1$ distinct points ${x}_{k}$ in the real interval $[a,b]$, with
($a=$)$$($=b$), on each subinterval
$[{x}_{k},{x}_{k+1}]$, $k=0,1,\mathrm{\dots},n-1$, a low-degree polynomial is
defined with coefficients determined by, for example, values ${f}_{k}$ and
${f}_{k}^{\prime}$ of a function $f$ and its derivative at the nodes
${x}_{k}$ and ${x}_{k+1}$. The set of all the polynomials defines a function,
the *spline*, on $[a,b]$. By taking more derivatives into account,
the smoothness of the spline will increase.

For splines based on Bernoulli and Euler polynomials, see §24.17(ii).

For many applications a spline function is a more adaptable approximating tool than the Lagrange interpolation polynomial involving a comparable number of parameters; see §3.3(i), where a single polynomial is used for interpolating $f(x)$ on the complete interval $[a,b]$. Multivariate functions can also be approximated in terms of multivariate polynomial splines. See de Boor (2001), Chui (1988), and Schumaker (1981) for further information.

In computer graphics a special type of spline is used which produces a
*Bézier curve*.
A cubic Bézier curve is defined by four points. Two are endpoints:
$({x}_{0},{y}_{0})$ and $({x}_{3},{y}_{3})$; the other points $({x}_{1},{y}_{1})$ and $({x}_{2},{y}_{2})$
are control points. The slope of the curve at $({x}_{0},{y}_{0})$ is tangent to the
line between $({x}_{0},{y}_{0})$ and $({x}_{1},{y}_{1})$; similarly the slope at
$({x}_{3},{y}_{3})$ is tangent to the line between ${x}_{2},{y}_{2}$ and ${x}_{3},{y}_{3}$.
The curve is described by $x(t)$ and $y(t)$, which are cubic polynomials with
$t\in [0,1]$. A complete spline results by composing several Bézier curves.
A special applications area of Bézier curves is mathematical typography
and the design of type fonts. See Knuth (1986, pp. 116-136).