Inverse Interpolation

The inverse interpolation method can be used to found a zero of a function f (x) by
computing f 1 (0) , where x  f 1 ( y ) is a polynomial that interpolates a table of
( f 1 ( x), x) values.

For example, let's suppose that we want to calculate a zero of the function
f ( x)  x 3  15x  4 knowing that this is close to x  0.3 .
We then first evaluate the function in three points close to x  0.3 :

x 0.2 0.3 0.4

y = f(x) 1.008 -0.473 -1.936

Second step: inverse polynomial interpolation

The following points are used to generate a quadratic interpolating polynomial

y 1.008 -0.473 -1.936

x 0.2 0.3 0.4

The polynomial can be generated in a number of ways including: Newton’s Divided

Difference Interpolating Polynomials .

xl f(xl) First Second
1.008 0.2 -0.0675 0.00028963
-0.473 0.3 -0.0684
-1.936 0.4

Thus the interpolating polynomial is:

f2(y) = 0.2+(y-1.008)( -0.0675)+(y-1.008)(y+0.473)( 0.00028963)

So an approximation of the zero of the function is:

f2(0) = 0.2+(0-1.008)( -0.0675)+(0-1.008)(0+0.473)( 0.00028963)

f2(0) = 0.2679.

Given x0 , f ( x0 ) , x1 , f ( x1 ) ,......, xn 1 , f ( xn 1 ) , xn , f ( xn )  and a function value f ( xk ) , how do
you determine x k ?
 1. Using one of the interpolation methods (Newton or Lagrange), obtain the function
f(x) :
f ( x)  a0  a1 x  a2 x 2    an x n
2. Plug in the value of f ( xk ) in the above equation and find one of the real roots of the
resultant equation using a root finding method:
a0  a1 x  a2 x 2    an x n  f ( xk )  0
Example:
Employ inverse interpolation using a cubic interpolating polynomial and root finding to
determine the value of x that corresponds to f(x) = 0.23 for the following data:

xi f(xi)
3 0.3333
4 0.25
5 0.2
6 0.1667

Solution:
First step: polynomial interpolation
The following points are used to generate a cubic interpolating polynomial

x0 = 3 f(x0) = 0.3333
x1 = 4 f(x1) = 0.25
x2 = 5 f(x2) = 0.2
x3 = 6 f(x3) = 0.1667
The polynomial can be generated in a number of ways including: Newton’s Divided

Difference Interpolating Polynomials .

xl f(xl) First Second Third
3 0.3333 -0.0833 0.01665 -0.0027
4 0.25 -0.05 0.00835
5 0.2 -0.0333
6 0.1667

The result is:

f3(x) = 0.3333+(x-3)( -0.0833)+(x-3)(x-4)( 0.01665)+(x-3)(x-4)(x-5)( -0.0027)

If we simplify the above expression, we get:

f3 ( x)  0.943 0.3261833x  0.0491x 2  0.00271667x3

Second step: Roots finding
The roots problem can then be developed by setting this polynomial equal to the desired value
of 0.23:
f3 ( x)  0.943 0.3261833x  0.0491x 2  0.00271667x3  0.23
0  0.713  0.3261833x  0.0491x2  0.00271667x3

Root finding in MATLAB

 fzero(f,x0): tries to find a point x where f(x) = 0.

% Define the function.

>> func=inline('  0.00271667*x^3  0.0491*x^2  0.3261833*x + 0.713')
% Solve the equation f(x) = 0.
>> fzero(func, 0)

Spline Method of Interpolation

Spline method was introduced to solve one of the drawbacks of the polynomial interpolation.
In fact, when the order (n) becomes large, in many cases, oscillations appear in the resulting
polynomial.

1. Linear spline interpolation

Given x0 , y 0 , x1 , y1 ,......, x n 1 , y n 1  x n , y n  , fit linear splines to the data. This simply
involves forming the consecutive data through straight lines. So if the above data is given in
an ascending order, the linear splines are given by  yi  f ( xi ) 

Linear splines
Then, the linear spline can be expressed as:
f ( x1 )  f ( x0 ) x0  x  x1
f ( x)  f ( x0 )  ( x  x0 ),
x1  x0

f ( x2 )  f ( x1 ) x1  x  x2
 f ( x1 )  ( x  x1 ),
x2  x1
 ......
f ( xn )  f ( xn1 ) xn 1  x  xn
 f ( xn1 )  ( x  xn1 ),
xn  xn1

f ( xi )  f ( xi 1 )
where : the first divided differenceof f(x).
xi  xi 1
2. Quadratic Splines
In these splines, a quadratic polynomial approximates the data between two consecutive data
points. The splines are given by

Given (n+1) points, there are 3(n) coefficients to find and hence, need to determine 3n
equations to solve:

1. Each quadratic spline goes through two consecutive data points:

ai xi2  bi xi  ci  f ( xi ) i  1,2,, n
ai xi21  bi xi 1  ci  f ( xi 1 ) i  1,2,, n

This condition gives 2n equations as there are ‘n’ quadratic splines going through two
consecutive data points.

2. The first derivatives at the interior points must be equal:

2ai xi  bi  2ai 1 xi  bi 1 i  1,2,, n  1

Since there are (n-1) interior points, we have (n-1) such equations.
 3. Now, the total number of equations is (2n)  (n  1)  (3n  1) equations. We can
assume that the first spline (or the last spline) is linear, that is:

a0  0 (or an  0 )

Example:
Construct a quadratic spline interpolating the following data:
xi -1 0 1

 -1 0 1
yi  sin( xi )
2
Solution: Assume that the quadratic spline S2 x  is of the form:

P x   a1 x 2  b1 x  c1 , ha x   1; 0
S 2 x    1
P2 x   a2 x  b2 x  c2 , ha x  0; 1
2

Since S2 x  interpolates the data points, we have

1) S 2  1  P1  1  f  1  a1  b1  c1  1
2) S 2 0  P1 0  f 0  c1  0
3) S 2 0  P2 0  f 0  c2  0
4) S 2 1  P2 1  f 1  a2  b2  c2  1 .
In addition the continuity of S 2 at the interior point, P10  P20 gives:
Pix   2ai x  bi , i  1,2
5) P10  P20  b1  b2
Above are the five linear equations about the six unknown coefficients a1, b1, c1, a2 , b2 , c2 ,
which have infinitely many solutions. If we impose one more constrain, for example,
S2  1  0 ,
6) P1 1  2a1   1  b1  0
Then there are six linear equations for six unknown coefficients, from which it can be solved
that
6)  b1  2a1
1)  a1  2a1  a1  1

We have a1  1 and b1  2 .
5)  b1  b2  2
4) a2  2  1  a2  1 .
Therefore the quadratic spline function is:
 x 2  2 x, ha x   1; 0
S 2 x    2
 x  2 x, ha x  0; 1 .

3. Cubic Interpolating Splines

To derive a mathematical model of a cubic spline, suppose the data are:

x1 , f ( x1 ),......, xn1 , f ( xn1 ), xn , f ( xn ) .

A cubic spline S 3 ( x ) is a piecewise cubic polynomial. This means that:

s1 ( x)  a1  b1 ( x  x1 )  c1 ( x  x1 ) 2  d1 ( x  x1 )3 , x1  x  x2

 

S3 ( x)si ( x)  ai  bi ( x  xi )  ci ( x  xi ) 2  di ( x  xi )3 , xi  x  xi 1
 

sn 1 ( x)  an 1  bn 1 ( x  xn 1 )  cn 1 ( x  xn 1 ) 2  d n 1 ( x  xn 1 )3 , xn 1  x  xn

For S 3 ( x ) to be a cubic interpolating spline, we have a set on interpolation conditions: that is,

at an interval [ xi , xi 1 ]

1. si ( xi )  yi , i  1,, n  1

2. si ( xi1 )  yi1 , i  1,, n  1

at each interior point, we continuity of S 3

3. si( xi 1 )  si1 ( xi 1 ), i  1,, n  2

at each interior point, we continuity of S 3

4. si( xi1 )  si1 ( xi1 ), i  1,, n  2

Since each of the (n  1) cubic pieces has four unknown coefficients, our description of the

function S 3 ( x ) involves 4(n  1) unknown coefficients. Assuring continuity of this function

and its first and second derivatives imposes 2(n  1)  2(n  2) linear constraints on its
coefficients. Therefore, there are a total of 4n  6 linear constraints on the 4(n  1) unknown
coefficients. In order that we have the same number of equations as unknowns, we need 2
more (linear) constraints. There are various ways of specifying two additional constraints such
as the following:
 1. The so–called natural cubic spline results when we impose the conditions:
s1'' ( x1 )  0, sn'' 1 ( xn )  0.
2. Instead of imposing the natural cubic spline conditions, we could use the correct
second derivative values:
s1'' ( x1 )  c1, sn'' 1 ( xn )  c2 where c1 and c2 are given.
3. In the complete cubic spline, the slope conditions
s1 ( x1 )  c1 , sn 1 ( xn )  c2 where c1 and c2 are given.

are imposed. These first derivative values of the data may not be readily available but
they can be replaced by accurate approximations.
Example:
Construct a cubic spline interpolating the following data:
xi -1 0 1

 -1 0 1
yi  sin( xi )
2
Solution:
Assume that the cubic spline S  x is of the form:

Since S2 x  interpolates the data points, we have:

The first and second derivative (i = 1; 2) are:

The continuity of the first and second derivative at each interior point, we have:


From the data we known that the function is periodic ( yi  sin( xi ) ), then the periodicity
2
may be used to specify the boundary conditions:
Therefore the cubic spline function is:

Piecewise Interpolation in MATLAB

MATLAB has several built-in functions to implement piecewise interpolation. The first is
spline:

yy=spline(x, y, xx)

This performs cubic spline interpolation.

Example:

• Generate data:

>> x = linspace(-1, 1, 9);

>> y = 1./(1+25*x.^2);

• Calculate 100 model points and determine not-a-knot interpolation

>> xx = linspace(-1, 1,100);

>> yy = spline(x, y, xx);

 Calculate actual function values at model points and plot data points, the 9-point
(solid), and the actual function (dashed),

>> yr = 1./(1+25*xx.^2);
>> plot(x, y, 'o', xx, yy, '-', xx, yr, '--')
MATLAB’s interp1 function
 While spline can only perform cubic splines, MATLAB’s interp1 function can
perform several different kinds of interpolation:

yi = interp1(x, y, xi, ‘method’)

1. x & y contain the original data

2. xi contains the points at which to interpolate
3. ‘method’ is a string containing the desired method:

‘nearest’ nearest neighbor interpolation

‘linear’ connects the points with straight lines
‘spline’ not-a-knot cubic spline interpolation
‘pchip’ or ‘cubic’ piecewise cubic Hermite interpolation

Example: 1. Find the linear splines for the following data set

i 1 2 3 4 5
x 0 5 7 8 10
y 0 2 -1 -2 20
>> x=[0 5 7 8 10];
>> y=[0 2 -1 -2 20];
>> n=1000; % Number of points.
>> xi=linspace(0,10,n); % Create vector of x values
>> yi=interp1(x,y,xi,'linear'); % Calculate y
>> plot(x,y,'o',xi,yi);

Example: 2. Find the quadratic splines for the following data set

i 1 2 3 4 5 6 7 8
x 0 0.2 0.6 1 1.4 1.6 1.8 2
y 1 1 1.2 1.1 0.6 0.5 0.5 0.5

>> x=[0 0.2 0.6 1 1.4 1.6 1.8 2];

>> y=[1 1 1.2 1.1 0.6 0.5 0.5 0.5];
>> n=1000; % Number of points.
>> xi=linspace(0,2,n); % Create vector of x values
>> yi=interp1(x,y,xi,'spline'); % Calculate yi

>> plot(x,y,'o',xi,yi);