TUE Scientific Computing Group, 2N330 Numerical Methods and Algorithms Numerieke Methoden en Algoritmiek

Course notes by Dr. Joseph M. Maubach.

The Interpolation of Functions

Contents

• 1 Introduction
• 2 The Interpolant, using a Monomial Basis
• 3 The Interpolant, using a Newton basis
• 4 The Interpolant, using a Lagrangian basis
• 5 The Interpolant, Overview
• 6 The Interpolation Error
• 7 Different Types of Interpolation
• 8 Efficient Polynomial Evaluation
• 9 Exercises
• Index

See also the textbook by Kreyszig, 17.3.

1 Introduction

In all of the interpolation material, at times $t_0 < t_1 < \ldots < t_n\in[a, b]$ function values ${\bf f}= [f(t_0), \ldots, f(t_n)]$. are sampled (measured). Hence, in all sections, we use $n+1$ sample times and start with index 0. The reason is that a polynomial of degree $n$ (example: degree 2: $p(t) = c_0 + c_1 t + c_2 t^2$) has $n+1$ associated degrees of freedom, of which the first one (constant) has index 0.

This section discusses the problems described next. Problem: A function f is know in sampled form:

$i$	0	1	2	3	4
$t_i$	-2	-1	0	1	2
$f(t_i)$	18	0	0	0	6

Table with sampled data

Determine

$$f'(0)$$ (1)

and

$$\int_{-2}^{ 2} f(t) dt.$$ (2)

Problem: We do not know $f$. Thus, we can not determine the answers. But, we can approximate the answers.

Proposed solution: Find a degree 4 polynomial $p$ (5 sample points) which approximates (resembles) $f$. Because

$$f'(0) \approx p'(0),$$

approximate the integral with

$$\int_{-2}^{ 2} f(t) dt \approx \int_{-2}^{ 2} p(t) dt.$$

Definition 1   A function $g$ interpolates $f$ at points $t_i$ if $g(t_i ) = f(t_i)$ for all $i\geq0$.

Figure 1: A degree 8 interpolant of Runge's function $f(t) = 1 / (1 + 25t^2)$

2 The Interpolant, using a Monomial Basis

Definition 2   : The monomials $\{1, t, ...., t^n\}$ form a basis for all polynomials through degree $n$ (i.e., up to and including degree $n$).

Reconsider the problem of interpolating $f$

$i$	0	1	2	3	4
$t_i$	-2	-1	0	1	2
$f(t_i)$	18	0	0	0	6

with a polynomial of degree $4$. Write the polynomial to be determined as

$$p(t) = c_0 + c_1 t + ... + c_4 t^4$$

Determine $c_i$:

$$\begin{array}{rcrcrcrc} 18 &=& f(t_0) &=& p(t_0) &=& c_0 -2\c... ...cdot c_1 +4\cdot c_2 + 8\cdot c_3 + 16\cdot c_4; \end{array}$$ (3)

This yields the system of equations (matrix form):

$$\begin{bmatrix}1 &-2 & 4 & -8 &16 1 &-1 & 1 & -1 & 1 1 & 0 ... ... c_4 \end{bmatrix} = \begin{bmatrix}18 0 0 0 6 \end{bmatrix}$$ (4)

abbreviated by

$$A {\bf c}= {\bf b}.$$

Note that row $i+1$ of $A$ is of the form

$$[t_i^0, t_i^1, t_i^2, t_i^3, t_i^4],$$

and that entry ${\bf b}_{i+1} = f(t_i)$, for $i=0, \ldots, 4$. This system of equations can be solved with matlab:

 A = [1, -2, 4, -8, 16; 1, -1, 1, -1, 1; 1, 0, 0, 0, 0; 1, 1, 1, 1, 1; 1, 2, 4, 8, 16]
 b = [18; 0; 0; 0; 6]
 c = A \ b

This gives:

$$c = [0, 1, -1, -1, 1].$$

This means that the coefficients are $c_0 = 0$, etc. and that the polynomial approximation is

$$p(t) = t - t^2 - t^3 + t^4.$$ (5)

Solution to the original problems 1 and 2 are:

$$f'(t) \approx p'(t) = 1 - 2t - 3t^2 + 4t^3$$

which attains value 0 at $t = 1$. and

$$\int_{-2}^{ 2} f(t) dt \approx \int_{-2}^{ 2} p(t)dt = 112/15.$$

Instead of the separate steps, one can call a polynomial fit routine

polyfit([-2, -1, 0, 1, 2], [18, 0, 0, 0, 6], 4)

The returned answer is

    1.0000   -1.0000   -1.0000    1.0000    0.0000,

i.e., the coefficients are in the reverse order. The commands

c = polyfit([-2, -1, 0, 1, 2], [18, 0, 0, 0, 6], 4)
t = -2:0.01:2
f = polyval(c,t)
plot(t, f)

also plot the polynomial. The arguments of polyval are the computed sequence of coefficients, and a sequence of sample data $t = [-2, -1.99, \ldots, 1.99, 2]$.

Example 1   The roots of a polynomial with rational coefficients can be found with the maple command solve (in matlab):

 syms x
 solve(x^5 - 15*x^4 + 85*x^3 -225*x^2 + 274*x -120)

and be approximated with the command roots:

 c = [1, -15, 85, -225, 274, -120]
 roots(c)

The roots turn out to be $1, 2, 3, 4, 5$.

Theorem 1   Let $t_0, t_1, \ldots, t_n$ be the $n+1$ sample points, and let ${\bf f}= [f(t_0), \ldots, f(t_n)]$ be the measured values. The related matrix is a VanderMonde Matrix:

$$A = \begin{bmatrix}1 & t_0& t_0^2& \cdots & t_0^{n} 1 & t_1& t... ...{n} \vdots&& & & \vdots 1 & t_n& t_n^2& \cdots & t_n^{n} \end{bmatrix}$$ (6)

When all $t_i$ are different, one can show that

$$\det(A) = \prod_{{\begin{matrix}i, j = 0 i < j\end{matrix}}}^n (t_i - t_j) \neq 0,$$

which implies that the coefficients

$${\bf c}= A^{-1}{\bf f}$$

define the unique interpolant

$$p(t) = c_0 + c_1 t + c_2 t^2 + \ldots + c_n t^n.$$ (7)

Remark 1   Solving a system $A{\bf x}= {\bf b}$, with $A$ a VanderMonde matrix is a dangerous numerical undertaking. Computations with $A$ cause a lot of floating point round-off, because for $t\in(0, 1)$ and large $n$, $t^n$ approximates zero, whence $det(A) \approx 0$.

Remark 2   The interpolant based on $n+1$ interpolation points need not be of degree $n$: Assume we intend to interpolate the three data points $(-1, 1)$, $(0, 0)$ and $(1, -1)$, all of form $(t_i, f(t_i)$. Because there are three points, we use a polynomial of degree 2

$$c_0 + c_1t + c_2t^2.$$

We will find that $c_0 = 0$, $c_1 = -1$ and that $c_2 = 0$, i.e., the interpolant is of degree $1$.

Remark 3   In matlab, entries of the sequence

 x = [1, 2, 3, 4]

are indexed with

 x(1)
 x(2)
 x(3)
 x(4)

This notebook indexes $t_0, \ldots t_n$, starting from 0.

Example 2   Problem: The function

$$f(t) = 1 / (1 + 25t^2)$$ (8)

is sampled at points $t_i$, $i=0, \ldots, 4$. This leads to the data set:

$i$	0	1	2	3	4
$t_i$	-1	-1/2	0	1/2	1
$f(t_i)$	1/26	4/29	0	4/29	1/26

Determine (approximations of)

$$f'(0), f''(0)$$ (9)

as well as

$$\int_{-1}^{ 1} f(t) dt,$$ (10)

from the sampled data from equation (8).

In order to see how good polynomial interpolation works, because we know $f$ from equation 8, we compute the exact answers first. With matlab

 syms x
 int(1 / (1 + 25*x^2), x)

As an alternative, with mathematica

 F[x_] = 1 / (1 + 25x^2)
 Integrate[F[x], {x, -1, 1}] = 2/5 *ATan[5]

which is approximately $0.54936$. In a similar manner, we find

$$\begin{array}{lcr} f'(t) &=& \frac{-50 t }{ (1 + 25t^2)^2} \\... ...}{ (1 + 25t^2)^3} - \frac{- 50 }{ (1 + 25t^2)^2} \end{array}$$ (11)

whence

$$\begin{array}{lcr} f'(0) &=& 0 f''(0) &=& -50 \end{array}$$ (12)

Indeed, at $t=0$, $f$ has a horizontal tangent line, and it is concave. These answers will be used to verify the accuracy of the interpolants.

The interpolant based on all 5 sample points is unique. It doesn't matter whether it has been calculated with the use of the monomial, Lagrangian or Newton's basis: In matlab:

c = polyfit([-1, -1/2, 0, 1/2, 1], [1/26, 4/29, 1, 4/29, 1/26], 4)
x = -1:0.01:1
y = polyval(c,x)
plot(x, y)

In mathematica:

 Expand[InterpolatingPolynomial[{{-1, 1/26}, {-1/2, 4/29}, {0, 1}, {1/2, 4/29}, {1, 1/26}}, x]]

is:

$$p_4(t) = 1 - \frac{3225}{754} t^2 + \frac{1250}{377} t^4.$$ (13)

The interpolant based on 3 points $t = -1/2, 0, 1/2$ is the second order polynomial:

$$p_2(t) = 1 - \frac{100}{29} t^2$$ (14)

We find

$$\begin{array}{lclcl} \int_{-1}^{ 1} p_4(t) dt &=& \frac{179}{... ...p_2(t) dt &=& \frac{-26}{87} &\approx& -0.298851 \end{array}$$ (15)

with

 Integrate[p_4[x], {x, -1, 1}]
 Integrate[p_2[x], {x, -1, 1}]

Note: The mathematica command

 N[Solve[D[p[x], x] == 0, x]])

can be used to compute the minima of $p_4$.

3 The Interpolant, using a Newton basis

Instead of the monomial and Lagrange basis we can construct our interpolant with the use of a Newton basis.

Example 3   Look at the original example 1. Its interpolant was computed in (5). Here, we take a different basis

$$\begin{array}{lclcl} \displaystyle n_0(t) &=& 1 & & \\ \displ... ... - t_1)(t - t_2)(t - t_3) &=& (t+2)(t+1)t(t - 1) \\ \end{array}$$ (16)

That is, our interpolant is of the form:

$$p(t) = c_0 + c_1(t - t_0) + c_2(t - t_0)(t - t_1) + \cdots + c_4(t - t_0)(t - t_1)(t - t_2)(t - t_3),$$ (17)

i.e., of the form

$$p(t) = c_0 + c_1(t + 2) + c_2 (t+2)(t+1) + c_3 (t+2)(t+1)t + c_4 (t+2)(t+1)t(t - 1).$$ (18)

Now, we demand $p = f$ at all $t = t_i$:

$$\begin{array}{lclclclcl} \displaystyle 18 &=& f(t_0) &=& f(-2... ... &=& p(2) &=& c_0 + 4c_1 + 12c_2 + 24c_3 + 24c_4 \end{array}$$ (19)

Typical for a Newton basis, the resulting system can be solved with a forward substitution. We find

$$c = [18, -18, 9, \ldots ].$$ (20)

Note that $c_0 = 18$, and not 0 as for the monomial basis, though the interpolant - write it out - still is $p(t) = t - t^2 - t^3 + t^4$.

4 The Interpolant, using a Lagrangian basis

The computation of the interpolant with the use of a monomial basis costs a lot of computations, which must be redone for different function $f$. With the use of a special basis this can be avoided.

Definition 3   Let $t_0 < t_1 < \ldots < t_n$ be the $n+1$ sample points. Then the related Lagrange polynomial $l_j(t)$, $j=0,\ldots, n$ is defined by:

$$\begin{array}{lcr} \displaystyle l_j(t) &=& \prod_{k = 0, k \... ... a_j &=& \prod_{k = 0, k \neq j}^{n} (t_j - t_k) \end{array}$$ (21)

for all $j=0,\ldots, n$. These polynomials are of degree $n$, and form a basis for for all polynomials of degree $n$ because of:

Theorem 2  

$$\begin{array}{lcr} \displaystyle l_j(t_i) &=& 1 \quad \text{i... ...= j \displaystyle &=& 0 \quad \text{elsewise} \end{array}$$ (22)

for all $i,j=0,\ldots,n$.

Example 4   Let $[t_0, t_1, t_2] = [-1, 0, 1]$. Then

$$\begin{array}{lcr} l_0(t) &=& \frac{(t - 0)(t - 1)}{(-1 - 0)(... ...2(t) &=& \frac{(t - -1)(t - 0)}{(1 - -1)(1 - 0)}. \end{array}$$ (23)

The zeros of $l_j(t)$ follow from its factored from, and the fact that $l_j(t_j) = 1$ is due to the fact that $l_j(t)$'s enumerator and denominator are identical for $t_j$.

Remark 4   Because all Lagrange basis functions are of degree $n$, we should write $l^n_j$, instead of $l_j$.

Theorem 3   The polynomial interpolant $p(t)$ of $f$ based on $n+1$ interpolation points $t_0, \ldots t_n$, formulated with the use of the Lagrange basis, is

$$p(t) = \sum_{j = 0}^n f(t_j) l^n_j(t).$$ (24)

Remark 5   Of course, the interpolant $p(t)$ is identical to the interpolant in equation 5 found in section 2. However, no system of equations has solved.

Proof: Due to the Kronecker property, for all $i=0,\ldots,n$:

$$p(t_i) = \sum_{j = 0}^n f(t_j) l^n_j(t_i) = f(t_i).$$

This is a results of property (22).

Example 5   Look at example 1. For $5$ sample points (degree $n = 4$), consider the first Lagrange basis polynomial:

$$\begin{array}{lclcr} \displaystyle a_1 &=& (-2 - -1)(-2 - 0)(... ...\displaystyle &=& t(t + 1)(t - 1)(t - 2)/24. & & \end{array}$$ (25)

Its properties are verified:

$$\begin{array}{lclcr} \displaystyle l_0(t) &=& 0 \quad\text{fo... ... l_0(-2) &=& -1\cdot -2\cdot -3\cdot -4/24 &=& 1 \end{array}$$ (26)

Moreover (factors reordered):

$$\begin{array}{rcl} \displaystyle l_1(t) &=& (2 - t) (-1 + t) ... ...style l_4(t) &=& (-1 + t) t (1 + t) (2 + t) / 24 \end{array}$$ (27)

The interpolant $p$ of $f$ - computed using the Lagrangian basis - is:

$$\begin{array}{rcl} \displaystyle p &=& f(t_0)l_0(t) + \ldots ... ...+ t)^2 \displaystyle &=& t - t^2 - t^3 + t^4. \end{array}$$ (28)

Remark 6   The use of a Lagrange basis for the construction of the interpolant has its advantage: No system of equations has to be solved. However, when we need to evaluate or differentiate or integrate the interpolant, we have to evaluate/differentiate/integrate all Lagrangian basis functions.

5 The Interpolant, Overview

We have seen that the different manners (different bases) for the computation of the interpolant have different characteristics. Here is an overview:

Type of basis	monomial	Newton	Lagrange
Linear system of equations	full - bad	triangular	none
Basis functions	simple	less simple	complex
Polynomial evaluation	cheap	less cheap	expensive

Here, bad stands for badly conditioned. However, polynomial interpolants based on equi-distant grid-points can perform in a poor manner, as the next example shows.

Example 6   Let $f$ be the infamous function of Runge:

$$f(t) = 1 / (1 + 25t^2),$$

defined in $[-1, 1]$. This function is smooth, but its higher derivatives grow larger and larger in magnitude. A plot of $f$ is obtained with:

 x = -1:0.01:1;
 y = 1 ./ (1 + 25*x.^2);
 plot(x, y)

Figure 2: A degree 8 and 16 interpolant of Runge's function $f(t) = 1 / (1 + 25t^2)$

Figure 3: A degree 32 interpolant of Runge's function $f(t) = 1 / (1 + 25t^2)$

Section 7 presents different interpolation techniques (non-equidistributed points $t_i$, non-polynomial interpolants, etc.) which have better approximation properties.

6 The Interpolation Error

Theorem 4   Let $f\colon[a, b]\to{\mathbb{R}}$ be a smooth function. Let $t_0, \ldots, t_n\in[a, b]$ be the interpolation points, and $p$ be the lowest degree polynomial which interpolates $f$ at these points. Then the interpolation error $R(f)$ is of the form

$$f(t) - p(t) = R(f) = \left(\prod_{i=0}^n (t - t_i)\right) \cdot \frac{f^{(n+1)}(\theta)}{(n+1)!}.$$ (29)

The proof is based on a repeated application of Rolle's theorem.

7 Different Types of Interpolation

The interpolant can be computed with the use of bases different from the ones we have seen. To mention a few widely used ones:

• Hermite basis, which makes use of matching higher derivatives at $t_i$;
• Spline basis, which makes use of matching higher derivatives at $t_i$;
• Bezier basis, which makes use of so called control points in between $t_i$;
• Tchebyshev basis, which makes use of non-equidistant sample points $t_i$;

8 Efficient Polynomial Evaluation

Let $k = 4$. The computation of $t^i$ costs $i-1$ operations, the computation of $c_i t^i$ costs $i$ operations and the evaluation of

$$p(t) = c_0 + c_1 t + ... + c_4 t^4$$ (30)

costs $k$ additions and

$$1 + \cdots + k = \frac{{\scriptstyle 1}}{{\scriptstyle 2}}k(1+k).$$ (31)

multiplications. The total amount of operations is

$$k + \frac{{\scriptstyle 1}}{{\scriptstyle 2}}k(1+k) = \frac{{\scriptstyle 1}}{{\scriptstyle 2}}k^2 + \frac32k,$$ (32)

which is quadratic in the polynomial degree $k$.

But it is possible to eliminate common subexpressions, or factor out common subexpressions. For instance, after $t^2$ is computed with the use of one multiplication, $t^3$ can be computed with one additional multiplication:

$$t^3 = t^2 \cdot t.$$

This technique is called Horner's method: Use nested evaluation of $p$:

$$\begin{array}{rcl} p(t) &=& (((c_4t + c_3)t + c_2)t + c_1)t +... ... c_0\\ &=& c_4t^4 + c_3t^3 + c_2t^2 + c_1t + c_0.\\ \end{array}$$ (33)

In this manner, the evaluation costs $2k$ operations: $k$ multiplications (with $t$) and $k$ additions. Hence the total amount of operations is linear with respect to the degree $k$.

Example 7   An example of a polynomial $p$ and its nested form:

$$\begin{array}{rcl} p(t) &=& 1 - t + 2t^2 + 4t^3 &=& ((4t +2)t - 1)t + 1. \end{array}$$ (34)

Example 8   Another example of a polynomial $p$ and its nested form:

$$\begin{array}{rcl} p(t) &=& 1 - t^3 + t^5 &=& (t^2 + 1)t^3 + 1. \end{array}$$ (35)

The nested form is computed just as above, assuming the coefficients in front of $t$, $t^2$ and $t^4$ are zero.

9 Exercises

EXAM 1   Answer true or false.

a)

Using a polynomial interpolant of high degree based on equi-distant interpolation points in $[a, b]$ a smooth function is always approximated well near the limit points $a$ and $b$.

b)

An evaluation of $p(t) = 1 + t^3 + t^5$ can be performed with eight (8) floating point operations.

c)

All polynomials of degree $n$ can be evaluated with minimally $n-1$ multiplications and $n$ additions.

d)

Let $p$ be a polynomial of degree $n$ and $t$ be a real number. The fastest possible evaluation $p(t)$ costs $\frac{{\scriptstyle 1}}{{\scriptstyle 2}}n (n + 1)$ floating point operations.

EXERCISE 1   Make a file called reacdata.prn which contains:

    0.57    0.91    0.12
    1.01    0.79    0.16
    1.57    0.71    0.25
    2.02    0.69    0.38
    2.48    0.60    0.41
    3.05    0.51    0.52
    4.01    0.41    0.59
    4.98    0.34    0.62
    5.50    0.31    0.65
    6.02    0.29    0.69
    7.04    0.25    0.71
    7.49    0.23    0.76
    8.03    0.19    0.78
    8.49    0.16    0.77

This file contains two concentrations, measured in time. The first column is the time (seconds), the other two columns are the concentrations. Visualize the measured concentrations in matlab:

  load reacdata.prn
   t = reacdata(:,1);
  c1 = reacdata(:,2);
  c2 = reacdata(:,3);
  plot(t,c1,'o',t,c2,'x')

Interpolate the first concentration data, first with a polynomial of degree $4$:

c = polyfit(t, c1, 4)
x = 0:0.01:9
y = polyval(c,x)
plot(x, y, '-', t, c1, 'o')

Next, interpolate it with a polynomial of degree 13 (fourteen samples):

c = polyfit(t, c1, 13)
x = 0:0.01:9
y = polyval(c,x)
plot(x, y, '-', t, c1, 'o')

Also interpolate the second concentration with a polynomial of degree $4$ and degree $13$.

EXERCISE 2   Interpolate the following data with a degree 4 polynomial:

$i$	0	1	2	3	4
$t_i$	-2	-1	0	1	2
$f(t_i)$	4	1	0	1	4

Is the (unique) interpolant of degree 4?

Table with sampled data

EXAM 2   Answer true or false.

a)

The polynomial interpolant through $5$ samples $(t_i, f(t_i))$ must be of degree 4.

b)

The polynomial interpolant computed with Newton's and Lagrange's methods must be identical.

EXAM 3   Van een functie $f$ is de onderstaande tabel met functiewaarden gegeven:

$i$	0	1	2
$t_i$	1	-1	0
$f(t_i)$	0	0	1

a)

Bereken de interpolant $p$ van $f$ met behulp van de monomiale basis $1, t$ en $t^2$.

b)

Geef de drie Lagrange basis polynomen van de graad 2 waarmee de interpolant $p$ van $f$ kan worden opgebouwd en bereken deze interpolant met behulp van Lagrange's methode.

c)

Geef de drie Newton basis polynomen van de graad 2 waarmee de interpolant $p$ van $f$ kan worden opgebouwd en bereken deze interpolant met behulp van Newton's methode.

d)

Geef een formule voor de interpolatie fout

$$f(t) - p(t).$$

EXERCISE 3   Determine the coefficients $c_3$ and $c_4$ in (20). Insert $c_0, \ldots, c_4$ into the polynomial $p$ in (18). Show that indeed

$$p(t) = t - t^2 - t^3 + t^4.$$

Index

Bezier

7

eliminate common subexpressions

8

equi-distant grid-points

5

factor out common subexpressions

8

function of Runge

5

Hermite

7

Horner's method

8

interpolation error

6

Lagrange

4

monomials

2

Spline

7

Tchebyshev

7

VanderMonde

2

Back to Top.