Chapter 6 Data Fitting

6.1 Interpolation:

6.1.1 The Basics:

The basic idea here is the following. We are given a set of data points

x₁

···

x_n

y₁

···

y_n

for some positive integer n. This data will be interpreted as the ordered pairs corresponding to some function f which must pss through each of these data point pairs. Hence, our basic condition is that for all relevant indices i, we have

f(x_i)

y_i

6.1.2 The Vandermonde Approach:

The first choice for f is to use a polynomial of of degree one less than the number of data points. Hence, we seek a polynomial p of the form

p(x)

a₁ + a₂ x + a₃ x² + ... + a_n x^n-1

satisfying

a₁ + a₂ x₁ + a₃ x₁² + ... + a_n x₁^n-1	=	y₁
a₁ + a₂ x₂ + a₃ x₂² + ... + a_n x₂^n-1	=	y₂
	···
a₁ + a₂ x_n + a₃ x_n² + ... + a_n x_n^n-1	=	y_n

this can be reorganized into a matrix equation

1	x₁	x₁²	...	x₁^n-1
1	x₂	x₂²	...	x₂^n-1
···	···	···	···	···
1	x_n	x_n²	...	x_n^n-1

a₁

···

a_n

y₁

···

y_n

The Matrix above is the Vandermonde matrix which is generally called V.

VanderMonde Matlab Code:

Now to see an easy way to set up the Vandermonde matrix, consider the 3 × 3 case. Start with all ones:

1	1	1
1	1	1
1	1	1

Note that we can get column 2 of the Vandermonde matrix V by pointwise column vector multiplication of the data vector x and the first column of V:

x₁

x₂

x₃

x₁

x₂

x₃

This sets the second column of V. Using it, we do another pointwise column vector multiplication of the data x and the new second column of V to create the third column of V.

x₁

x₂

x₃

x₁

x₂

x₃

x₁²

x₂²

x₃²

which completes the calculation of the Vandermonde 3 × 3 matrix. So to do this in Matlab is actually easy. We just use the new operator

.*

to indicate we want to do pointwise column vector multiplication. Consider the following code:

Vander.m

 1 function a = Vander(x,y)
 2 %
 3 % x  column vector of domain data points
 4 % y  column vector of range data points
 5 %
 6 % a  column vector of coefficients for the
 7 %    interpolating polynomial p
 8 %
 9 n = length(x);
10 V = ones(n,n);
11 for j=2:n
12   %Set up column j
13   V(:,j) = x.*V(:,j-1);
14 end
15 %Use the built in
16 c = V\y 
17 %Do Each Step Explicitly
18 [L,U,P] = lu(V);
19 z = LTriSol(L,P*y);
20 a = UTriSol(U,z);
21 
22

syntax highlighted by Code2HTML, v. 0.8.8b

In this code, we have inserted two ways to do the required linear system calculations. These are the built in

%Use the built in
c = V\y

which is the same as the explicit steps

%Do Each Step Explicitly
[L,U,P] = lu(V);
z = LTriSol(L,P*y);
a = UTriSol(U,z);

You'll note that c and a are the same.

Run Time Output:

Vander.prompts

 1 >> x = [0.0; 0.5; 1.0; 6.0; 7.0; 9.0]
 2 
 3 x =
 4 
 5          0
 6     0.5000
 7     1.0000
 8     6.0000
 9     7.0000
10     9.0000
11 
12 >> y = [0.0; 1.6; 2.0; 2.0; 1.5; 0.0]
13 
14 y =
15 
16          0
17     1.6000
18     2.0000
19     2.0000
20     1.5000
21          0
22 
23 >> a = Vander(x,y);
24 
25 c =
26 
27          0
28     4.8643
29    -3.8559
30     1.1208
31    -0.1348
32     0.0057
33 
34 >> a
35 
36 a =
37 
38          0
39     4.8643
40    -3.8559
41     1.1208
42    -0.1348
43     0.0057

syntax highlighted by Code2HTML, v. 0.8.8b

Evaluating the Interpolating Polynomial:

Now that we have the desired coefficients a, we need efficient ways to evaluate the polynomial p. We'll use Horner's method.

Horner's Method in Matlab:

HornerV.m

 1 function pval = HornerV(a,z)
 2 %
 3 % a    column vector of polynomial coefficients
 4 % z    column vector of domain points to evaluate
 5 %      the polynomial at.
 6 %
 7 % pval column vector of same size at z.
 8 %      pval(i) = polynomial(z(i))
 9 %
10 n = length(a);
11 m = length(z);
12 pval = a(n)*ones(m,1);
13 for k=n-1:-1:1
14   pval = z.*pval + a(k);
15 end

syntax highlighted by Code2HTML, v. 0.8.8b

Run Time Output:

Here we show that if evaluate our newly created Vandermonde polynomial at the original data x we get back y.

HornerV.prompts

 1 >> w = HornerV(a,x)
 2 
 3 w =
 4 
 5          0
 6     1.6000
 7     2.0000
 8     2.0000
 9     1.5000
10    -0.0000

syntax highlighted by Code2HTML, v. 0.8.8b

Now to see what our polynomial looks like, let's plot it on the interval -1,10].

Plot.prompts

1 >> close all
2 >> x0 = linspace(-1,10,110)';
3 >> y0 = HornerV(a,x0);
4 >> plot(x0,y0);
5 >> print -deps plot1

syntax highlighted by Code2HTML, v. 0.8.8b

We saved the plot to a file called plot1.eps so we can include it in both our web page and our printed document. The plot looks like this:

6.1.3 Cubic Hermite Interpolation:

This special cubic intepolates both the functions values and the derivative values of the function.

Just One Interval:

Given an interval [[x_L, x_R], we wish to fit a cubic polynomial to the data (x_L,y_L,S_L) and (x_R,y_R,s_R) where the data represents the domain values, x_L and x_R; the function values, f(x_L) = y_L and f(x_R) = y_R; and derivative values f^'(x_L) = s_L and f^'(x_R) = s_R for some function f we wish to interpolate. Write the general cubic we seek in the following form:

p(x)

a + b(x-x_L) + c(x-x_L)² + d(x-x_L)²(x-x_R)

This is the Hermite cubic polynomial form. Note that this then implies that

p^'(x)

b + 2c(x-x_L) + 2d(x-x_l)(x-x_R) + d(x-x_l)²

Now let h denote the difference x_R - x_L and apply the constraints:

p(x_L)	=	y_L ® a = y_L
p' (x_L)	=	s_L ® b = s_L
p(x_R)	=	y_R ® a + bh + ch² = y_L
p' (x_R)	=	s_R ® b + 2ch + dh² = s_L

In matrix form this gives

1	0	0	0
0	1	0	0
1	h	h²	0
0	1	2h	h²

y_L

s_l

y_R

s_R

Letting D be y_r - y_L/h, we see this can easily be solved to give:

y_L

s_L

y_R - y_L - s_Lh

h²

D - s_L

s_L + s_R - 2 D

h²

The General Problem: n Intervals:

So how do we do the general problem? We start with the data triples

{(x₁,y₁,s₁), ..., (x_n,y_n,s_n) }

and define a function p which is a cubic of the form we used above on each subinterval [x_i, x_i+1] for all indices i from 1 to n-1. Hence, p is defined piecewise as

p(x)

q₁(x)	x₁ £ x £ x₂
q₂(x)	x₂ £ x £ x₃
···	···
q_n-1(x)	x_n-1 £ x £ x_n

and each q_i has the form

q_i(x)

a_i + b_i(x-x_i) + c_i(x-x_i)² + d_i(x-x_i)²(x-x_i+1)

To find these n-1 Hermite cubic polynomials, we apply the results from above for the case of one interval. Instead of using h and D like before, we will use

h_i

x_i+1 - x_i

D_i

y_i+1 - y_i

h_i

We know that the solution is given by

a_i

y_i

b_i

s_i

c_i

D_i - s_i

h_i

d_i

s_i + s_i+1 - 2 D_i

h_i²

The Basic Matlab Code:

First, we start with the code which implements the Hermite Cubic Spline idea:

pwC.m

 1 function [a,b,c,d] = pwC(x,y,s)
 2 %
 3 % x,y,s      column vectors of data of size n
 4 % a,b,c,d    column vectors of size n-1 that hold
 5 %            hermite polynomial coeficients
 6 %
 7 n = length(x);
 8 a = y(1:n-1);
 9 b = s(1:n-1);
10 Dx = diff(x);
11 Dy = diff(y);
12 D = Dy ./ Dx;
13 c = (D - s(1:n-1)) ./ Dx;
14 d = (s(2:n) + s(1:n-1) - 2*D) ./ (Dx .* Dx);

syntax highlighted by Code2HTML, v. 0.8.8b

Locating a Point in the Interval:

Now that we can compute the piecewise Hermite polynomial interploant, we need to evaluate the resulting funtion p at a point z: here is the code to locate where z is in the interval [x₁,x_n] using binary search:

Locate.m

 1 function i = Locate(x,z,g)
 2 %
 3 % x     column vector of size n which is the 
 4 %       partition x_1 < x_2 < ... < x_n
 5 % z     scalar value which is in the interval
 6 %       [x_1,x_n]
 7 % g     optional third argument which is an
 8 %       index between 1 and n-1
 9 %
10 % i     the interval that z belongs to
11 %
12 if nargin==3
13   %Try initial guess
14   if(x(g)<=z) & (z<=x(g+1))
15     i=g;
16     return;
17   end
18 end
19 
20 n = length(x);
21 if z==x(n)
22   i = n-1;
23 else
24   % Use Binary Search
25   Left = 1;
26   Right = n;
27   while Right > Left +1
28     % x(Left) <= z <= x(Right)
29     mid = floor((Left + Right)/2);
30     if z < x(mid)
31      Right = mid;
32     else
33       Left = mid;
34     end
35   end
36   i = Left;
37 end

syntax highlighted by Code2HTML, v. 0.8.8b

Evaluating a Piecewise Cubic:

Here is the code to evaluate p(z) once z has been located:

pwCEval.m

 1 function CVals = pwCEval(a,b,c,d,x,z)
 2 %
 3 % a,b,c,d     column vectors of size n-1 that store
 4 %             Hermite cubic polynomial interpolation
 5 % x           column vector of size n which is
 6 %             the partition x_1 < x_2 < ... < x_n
 7 % z           column vector of size m with each
 8 %             value $z_i$ lying in the interval
 9 %             [x_1,x_n]
10 %
11 % CVals       column vector of size m which stores the
12 %             m values of the interpolating fucntion  p,
13 %             p(z_1),...,p(z_m)
14 %
15 m = length(z);
16 CVals = zeros(m,1);
17 g = 1;
18 for j = 1:m
19   i = Locate(x,z(j),g);
20   %Use a Hormer Mutliply Method
21   CVals(j) = d(i)*(z(j)-x(i+1)) + c(i);
22   CVals(j) = CVals(j)*(z(j)-x(i)) + b(i);
23   CVals(j) = CVals(j)*(z(j) - x(i)) + a(i);
24   g = i;
25 end

syntax highlighted by Code2HTML, v. 0.8.8b

Some Run Time Results:

A function that occurs in many places in the sigmoid curve given by

s(x)

0.5(1 + tanh(

x-x₀

)

where x₀ determines the location of the inflection point of s and g determines the slope at x₀. We will use the piecewise Hermite polynomials to approximate s for x₀ = 0 and g = 0.5 on the interval [-2,2].

Here is the run time code:

Sigmoid.prompts

 1 >> z = linspace(-2,2,200)';
 2 >> fvals = 0.5+0.5*tanh(2*z);
 3 >> n = 4;
 4 >> x = linspace(-2,2,n)';
 5 >> y = 0.5+0.5*tanh(2*x);
 6 >> s = sech(2*x).*sech(2*x);
 7 >> [a,b,c,d] = pwC(x,y,s);
 8 >> Cvals = pwCEval(a,b,c,d,x,z);
 9 >> plot(z,fvals,z,Cvals,x,y,'*');
10 >> print -deps sigmoid1

syntax highlighted by Code2HTML, v. 0.8.8b

and you can see the results in the plot

We saved the plot to a file called plot1.eps so we can include it in both our web page and our printed document. The plot looks like this:

Now we will use the piecewise Hermite polynomials to approximate s for x₀ = 0 and g = 0.2 on the interval [-2,2].

Here is the run time code:

Sigmoid2.prompts

1 >> fvals = 0.5 + 0.5*tanh(5*z);
2 >> y = 0.5+0.5*tanh(5*x);
3 >> s = 2.5*sech(5*x).*sech(5*x);
4 >> [a,b,c,d] = pwC(x,y,s);
5 >> Cvals = pwCEval(a,b,c,d,x,z);
6 >> plot(z,fvals,z,Cvals,x,y,'*');
7 >> print -deps sigmoid2

syntax highlighted by Code2HTML, v. 0.8.8b

and you can see the results in the plot:

This one is not so good!

Exercises:

Consider the following function which has two rather high humps:

f(x)

(x - 0.3)² + 0.01

(x - 0.9)² + 0.04

- 6

Let's construct the piecewise Hermite cubic polynomial to this function on the interval [0,3] for the cases n = 4, n = 8 and n = 16. Use uniformly spaced points like we have done in our Matlab code examples in this section.
Let's construct the piecewise Hermite cubic polynomial to this function on the interval [0,3] using nonuniformly spaced points n = 4, n = 8 and n = 16. Use your judgement on how to space the points and comment on how and why you are choosing them. Discuss any thoughts you might have on how to automate this procedure: how could you infer from data that you are trying to interpolate a function with high humps? Provide a plot for each case of course!

6.1.4 General Cubic Splines:

The Hermite cubic splines interpolate both function and function derivative data. This makes the resulting interpolating function continuous and makes its derivative continuous too. But we don't know that the function's second derivative is continous. Another approach is to think of the points {s₁, ..., s_n} as parameters to be chosen so that the interpolating function is continuous through its second derivative. Of course, in this approach we don't match derivative data.

The Basics:

The general problem still concerns the triples

{(x₁,y₁,s₁), ..., (x_n,y_n,s_n) }

and we still define a function p which is a Hermite cubic on each subinterval [x_i, x_i+1] for all indices i from 1 to n-1. Hence, p is defined piecewise as

p(x)

q₁(x)	x₁ £ x £ x₂
q₂(x)	x₂ £ x £ x₃
···	···
q_n-1(x)	x_n-1 £ x £ x_n

and each q_i has the form

q_i(x)

a_i + b_i(x-x_i) + c_i(x-x_i)² + d(x-x_i)²(x-x_i+1)

The solutions we have derived for the Hermite cubics then give us

a_i

y_i

b_i

s_i

c_i

D_i - s_i

h_i

d_i

s_i + s_i+1 - 2 D_i

h_i²

The only difference here is that the points s_i are now being though of as parameters to be chosen. We still find the form of q_i is as follows:

q_i(x)

y_i + s_i(x-x_i) + (

D_i - s_i

h_i

)(x-x_i)²

+ (

s_i + s_i+1 - 2 D_i

h_i²

) (x-x_i)²(x-x_i+1)

q_i^'(x)

s_i + 2(

D_i - s_i

h_i

)(x-x_i) + 2(

s_i + s_i+1 - 2 D_i

h_i²

) (x-x_i)(x-x_i+1)

+ (

s_i + s_i+1 - 2 D_i

h_i²

)(x-x_i)²

q_i^{'
'}(x)

D_i - s_i

h_i

) + 2(

s_i + s_i+1 - 2 D_i

h_i²

)(x-x_i+1)

+ 4(

s_i + s_i+1 - 2 D_i

h_i²

) (x-x_i)

D_i - s_i

h_i

)

+ (

s_i + s_i+1 - 2 D_i

h_i²

) (4(x-x_i)) + 2(x-x_i+1)

A similar calculation gives us the second derivative of q_i+1

q_i+1(x)

y_i+1 + s_i+1(x-x_i+1) + (

D_i+1 - s_i+1

h_i+1

)(x-x_i+1)²

+ (

s_i+1 + s_i+2 - 2 D_i+1

h_i+1²

) (x-x_i+1)²(x-x_i+2)

q_i+1^'(x)

s_i+1 + 2(

D_i+1 - s_i+1

h_i+1

)(x-x_i+1) + 2(

s_i+1 + s_i+2 - 2 D_i+1

h_i+1²

) (x-x_i+1)(x-x_i+2)

+ (

s_i+1 + s_i+2 - 2 D_i+1

h_i+1²

) (x-x_i+1)²

q_i+1^{'
'}(x)

D_i+1 - s_i+1

h_i+1

)

+ (

s_i+1 + s_i+2 - 2 D_i+1

h_i+1²

) (4(x-x_i+1) + 2(x-x_i+2)

To get second continuity of these two cubics, it is enough to force continuity at the point that they share, x_i+1. This gives

q_i^{'
'}(x_i+1)

D_i - s_i

h_i

) + 4h_i(

s_i + s_i+1 - 2 D_i

h_i²

)

q_i+1^{'
'}(x_i+1)

D_i+1 - s_i+1

h_i+1

) + -2h_i+1(

s_i+1 + s_i+2 - 2 D_i+1

h_i+1²

)

Equating we get

D_i - s_i

h_i

) + 4h_i(

s_i + s_i+1 - 2 D_i

h_i²

)

D_i+1 - s_i+1

h_i+1

) - 2h_i+1(

s_i+1 + s_i+2 - 2 D_i+1

h_i+1²

)

h_i

(s_i + 2s_i+1 - - 3D_i)

h_i+1

(-2s_i+1 - s_i+2 - + 3D_i+1)

This simplifies to

h_i+1s_i + 2(h_i + h_i+1)s_i+1 + h_i s_I+1

3(h_i D_i+1 + h_i+1D_i)

Now these equations are satisfied for n = 1 to n-2. If we choose {s₁, ..., s_n } so satisfy the above equations, then the resulting interpolant function will have a continuous second derivative. Note we have n-2 equations in n unknowns, so there are two unresolved degress of freedom.

We get these equations:

h₂ s₁ + 2(h₁ + h₂) s₂ + h₁ s₃	=	3(h₁ D₂ + h₂ D₁)
h₃ s₂ + 2(h₂ + h₃) s₃ + h₂ s₄	=	3(h₂ D₃ + h₃ D₂)
···	=	···
h_n-1 s_n-2 + 2(h_n-2 + h_n-1) s_n-1 + h_n-2 s_n	=	3(h_n-2 D_n-1 + h_n-1 D_n-2)

It is traditional to treat the parameters s₁ and s_n are the free parameters. Move these over to the right hand side and we get

2(h₁ + h₂) s₂ + h₁ s₃	=	3(h₁ D₂ + h₂ D₁) - h₂ s₁
h₃ s₂ + 2(h₂ + h₃) s₃ + h₂ s₄	=	3(h₂ D₃ + h₃ D₂)
···	=	···
h_n-1 s_n-2 + 2(h_n-2 + h_n-1) s_n-1	=	3(h_n-2 D_n-1 + h_n-1 D_n-2) - h_n-2 s_n

The resulting n-2 by n-2 system then looks like this in matrix form:

2(h₁ + h₂)	h₁	0	0	0	···	0
h₃	2(h₂ + h₃)	h₂	0	···	0
0	h₄	2(h₃ + h₄)	h₃	0	···	0
···	···	···	···	···	···	···
0	0	···	···	0	h_n-1	2(h_n-2 + h_n-1)

s₂

s₃

···

s_n-1

3(h₁ D₂ + h₂ D₁) - h₂ s₁

3(h₂ D₃ + h₃ D₂)

···

3(h_n-3 D_n-2 + h_n-2 D_n-3)

3(h_n-2 D_n-1 + h_n-1 D_n-2) - h_n-2 s_n

EndPoint Conditions:

We can set the values of s₁ and s_n a lot of ways and the resulting interpolants we find will have different sorts of properties.

Set s₁ to be A and s_n to be B. The top and bottom equations then become

2(h₁ + h₂) s₂ + h₁ s₃ = 3(h₁ D₂ + h₂ D₁) - h₂ A

h_n-1 s_n-2 + 2(h_n-2 + h_n-1) s_n-1 = 3(h_n-2 D_n-1 + h_n-1 D_n-2) - h_n-2 B
The endpoint values are both set to be zero. These interpolants are called natural splines . This says A and B are 0 and gives the equations:

2(h₁ + h₂) s₂ + h₁ s₃ = 3(h₁ D₂ + h₂ D₁)

h_n-1 s_n-2 + 2(h_n-2 + h_n-1) s_n-1 = 3(h_n-2 D_n-1 + h_n-1 D_n-2)

Exercises:

Let's repeat the exercises for the Hermite Cubic interpolants using these new splines.

Again, we consider the following function which has two rather high humps:

f(x)

(x - 0.3)² + 0.01

(x - 0.9)² + 0.04

- 6

Let's construct the Natural Spline interpolants to this function on the interval [0,3] for the cases n = 4, n = 8 and n = 16. Use uniformly spaced points like we have done in our Matlab code examples in this section. You will need to look at the Matlab code provided in Van Loan on page 115 and on page 116. The Matlab fragments will set up the matrix you need to solve to find the natural cubic spline interpolant. You will need to modify, if necessary, the Matlab code for Location and Evaluation of the resulting function. Then plot the true function and the interpolants on the same graph like before.
This is the same except let's set A to be 1 and B to be -1 just for grins. Do the same thing as the first part. What do you think of this new choice of A and B? In general, do you think you could find a better choice of A and B than the natural spline choice?

2(h₁ + h₂) s₂ + h₁ s₃	=	3(h₁ D₂ + h₂ D₁) - h₂ A
h_n-1 s_n-2 + 2(h_n-2 + h_n-1) s_n-1	=	3(h_n-2 D_n-1 + h_n-1 D_n-2) - h_n-2 B