AMTH247 Lecture 11

Interpolation I

Reading: Heath §7.1, §7.2, §7.3.1

Contents

• Interpolation
  • The General Problem
  • Linear Interpolation
  • Reduction to Linear Equations
  • Basis Functions
  
  
• Polynomial Interpolation
  • Monomial Basis
  
  
• About this document ...

All the Scilab functions used in this lecture are in the file l11.sci in the directory for the lecture. The same thing will apply to future lectures.

Interpolation

The General Problem

The general problem of interpolation is the following:

Given data

$$(t_i, y_i), \quad i = 1, \dots n$$ (1)

Find a function $f(t)$ such that

$$y_i = f(t_i) , \quad i = 1, \dots n$$ (2)

The equation (2) is called the interpolation equation or interpolation condition. It says that the function $f(t)$ passes through the data points. A function $f(t)$ satisfying the interpolation condition is called an interpolating function for the data. There are numerous applications calling for the construction of an interpolating function for a set of data. Some of these are:

1. Plotting a smooth curve through discrete data points.
2. Reading between lines of a table.
3. Numerical manipulation, e.g. integration, of tabular data.

Interpolation, as we will see later, is the basis of variety of other numerical algorithms. Heath, §7.1, gives extended account of the uses of interpolation.

Linear Interpolation

Of course, for any set of data there is an unlimited number of interpolating functions. We will restrict attention to linear interpolation where the interpolating function can be expressed as a linear combination

$$f(t) = \sum_{j=1}^n x_j \phi(t) = x_1 \phi_1(t) + \dots + x_n \phi_n(t)$$ (3)

of $n$ basis functions $\phi_j(t)$. We will see shortly, that for the linear interpolation problem to be well posed, the number of basis functions must be equal to the number of data points. The problem of determining the interpolating function then becomes one of determining the coefficients $x_j$ of the basis functions $\phi_j(t)$. These coefficients are determined by the interpolation conditions (2).

Example:

We will interpolate the data

$t$	-1	0	1
$y$	1.4	0.8	1.7

by a function of the form

$$f(t) = x_1 e^{-t} + x_2 + x_3 e^{t}$$

In this case the basis functions are

$$\phi_1(t) = e^{-t}, \quad \phi_2(t) = 1, \quad \phi_3(t) = e^{t} .$$

Our problem is to determine the coefficients $x_1$, $x_2$, $x_3$. Applying the interpolation conditions, equation (2), to the function $f(t)$ at each of the three data points gives the system of equations

$$1.4 = x_1 e^{1} + x_2 + x_3 e^{-1}$$

$$0.8 = x_1 e^{0} + x_2 + x_3 e^{0}$$

$$1.7 = x_1 e^{-1} + x_2 + x_3 e^{1}$$

This is a set of three linear equations for the three unknown coefficients, $x_1$, $x_2$, $x_3$. It is easily solved in Scilab

-->a = [exp(1) 1 exp(-1)
-->1 1 1
-->exp(-1) 1 exp(1)]
 a  =
 
!   2.7182818    1.    0.3678794 !
!   1.           1.    1.        !
!   0.3678794    1.    2.7182818 !
 
-->b = [1.4 0.8 1.7]'
 b  =
 
!   1.4 !
!   0.8 !
!   1.7 !
 
-->x = a\b
 x  =
 
!   0.6266863 !
! - 0.5810104 !
!   0.7543241 !

So our interpolating function is

$$f(t) = 0.627 e^{-t} - 0.581 + 0.754 e^{t}$$

It is easy to graph the interpolating function:

-->t = (-1:0.01:1)';
 
-->y = x(1)*exp(-t) + x(2) + x(3)*exp(t);
 
-->plot2d(t,y)

Reduction to Linear Equations

Just as in the example above, the determination of the coefficients in a general linear interpolation problem reduces to solving a system of linear equations. Given data

$$(t_i, y_i), \quad i = 1, \dots n ,$$

and basis functions $\phi_1(t), \dots, \phi_n(t)$, our problem is to determine coefficients $x_j$ such that the function

$$f(t) = x_1 \phi_1(t) + \dots + x_n \phi_n(t)$$

satisfies the interpolation conditions

$$y_i = f(t_i) , \quad i = 1, \dots n .$$

In terms of the basis functions, the interpolation conditions are

$$y_i = f(t_i) = x_1 \phi_1(t_i) + \dots + x_n \phi_n(t_i) , \quad i = 1, \dots n ,$$

or, writing the equation for each data point,

$$y_1 = x_1 \phi_1(t_1) + \dots + x_n \phi_n(t_1)$$

$$y_2 = x_1 \phi_1(t_2) + \dots + x_n \phi_n(t_2)$$

$$\vdots \vdots$$

$$y_n = x_1 \phi_1(t_n) + \dots + x_n \phi_n(t_n) .$$

This is a system of $n$ linear coefficients for the $n$ coefficients $x_1, \dots , x_n$. In matrix form, these equations are

$$\begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} = \b... ...n) \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}$$ (4)

The $n \times n$ matrix in this equation is formed by evaluating the $n$ basis functions $\phi_j(t)$ at the $n$ data points $t_i$. It is called the basis matrix.

Example:

Again we will interpolate the data

$t$	-1	0	1
$y$	1.4	0.8	1.7

this time by a quadratic polynomial

$$f(t) = x_1 + x_2 t + x_3 t^2 .$$

In this case the basis functions are

$$\phi_1(t) = 1, \quad \phi_2(t) = t, \quad \phi_3(t) = t^2.$$

The basis matrix is

$$\begin{bmatrix} 1 & t_1 & t_1^2 \\ 1 & t_2 & t_2^2 \\ 1 & t_... ...x} = \begin{bmatrix} 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \end{bmatrix}$$

and we need to solve

$$\begin{bmatrix} 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \end{bm... ... \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix}1.4 \\ 0.8 \\ 1.7 \end{bmatrix}$$

Solving the problem in Scilab:

-->a = [1 -1 1
-->1 0 0
-->1 1 1]
 a  =
 
!   1.  - 1.    1. !
!   1.    0.    0. !
!   1.    1.    1. !
 
-->y = [1.4 0.8 1.7]'
 y  =
 
!   1.4 !
!   0.8 !
!   1.7 !
 
-->x = a\b
 x  =
 
!   0.8  !
!   0.15 !
!   0.75 !

giving the interpolating polynomial

$$f(t) = 0.8 + 0.15 t + 0.75 t^2 .$$

We can evaluate and plot this polynomial

-->t = (-1:0.01:1)';                     
 
-->y = x(1) + x(2)*t + x(3)*t.^2;
 
-->plot2d(t,y)

The two different interpolating functions we have constructed for our data give very similar, but not identical, graphs.

Basis Functions

Choosing different sets of basis functions will, for a given set of data, produce interpolating functions from different classes of functions. Some common classes of functions are

1. Polynomials
2. Trigonometric functions
3. Splines

The choice of class of interpolating functions to use for any particular set of data, depends on the data; we choose a class of functions which seems to give a good representation of the data. Any such choice must be somewhat subjective, but sometimes we have theoretical reasons for choosing a particular class of interpolating functions. Even for a single class of functions there can be different choices of basis functions. We will give several sets of basis functions for polynomials later. Different bases, in this case, do not give different interpolating functions, just different representations of that function. The choice of basis has relevance to the computational problem

1. Determining the coefficients $x_j$ from the linear equations (4):
   the choice of basis functions determines the structure of the basis matrix.
   1. For some choice of basis functions the basis matrix may take a simple form, e.g. a triangular matrix, for which determining the coefficients is easy.
   2. More generally, we would like the basis matrix to be well-conditioned. This is particularly important if the coefficients $x_j$ are of interest in themselves and need to be determined accurately.
   3. Even when the basis matrix is ill-conditioned, the solution for the coefficients $x_j$ may still give a small residual, although it is quite inaccurate. A small residual means that the interpolating function still does a good job of interpolating the data.
2. Evaluating the interpolating function:
   for most purposes the interpolating function will be evaluated many times. Therefore it is desirable to choose basis functions which can be evaluated easily.

Polynomial Interpolation

Given any set of data

$$(t_i, y_i), \quad i = 1, \dots n ,$$

there is a unique polynomial of degree at most $n-1$ which interpolates the data. Note that a polynomial of degree $n-1$ has $n$ coefficients, the same as the number of data points.

Monomial Basis

The most obvious basis for the polynomials of degree at most $n-1$ is the monomial basis of powers of $t$

$$\phi_1(t) = 1, \quad \phi_2(t) = t, \quad \phi_3(t) = t^2, \quad \dots, \quad \phi_n(t) = t^{n-1}$$

which gives the familiar representation of the interpolating polynomial

$$p(t) = x_1 + x_2 t + x_3 t^2 + \dots + x_n t^{n-1} .$$

The Basis Matrix

For data $t_1, \dots, t_n$ the basis matrix is

$$\mathbf{A} = \begin{bmatrix} 1 & t_1 & t_1^2 & \dots & t_1^{n-1}... ...ots & \ddots & \vdots \\ 1 & t_n & t_n^2 & \dots & t_n^{n-1} \end{bmatrix}$$

Note that the number of data points determines the number of basis functions. Below is a Scilab function to compute the basis matrix for the monomial basis at a given set of data values. The key point is to note that the columns of the basis matrix are just powers of the vector of data points considered as a column vector. In the function monbasis (t) the vector t of data points is assumed to be column vector.

function a = monbasis (t)
  n = length(t)
  a = zeros(n,n)
  for j = 1:n
    a(:,j) = t.^{j-1}
  end
endfunction

Now, to compute the coefficients of the interpolating polynomial we need to solve the linear system

$$\mathbf{A} \mathbf{x} = \mathbf{y}$$

where $\mathbf{A}$ is the basis matrix given above, $\mathbf{x}$ is the column vector of coefficients $x_j$, and $\mathbf{y}$ is the column vector of data values $y_j$. We can write a Scilab function to compute the coefficients of the interpolating polynomial in the monomial bases which interpolates a given set of $(t,y)$ data. It is very simple:

function x = moninterp (t, y)
  a = monbasis(t)
  x = a\y
endfunction

Example:

We will do the same polynomial interpolation as we did earlier, this time using moninterp:

-->t = [-1 0 1]'
 t  =
 
! - 1. !
!   0. !
!   1. !
 
-->y = [1.4 0.8 1.7]'
 y  =
 
!   1.4 !
!   0.8 !
!   1.7 !
 
-->x = moninterp(t,y)
 x  =
 
!   0.8  !
!   0.15 !
!   0.75 !

Horner's Rule

Once we have the coefficients of the interpolating polynomial we need a procedure to evaluate it. There is a clever method for evaluating polynomials. The standard method writes a polynomial

$$p(t) = x_1 + x_2 t + x_3 t^2 + \dots + x_n t^{n-1} ,$$

evaluates the powers, multiplies by the coefficients and adds. If we evaluate $p(t)$ by writing it

$$p(t) = x_1 + t(x_2 + t(x_3 + t(\dots (x_{n-1} + x_n t) \dots )))$$

we get Horner's rule for evaluating a polynomial. Horner's rule requires less arithmetic and is less sensitive to rounding errors than the standard method. It is easy modify Horner's method to evaluate the derivatives of $p(t)$ at the same time the polynomial is evaluated. Here it is in Scilab function moneval which takes a vector x of coefficients of a polynomial and a vector t of points at which to evaluate the polynomial. It returns the vector y of the values of the polynomial at the points t evaluated by Horner's rule.

function y = moneval(x, t)
  n = length(x);
  y = zeros(t);
  for i = n:-1:1
    y = y.*t + x(i);
  end
endfunction

An Example

We will use the data from Heath, p 309.

$t$	1.0	2.0	3.0	4.0	5.0	6.0
$y$	1.9	2.7	4.8	5.3	7.1	9.4

First we enter and plot the data

-->t =(1:6)'
 t  =
 
!   1. !
!   2. !
!   3. !
!   4. !
!   5. !
!   6. !
 
-->y = [1.9 2.7 4.8 5.3 7.1 9.4]'
 y  =
 
!   1.9 !
!   2.7 !
!   4.8 !
!   5.3 !
!   7.1 !
!   9.4 !

-->plot2d(t,y, style = -1)

Now we will compute the interpolating polynomial of degree 5 for the data

-->x = moninterp(t, y)
 x  =
 
!   20.6      !
! - 40.241667 !
!   29.820833 !
! - 9.6291667 !
!    0.0791667 !

and evaluate and plot the interpolating polynomial

-->tt = (0:.01:7)';       
 
-->yy = moneval(x, tt);   
 
-->plot2d(tt, yy)

Finally we can compute the residual, the difference between the data values and the polynomial values at the data points.

-->y1 = moneval(x, t)
 y1  =
 
!   1.9 !
!   2.7 !
!   4.8 !
!   5.3 !
!   7.1 !
!   9.4 !
 
-->res = y1-y        
 res  =
 
   1.0E-13 *
 
! - 0.0133227 !
!   0.0976996 !
! - 0.0976996 !
! - 0.0976996 !
!   0.5861978 !
!   0.2664535 !

These residuals are rather large, suggesting that the solution may be inaccurate. We can test this by looking at the condition number of the basis matrix:

-->a = monbasis(t)
 a  =
 
!   1.    1.    1.     1.      1.       1.    !
!   1.    2.    4.     8.      16.      32.   !
!   1.    3.    9.     27.     81.      243.  !
!   1.    4.    16.    64.     256.     1024. !
!   1.    5.    25.    125.    625.     3125. !
!   1.    6.    36.    216.    1296.    7776. !
 
-->cond(a)        
 ans  =
 
    731200.94  
 
-->cond(a)*%eps   
 ans  =
 
    1.624E-10

and we see that while the residual is large, the expected error in the coefficients is larger still. This problem of the basis matrix having a large condition number is typical of monomial bases.

About this document ...

AMTH247 Lecture 11

Interpolation I

This document was generated using the LaTeX2HTML translator Version 2K.1beta (1.61)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split +0 lecture.tex

The translation was initiated by amth247 on 2003-05-05