Difference between revisions of "Interpolation"

From Sutherland_wiki
Jump to: navigation, search
m
 
(21 intermediate revisions by the same user not shown)
Line 1: Line 1:
 
== Linear Interpolation ==
 
== Linear Interpolation ==
−
A linear function may be written as
+
 
−
<math>f(x) = a_0 + a_1 x.</math>
+
A linear function may be written as <math>f(x) = a_0 + a_1 x.</math>
−
Given any two data points, <math>( x_1, y_1 )</math>, and <math>(x_2,y_2)</math>, we can determine a linear function that exactly passes through these points. We can do this by solving for <math>a_0</math> and <math>a_1</math>. Since there are two unknowns we require two equations. They are:
+
Given any two data points, <math>( x_1, y_1 )</math>, and
−
<center><math>\begin{align}
+
<math>(x_2,y_2)</math>, we can determine a linear function that
−
y_1 &= a_0 + a_1 x_1 \\
+
exactly passes through these points. We can do this by solving for
−
y_2 &= a_0 + a_1 x_2
+
<math>a_0</math> and <math>a_1</math>. Since there are two unknowns we
−
\end{align}</math></center>
+
require two equations. They are:
−
These may be written in [[LinearAlgebra#Linear_Systems_of_Equations|matrix form]] as
+
:<math>
−
<center><math>
+
\begin{align}
−
\left[ \begin{array}{cc} 1 & x_1 \\ 1 & x_2 \end{array} \right]
+
  y_1 &= a_0 + a_1 x_1 \\
−
\left( \begin{array}{c} a_0 \\ a_1 \end{array} \right)
+
  y_2 &= a_0 + a_1 x_2
−
=
+
\end{align}
−
\left( \begin{array}{c} y_1 \\ y_2 \end{array} \right)
+
</math>
−
</math></center>
+
 
−
Given values for <math>( x_1, y_1 )</math>, and <math>(x_2,y_2)</math>, these equations may be easily [[Linear_Systems_in_Matlab|solved in MATLAB]]. But since they are so simple, we can easily by hand to obtain a general solution,
+
These may be written in
−
<center><math>\begin{align}
+
[[Linear_Algebra#Linear_Systems_of_Equations|matrix form]] as
 +
:<math>
 +
\left[ \begin{array}{cc}
 +
  1 & x_1 \\
 +
  1 & x_2
 +
\end{array} \right]
 +
\left( \begin{array}{c}
 +
  a_0 \\ a_1
 +
\end{array} \right) =
 +
\left( \begin{array}{c}
 +
  y_1 \\ y_2
 +
\end{array} \right)
 +
</math>
 +
 
 +
Given values for <math>( x_1, y_1 )</math>, and
 +
<math>(x_2,y_2)</math>, these equations may be easily
 +
[[Linear_Systems_in_Matlab|solved in MATLAB]]. But since they are so
 +
simple, we can solve them by hand to obtain a general solution,
 +
:<math>
 +
\begin{align}
 
  a_0 &= \frac{x_1 y_2 -x_2 y_1}{x_1-x_2} \\
 
  a_0 &= \frac{x_1 y_2 -x_2 y_1}{x_1-x_2} \\
−
  a_1 &= \frac{y_1-y_2}{x_1-x_2}
+
  a_1 &= \frac{y_1-y_2}{x_1-x_2} \end{align}
−
\end{align}</math></center>
+
</math>
−
This gives the coefficients, which may then be substituted into the original equation <math>f(x) = a_0 + a_1 x</math> and simplified to obtain
+
 
−
<center><math>
+
This gives the coefficients, which may then be substituted into the
 +
original equation <math>f(x) = a_0 + a_1 x</math> and simplified to
 +
obtain
 +
:<math>
 
   f(x) = \left( \frac{y_2-y_1}{x_2-x_1} \right) \left( x-x_1 \right) + y_1.
 
   f(x) = \left( \frac{y_2-y_1}{x_2-x_1} \right) \left( x-x_1 \right) + y_1.
−
</math></center>
+
</math>
 
This is a very convenient equation for performing linear interpolation.
 
This is a very convenient equation for performing linear interpolation.
  
−
=== Example ===
+
 
−
The density of air at various temperatures and atmospheric pressure is given in the following table
+
=== Example ===  
 +
 
 +
The density of air at various temperatures and atmospheric pressure is
 +
given in the following table
  
 
{| border="1" cellpadding="3" cellspacing="0" align="center" style="text-align:center"
 
{| border="1" cellpadding="3" cellspacing="0" align="center" style="text-align:center"
 
|-
 
|-
 
! Temperature (&deg;F)  
 
! Temperature (&deg;F)  
−
| 1.39 || 1.16 || 0.99 || 0.87 || 0.78 || 0.69 || 0.63 || 0.58 || 0.54 || 0.50 || 0.46 || 0.44 || 0.41
+
| -9.7 || 80 || 170 || 260 || 350 || 440 || 530 || 620 || 710 || 800 || 890 || 980 || 1070
 
|-
 
|-
 
! Density (kg/m<sup>3</sup>)
 
! Density (kg/m<sup>3</sup>)
−
| -10 || 80 || 170 || 260 || 350 || 440 || 530 || 620 || 710 || 800 || 890 || 980 || 1070
+
| 1.39 || 1.16 || 0.99 || 0.87 || 0.78 || 0.69 || 0.63 || 0.58 || 0.54 || 0.50 || 0.46 || 0.44 || 0.41
 
|}
 
|}
  
−
Estimate the density of air at 32 &deg;F. Using linear interpolation, we have <math>\rho(T)=a_0 + a_1 T</math>.
+
Estimate the density of air at 32 &deg;F. Using linear interpolation,
−
We define <math>(T_1,\rho_1)=(-9.7,1.39)</math> and <math>(T_2,\rho_2)=(80,1.161)</math>. Now we can calculate  
+
we have <math>\rho(T)=a_0 + a_1 T</math>. We define
−
 
+
(T<sub>1</sub>, &rho;<sub>1</sub>) = (-9.7, 1.39) and
−
:<math>\begin{align}
+
(T<sub>2</sub>, &rho;<sub>2</sub>) = (80, 1.161). Now we can calculate
 +
:<math>
 +
\begin{align}
 
   \rho(32) &\approx \frac{\rho_2 -\rho_1}{T_2-T_1}\left( T - T_1 \right) + \rho_1 \\
 
   \rho(32) &\approx \frac{\rho_2 -\rho_1}{T_2-T_1}\left( T - T_1 \right) + \rho_1 \\
 
           &= \frac{1.16-1.39}{-9.7-80} \left( 32 + 9.7 \right) + 1.39 \\
 
           &= \frac{1.16-1.39}{-9.7-80} \left( 32 + 9.7 \right) + 1.39 \\
 
           &= 1.28.
 
           &= 1.28.
−
  \end{align}</math>
+
  \end{align}
 +
</math>
 +
 
  
 
==== Implementation in Matlab ====
 
==== Implementation in Matlab ====
−
There are three ways to do interpolation in MATLAB.
+
 
 +
There are three ways to do linear interpolation in MATLAB.
 
<ol>
 
<ol>
−
 
+
<li> Form and solve the linear system. This is a simple 2x2 system since we only need to solve for the slope and intercept (two coefficients)
−
<li> Solve the linear system.
 
 
<source lang="matlab">
 
<source lang="matlab">
 
  T1=-9.7;  T2=80;
 
  T1=-9.7;  T2=80;
Line 61: Line 90:
 
</source>
 
</source>
  
−
<li> Use the general equation we derived above:
+
<li> Use the general equation we derived above.  This is basically solving the linear system analytically and using the result to do the interpolation.
 
<source lang="matlab">
 
<source lang="matlab">
 
  T1=-9.7; T2-80;
 
  T1=-9.7; T2-80;
Line 69: Line 98:
 
</source>
 
</source>
  
−
<li> Use MATLAB's built-in interpolation tool
+
<li> Use MATLAB's built-in interpolation tool, '''<tt>interp1</tt>''':
 
<source lang="matlab">
 
<source lang="matlab">
 
  Ti = [-9.7 80];
 
  Ti = [-9.7 80];
Line 77: Line 106:
 
</ol>
 
</ol>
  
−
All three of these methods provide exactly the same answer. However, using MATLAB's '''interp1''' function allows you to use a vector for the points that you want to interpolate to. For example, if we wanted the density at 30, 50, 70, 10, and 110 &deg;F, we could accomplish this by:
+
All three of these methods provide exactly the same answer. However,
 +
using MATLAB's '''interp1''' function allows you to use a vector for
 +
the points that you want to interpolate to. For example, if we wanted
 +
the density at 30, 50, 70, 10, and 110 &deg;F, we could accomplish
 +
this by:
 
<source lang="matlab">
 
<source lang="matlab">
 
   Ti = [-9.7 80 170 ];
 
   Ti = [-9.7 80 170 ];
Line 85: Line 118:
 
</source>
 
</source>
  
−
==== Implementation in Excel ====
 
−
{{Stub}|section}}
 
  
 +
== Polynomial Interpolation ==
 +
 +
Polynomial interpolation is a simple extension of linear
 +
interpolation. In general, an <var>n<sup>th</sup></var> degree
 +
polynomial is given as
 +
:<math>
 +
f(x)=\sum_{i=0}^n a_i x^i.
 +
</math>
  
 +
If <var>n=1</var> then we recover a first-degree polynomial, which is
 +
linear. The formula gives <math>f(x)=\sum_{i=0}^1 a_i x^i = a_0 + a_1
 +
x</math>.
  
 +
In general, if we want to interpolate a set of data using an
 +
<var>n<sup>th</sup></var> degree polynomial, then we must determine
 +
<var>n+1</var> coefficients, <math>a_0 \ldots a_n</math>. Therefore,
 +
we require <var>n+1</var> points to interpolate using an
 +
<var>n<sup>th</sup></var> degree polynomial.
  
−
== Polynomial Interpolation ==
+
The equations that must be solved are given as
−
Polynomial interpolation is a simple extension of linear interpolation.  In general, an <var>n<sup>th</sup></var> degree polynomial is given as <center><math>f(x)=\sum_{i=0}^n a_i x^i.</math></center>
+
:<math>
 +
\left[\begin{array}{ccccc}
 +
  1 & x_{1} & x_{1}^2 & \cdots & x_{1}^{n}\\
 +
  1 & x_{2} & x_{2}^2 & \cdots & x_{2}^{n}\\
 +
  \vdots & \vdots & \vdots & \cdots & \vdots \\
 +
  1 & x_{n+1} & x_{n+1}^2 & \cdots & x_{n+1}^{n}
 +
\end{array}\right]
 +
\left(\begin{array}{c}
 +
  a_{0}\\ a_{1}\\ \vdots\\ a_{n}
 +
\end{array}\right)
 +
=
 +
\left(\begin{array}{c}
 +
  y_{1}\\ y_{2}\\ \vdots\\ y_{n+1}
 +
\end{array}\right)
 +
</math>
  
−
If <var>n=1</var> then we recover a first-degree polynomial, which is linear.  The formula gives <math>f(x)=\sum_{i=0}^1 a_i x^i = a_0 + a_1 x</math>.
+
Solving these equations provides the values for the coefficients,
 +
<math>a_0 \ldots a_n</math>. Once we know these coefficients, we can
 +
evaluate the polynomial at any point.
  
−
In general, if we want to interpolate a set of data using an <var>n<sup>th</sup></var> degree polynomial, then we must determine <var>n+1</var> coefficients, <math>a_0 \ldots a_n</math>.  Therefore, we require <var>n+1</var> points to interpolate using an <var>n<sup>th</sup></var> degree polynomial.
 
  
−
The equations that must be solved are given as
+
=== Example ===
−
<center><math>
 
−
\left[\begin{array}{ccccc}
 
−
1 & x_{1} & x_{1}^2 & \cdots & x_{1}^{n}\\
 
−
1 & x_{2} & x_{2}^2 & \cdots & x_{2}^{n}\\
 
−
\vdots & \vdots & \vdots & \cdots & \vdots \\
 
−
1 & x_{n+1} & x_{n+1}^2 & \cdots & x_{n+1}^{n}\end{array}\right]\left(\begin{array}{c}
 
−
a_{0}\\
 
−
a_{1}\\
 
−
\vdots\\
 
−
a_{n}\end{array}\right)
 
−
=
 
−
\left(\begin{array}{c}
 
−
y_{1}\\
 
−
y_{2}\\
 
−
\vdots\\
 
−
y_{n+1}\end{array}\right)
 
−
</math></center>
 
  
−
Solving these equations provides the values for the coefficients, <math>a_0 \ldots a_n</math>.  Once we know these coefficients, we can evaluate the polynomial at any point.
+
Using a third degree (cubic) polynomial and the data in the table above,
 +
approximate the density of air at 50 &deg;F and 300 &deg;F.
  
−
=== Example ===
+
A third degree polynomial (cubic polynomial) has four coefficients,
 +
<math>a_0, \, a_1, \, a_2, \, a_3</math>. To determine these
 +
coefficients we require four equations. Assuming that we have four
 +
observations for <math>\rho(T)</math>, the equations to be solved are:
 +
:<math>
 +
\begin{align}
 +
\rho_1 &= a_0 + a_1 T_1 + a_2 T_1^2 + a_3 T_1^3 \\
 +
\rho_2 &= a_0 + a_1 T_2 + a_2 T_2^2 + a_3 T_2^3 \\
 +
\rho_3 &= a_0 + a_1 T_3 + a_2 T_3^2 + a_3 T_3^3 \\
 +
\rho_4 &= a_0 + a_1 T_4 + a_2 T_4^2 + a_3 T_4^3
 +
\end{align}
 +
</math>
  
−
Using a third degree polynomial, approximate the density of air at the following temperatures: 50, 300, 500 &deg;F.
+
Once we have solved these equations we have the polynomial
 +
coefficients and we can then evaluate the resulting polynomial at any
 +
temperature that we wish.
  
−
We require 4 points to fit a third degree polynomial, <math>a_0 \, a_1 \, a_2 \, a_3</math>. We will choose the 4 points closest to the temperature that we wish to interpolate to.
+
We still must choose what 4 points in the table to use in the above
 +
equations. We will choose the 4 points in the table that are closest
 +
to the temperature that we wish to interpolate to. Note that if we
 +
want to evaluate the density at different temperatures, we may want to
 +
use different points in the table to determine the coefficients!
  
−
==== Implementation in Matlab ====
+
Our procedure to perform the interpolation is:
−
Our procedure is:
 
 
# Determine which points we need to use.
 
# Determine which points we need to use.
 
# Set up and solve the linear system to obtain the polynomial coefficients <math>a_0 \ldots a_3</math>.
 
# Set up and solve the linear system to obtain the polynomial coefficients <math>a_0 \ldots a_3</math>.
 
# Evaluate the polynomial to obtain the temperature.
 
# Evaluate the polynomial to obtain the temperature.
−
<ul><ul>
 
−
<source lang="matlab">
 
−
Ti = [ -9.7  80  170  260  350  440  530  620  710  800  890  980 1070 ]';
 
−
rhoi=[ 1.39 1.16 0.99 0.87 0.78 0.69 0.63 0.58 0.54 0.50 0.46 0.44 0.41 ]';
 
  
 +
 +
==== Implementation in Matlab ====
 +
 +
The MATLAB code to obtain the density at 50 &deg;F is
 +
:{| border="1" cellpadding="5" cellspacing="0"
 +
|-
 +
|<source lang="matlab">
 +
% Set the data points from the table.
 +
Ti = [ -9.7  80  170  260  350  440  530 ...
 +
        620  710  800  890  980 1070 ]';
 +
rhoi=[ 1.39 1.16 0.99 0.87 0.78 0.69 0.63 ...
 +
      0.58 0.54 0.50 0.46 0.44 0.41 ]';
 +
 +
% Set the temperature we want to evaluate the density at.
 
T = 50;
 
T = 50;
 +
 +
% determine the set of points that we will use to construct
 +
% the interpolant. We want to use points that are closest
 +
% to the temperature that we are interested in.
 
points = 1:4;
 
points = 1:4;
 +
 +
% set up the linear system
 
A = [ ones(4,1), ones(4,1).*Ti(points), ones(4,1)*Ti(points).^2, ones(4,1)*Ti(points).^3 ];
 
A = [ ones(4,1), ones(4,1).*Ti(points), ones(4,1)*Ti(points).^2, ones(4,1)*Ti(points).^3 ];
 
b = rhoi(points);
 
b = rhoi(points);
 +
 +
% solve for the polynomial coefficients
 
a = A\b;
 
a = A\b;
  
 +
% evaluate the density.
 
rho = a(0) + a(1)*T + a(2)*T^2 + a(3)*T^3;
 
rho = a(0) + a(1)*T + a(2)*T^2 + a(3)*T^3;
 +
</source>
 +
|}
 +
 +
For <var>T=300</var> we would probably want to choose a different set
 +
of points to determine the polynomial coefficients. Therefore, we
 +
would repeat the above with the only change being:
 +
<ul>
 +
<source lang="matlab">
 +
T=300;
 +
points=4:8;
 
</source>
 
</source>
 
</ul>
 
</ul>
−
For <var>T=300</var> we would repeat the above with the only change being:
+
 
 +
=== Polynomial Interpolation in Matlab ===
 +
 
 +
Matlab has a command '''<tt>polyfit</tt>''' that can be used to fit a
 +
polynomial to data.  Its syntax is
 +
<source lang="matlab">
 +
  a = polyfit( x, y, order );
 +
</source>
 +
: '''<tt>x</tt>''' - the independent variable vector.
 +
: '''<tt>y</tt>''' - the dependent variable vector.
 +
: '''<tt>order</tt>''' - the order of the polynomial to fit the data to.
 +
 
 +
<u>'''Note:'''</u> for <math>n</math> data points, we can fit an
 +
<math>n-1^\mathrm{th}</math> order polynomial. However, Matlab will
 +
allow you to use an order that is less than <math>n-1</math>. In the
 +
case where you use an order less than <math>n-1</math>, Matlab will
 +
apply [[Regression|regression]] to determine the polynomial
 +
coefficients. In this case, the resulting polynomial will not
 +
necessarily interpolate the data points!
 +
 
 +
Once you have the polynomial coefficients from '''<tt>polyfit</tt>''',
 +
you can evaluate the polynomial easily using the
 +
'''<tt>polyval</tt>''' command:
 +
<source lang="matlab">
 +
  yinterp = polyval( a, xinterp );
 +
</source>
 +
: '''<tt>a</tt>''' - the polynomial coefficients obtained from <tt>polyfit</tt>
 +
: '''<tt>xinterp</tt>''' - the independent variable values we want to interpolate to.
 +
: '''<tt>yinterp</tt>''' - the interpolated values.
 +
 
 +
 
 +
==== Example ====
 +
 
 
<ul><source lang="matlab">
 
<ul><source lang="matlab">
−
T=300;
+
% Set the data points from the table.
−
points=4:8;
+
Ti = [ -9.7  80  170  260  350  440  530 620  710  800  890  980 1070 ]';
 +
rhoi=[ 1.39 1.16 0.99 0.87 0.78 0.69 0.63 0.58 0.54 0.50 0.46 0.44 0.41 ]';
 +
 
 +
% Set the temperature we want to evaluate the density at.
 +
T = 50;
 +
 
 +
% fit the data to a polynomial.  We will use a third order polynomial.
 +
% Therefore, we require four points.  We choose the points nearest T=50.
 +
% If we used all of the data and a third order polynomial, we would end
 +
% up with regression, not interpolation...
 +
points = 1:4;
 +
a = polyfit( x(points), y(points), 3 );
 +
 
 +
% interpolate to T=50:
 +
rho = polyval( T, a );
 +
 
 +
% print out the result:
 +
fprintf('The density at %.1f F is: %.2f kg/m^3\n',T,rho);
 
</source></ul>
 
</source></ul>
−
For <var>T=500</var> we would change
+
 
 +
 
 +
== Cubic Spline Interpolation ==
 +
 
 +
Cubic spline interpolation uses cubic polynomials to interpolate
 +
datasets. The concept is illustrated in the following figure:
 +
[[image:spline1.png|center|500px]]
 +
 
 +
The data points are connected with cubic functions, and on each
 +
interval the coefficients must be determined. If we have
 +
<math>n</math> data points then we can define <math>n-1</math> cubic
 +
functions, one between each data point. For each spline we have four
 +
coefficients to determine. Therefore, we require <math>4(n-1)</math>
 +
constraints to determine the <math>4(n-1)</math> coefficients.
 +
 
 +
The constraints may be summarized in the following table:
 +
 
 +
{| border="1" cellpadding="3" cellspacing="0" align="center" style="text-align:center"
 +
|+ Constraints for Cubic Splines
 +
|-
 +
! width="175px"|Description
 +
! Mathematical Statement
 +
! width="75px"|Number of Constraints
 +
|- valign="top"
 +
| At each data point, two splines meet.  Here the splines must be equal to the data values.
 +
| <math>f_{i}(x_{i+1}) = y_{i+1} \quad 2 \le i \le n-2</math>
 +
<math>f_{i+1}(x_{i+1}) = y_{i+1} \quad 2 \le i \le n-2</math>
 +
| <math>2(n-2)</math>
 +
|- valign="top"
 +
| At the first and last point, the function must equal the data.
 +
| <math>f_1(x_1) = y_1</math>
 +
<math>f_{n-1}(x_n)=y_n</math>
 +
| <math>2</math>
 +
|- valign="top"
 +
| At each interior point the first derivatives of adjacent splines must be equivalent.
 +
| <math>\left. \frac{\partial f_{i}}{\partial x} \right|_{i+1} = \left. \frac{\partial f_{i+1}}{\partial x} \right|_{i+1}</math>
 +
| <math>n-2</math>
 +
|- valign="top"
 +
| At each interior point, the second derivatives of adjacent splines must be equivalent.
 +
| <math>\left. \frac{\partial^2 f_{i}}{\partial x^2} \right|_{i+1} = \left. \frac{\partial^2 f_{i+1}}{\partial x^2} \right|_{i+1}</math>
 +
| <math>n-2</math>
 +
|- valign="top"
 +
| We specify the curvature at the end points.  By choosing the curvature to be equal to zero, we obtain a natural spline.
 +
| <math>\left. \frac{\partial^2 f_1}{\partial x^2} \right|_{1} = 0</math>
 +
<math>\left. \frac{\partial^2 f_{n-1}}{\partial x^2} \right|_{n} = 0</math>
 +
| <math>2</math>
 +
|}
 +
 
 +
Adding up all of the constraints, we see that we have
 +
<math>4(n-1)</math> constraints - precisely the number that we need.
 +
We could formulate equations for all of these constraints and solve
 +
them all to determine the coefficients for the cubic function on each
 +
interval. <br clear="all" />
 +
 
 +
 
 +
=== Example 1 ===
 +
 
 +
Using the data in the [[#Example|table above]], determine the density
 +
of air at the following temperatures: 60, 120, 180, 240, 300, 360,
 +
420, 480, 540 &deg;F.
 +
 
 +
==== Implementation in Matlab ====
 +
 
 +
Matlab has built-in functions for cubic spline interpolation:
 +
<source lang="matlab">
 +
  y = interp1( xi, yi, x, 'spline' );
 +
</source>
 +
: <tt>(xi,yi)</tt> are the points at which we have data defined.
 +
: <tt>x</tt> is the point(s) where we want to interpolate.
 +
: <tt>'spline'</tt> tells Matlab to interpolate using cubic splines.
 +
 
 +
Using this built-in function, we may obtain the interpolated values as follows:
 +
 
 
<ul><source lang="matlab">
 
<ul><source lang="matlab">
−
T=500;
+
% Set the data points from the table.
−
points=5:9;
+
Ti = [ -9.7  80  170  260  350  440  530 ...
 +
        620  710  800  890  980 1070 ]';
 +
rhoi=[ 1.39 1.16 0.99 0.87 0.78 0.69 0.63 ...
 +
      0.58 0.54 0.50 0.46 0.44 0.41 ]';
 +
 
 +
% set the temperatures that we want to interpolate to.
 +
T = 60:60:540;
 +
 
 +
% determine the density at the desired temperatures.
 +
rho = interp1( Ti, rhoi, T, 'spline' );
 
</source></ul>
 
</source></ul>
−
</ul>
 
  
−
== Cubic Spline Interpolation ==
+
<!-- jcs Not sure if this is really helpful...
−
{{Stub|section}}
+
=== Example 2 ===
 +
 
 +
Assume that we have four points,
 +
<math>(x_1,y_1),\,(x_2,y_2),\,(x_3,y_3),\,(x_4,y_4)</math>. This means
 +
that we have three cubic equations that form the cubic spline. We will
 +
write these equations as
 +
 
 +
:<math> \begin{align}
 +
f_1(x) &= a_1 + b_1 x + c_1 x^2 + d_1 x^3 \\
 +
f_2(x) &= a_2 + b_2 x + c_2 x^2 + d_2 x^3 \\
 +
f_3(x) &= a_3 + b_3 x + c_3 x^2 + d_3 x^3
 +
\end{align}
 +
</math>
  
−
=== Example ===
+
From the table above, the equations that we must solve are given by
  
−
==== Implementation in Matlab ====
+
{| border="1" cellpadding="3" cellspacing="0" align="center" style="text-align:left"
 +
|+ Constraints for Cubic Splines
 +
|-
 +
! width="125px"|Description !! Equations
 +
|- valign="top"
 +
| Splines must pass through the data points.
 +
|<math>
 +
\begin{array}{lcl}
 +
f_{1}(x_{1})=y_{1} & \Rightarrow & a_{1}+b_{1}x_{1}+c_{1}x_{1}^{2}+d_{1}x_{1}^{3}=y_{1}\\
 +
f_{1}(x_{2})=y_{2} & \Rightarrow & a_{1}+b_{1}x_{2}+c_{1}x_{2}^{2}+d_{1}x_{2}^{3}=y_{2}\\
 +
f_{2}(x_{2})=y_{2} & \Rightarrow & a_{2}+b_{2}x_{2}+c_{2}x_{2}^{2}+d_{2}x_{2}^{3}=y_{2}\\
 +
f_{2}(x_{3})=y_{3} & \Rightarrow & a_{2}+b_{2}x_{3}+c_{2}x_{3}^{2}+d_{2}x_{3}^{3}=y_{3}\\
 +
f_{3}(x_{3})=y_{3} & \Rightarrow & a_{3}+b_{3}x_{3}+c_{2}x_{3}^{2}+d_{2}x_{3}^{3}=y_{3}\\
 +
f_{3}(x_{4})=y_{4} & \Rightarrow & a_{3}+b_{3}x_{4}+c_{2}x_{4}^{2}+d_{2}x_{4}^{3}=y_{4}
 +
\end{array}
 +
</math>
 +
|- valign="top"
 +
| First derivatives must match at data points (where splines meet).
 +
|<math>
 +
\begin{array}{lcl}
 +
f_1^\prime(x_2) = f_2^\prime(x_2) & \Rightarrow & b_1 + 2c_1 + 3 d_1 x_2^2 = b_2 + 2c_2 x_2 + 3d_2 x_2^2 \\
 +
f_2^\prime(x_3) = f_3^\prime(x_3) & \Rightarrow & b_2 + 2c_2 + 3 d_2 x_3^2 = b_3 + 2c_3 x_3 + 3d_3 x_3^2
 +
\end{array}
 +
</math>
 +
|- valign="top"
 +
| Second derivatives must match at data points (where splines meet).
 +
|<math>
 +
\begin{array}{lcl}
 +
f_1^{\prime\prime}(x_2) = f_2^{\prime\prime}(x_2) = 0 & \Rightarrow & 2c_1 + 6d_1 x_2 = 2c_2 + 6x_2 x_2 \\
 +
f_2^{\prime\prime}(x_3) = f_3^{\prime\prime}(x_3) = 0 & \Rightarrow & 2c_2 + 6d_2 x_3 = 2c_3 + 6x_3 x_3 \\
 +
\end{array}
 +
</math>
 +
|- valign="top"
 +
| Second derivatives are zero at end points
 +
|<math>
 +
\begin{array}{lcl}
 +
f_1^{\prime\prime}(x_1) = 0 & \Rightarrow & 2c_1 + 6d_1 x_1 = 0 \\
 +
f_3^{\prime\prime}(x_4) = 0 & \Rightarrow & 2c_3 + 6d_3 x_4 = 0
 +
\end{array}
 +
</math>
 +
|}
  
 +
These equations are linear with respect to the unknown coefficients and may be
 +
written in matrix form as
 +
:<math>
 +
\left[ \begin{array}{cccccccccccc}
 +
1 & x_1 & x_1^2 & x_1^3 & 0 & 0  & 0    & 0    & 0 & 0  & 0    & 0    \\
 +
1 & x_2 & c_2^2 & x_2^3 & 0 & 0  & 0    & 0    & 0 & 0  & 0    & 0    \\
 +
0 & 0  & 0    & 0    & 1 & x_2 & x_2^2 & x_2^3 & 0 & 0  & 0    & 0    \\
 +
0 & 0  & 0    & 0    & 1 & x_3 & x_3^2 & x_3^4 & 0 & 0  & 0    & 0    \\
 +
0 & 0  & 0    & 0    & 0 & 0  & 0    & 0    & 1 & x_3 & x_3^2 & x_3^3 \\
 +
0 & 0  & 0    & 0    & 0 & 0  & 0    & 0    & 1 & x_4 & x_4^2 & x_4^3 \\
 +
0 & 1  & 2x_2  & 3x_2^2& 0 & 1  & 2x_2  & 3x_2^2& 0 & 0  & 0    & 0    \\
 +
0 & 0  & 0    & 0    & 0 & 1  & 2x_3  & 3x_3^2& 0 & 1  & 2x_3  & 3x_3^2 \\
 +
0 & 0  & 2    & 6x_2  & 0 & 0  & 2    & 6x_2  & 0 & 0  & 0    & 0    \\
 +
0 & 0  & 0    & 0    & 0 & 0  & 2    & 6x_3  & 0 & 0  & 2    & 6x_3 \\
 +
0 & 0  & 2    & 6x_1  & 0 & 0  & 0    & 0    & 0 & 0  & 0    & 0    \\
 +
0 & 0  & 0    & 0    & 0 & 0  & 0    & 0    & 0 & 0  & 2    & 6x_3 \\
 +
\end{array} \right]
 +
\left( \begin{array}{c}
 +
a_1 \\ b_1 \\ c_1 \\ d_1 \\
 +
a_2 \\ b_2 \\ c_2 \\ d_2 \\
 +
a_3 \\ b_3 \\ c_3 \\ d_3
 +
\end{array} \right)
 +
=
 +
\left( \begin{array}{c}
 +
y_1 \\ y_2 \\ y_2 \\ y_3 \\ y_3 \\ y_4 \\
 +
0 \\ 0 \\
 +
0 \\ 0 \\
 +
0 \\ 0
 +
\end{array} \right)
 +
</math>
  
 +
-->
  
 +
<!-- jcs upload a matlab script to do cubic splines. -->
  
 
== Lagrange Polynomial Interpolation ==
 
== Lagrange Polynomial Interpolation ==
 +
{{Stub|section}}
  
−
Given <math>n_p</math> points <math>(x_k,y_k),</math> the <math>n^\mathrm{th}</math> order Lagrange polynomial that interpolates these function values, <math>f(x)</math> are expressed as  
+
Given <math>n_p</math> points <math>(x_k,y_k),</math> the
−
<center> <math> f(x)=\sum_{k=0}^{n} y_{k} L_{k}(x), </math> </center>
+
<math>n^\mathrm{th}</math> order Lagrange polynomial that interpolates
 +
these function values, <math>f(x)</math> are expressed as
 +
:<math>
 +
f(x)=\sum_{k=0}^{n} y_{k} L_{k}(x),
 +
</math>
 
where <math>L_{k}(x)</math> is the Lagrange polynomial given by
 
where <math>L_{k}(x)</math> is the Lagrange polynomial given by
−
<center> <math> L_{k}(x)=\prod_{\stackrel{i=0}{i\ne k}}^{n}\frac{x-x_{i}}{x_{k}-x_{i}}, </math> </center>
+
:<math>
−
and <math>n_p = n+1</math>. In other words, for an <math>n^\mathrm{th}</math> order interpolation, we require <math>n_p=n+1</math> points.
+
  L_{k}(x)=\prod_{\stackrel{i=0}{i\ne k}}^{n}\frac{x-x_{i}}{x_{k}-x_{i}},
 +
</math>
 +
and <math>n_p = n+1</math>. In other words, for an
 +
<math>n^\mathrm{th}</math> order interpolation, we require
 +
<math>n_p=n+1</math> points.
  
−
The <math>\prod</math> operator represents the continued product, and is analogous to the <math>\sum</math> operator for summations. For example, <math>\prod_{i=1}^{3} (x+i) = (x+1)(x+2)(x+3)</math>.
+
The <math>\prod</math> operator represents the continued product, and
 +
is analogous to the <math>\sum</math> operator for summations. For
 +
example, <math>\prod_{i=1}^{3} (x+i) = (x+1)(x+2)(x+3)</math>.
  
 
=== Example ===
 
=== Example ===

Latest revision as of 09:15, 26 August 2009

Linear Interpolation

A linear function may be written as f(x) = a_0 + a_1 x. Given any two data points, ( x_1, y_1 ), and (x_2,y_2), we can determine a linear function that exactly passes through these points. We can do this by solving for a_0 and a_1. Since there are two unknowns we require two equations. They are:

\begin{align} y_1 &= a_0 + a_1 x_1 \\ y_2 &= a_0 + a_1 x_2 \end{align}

These may be written in matrix form as

\left[ \begin{array}{cc} 1 & x_1 \\ 1 & x_2 \end{array} \right] \left( \begin{array}{c} a_0 \\ a_1 \end{array} \right) = \left( \begin{array}{c} y_1 \\ y_2 \end{array} \right)

Given values for ( x_1, y_1 ), and (x_2,y_2), these equations may be easily solved in MATLAB. But since they are so simple, we can solve them by hand to obtain a general solution,

\begin{align} a_0 &= \frac{x_1 y_2 -x_2 y_1}{x_1-x_2} \\ a_1 &= \frac{y_1-y_2}{x_1-x_2} \end{align}

This gives the coefficients, which may then be substituted into the original equation f(x) = a_0 + a_1 x and simplified to obtain

f(x) = \left( \frac{y_2-y_1}{x_2-x_1} \right) \left( x-x_1 \right) + y_1.

This is a very convenient equation for performing linear interpolation.


Example

The density of air at various temperatures and atmospheric pressure is given in the following table

Temperature (°F) -9.7 80 170 260 350 440 530 620 710 800 890 980 1070
Density (kg/m3) 1.39 1.16 0.99 0.87 0.78 0.69 0.63 0.58 0.54 0.50 0.46 0.44 0.41

Estimate the density of air at 32 °F. Using linear interpolation, we have \rho(T)=a_0 + a_1 T. We define (T1, ρ1) = (-9.7, 1.39) and (T2, ρ2) = (80, 1.161). Now we can calculate

\begin{align} \rho(32) &\approx \frac{\rho_2 -\rho_1}{T_2-T_1}\left( T - T_1 \right) + \rho_1 \\ &= \frac{1.16-1.39}{-9.7-80} \left( 32 + 9.7 \right) + 1.39 \\ &= 1.28. \end{align}


Implementation in Matlab

There are three ways to do linear interpolation in MATLAB.

  1. Form and solve the linear system. This is a simple 2x2 system since we only need to solve for the slope and intercept (two coefficients) 
     T1=-9.7;  T2=80;
 rho1=1.39; rho2=1.16;
 T = 32;
 A = [ 1 T1; 1 T2 ];
 b = [ rho1; rho2 ];
 a = A\b;
 rho = a(1) + a(2)*T;
  2. Use the general equation we derived above. This is basically solving the linear system analytically and using the result to do the interpolation. 
     T1=-9.7; T2-80;
 rho1=1.39; rho2=1.16;
 T = 32;
 rho = (rho2-rho1)/(T2-T1) * (T-T1) + rho1;
  3. Use MATLAB's built-in interpolation tool, interp1: 
     Ti = [-9.7 80];
 rhoi = [1.39 1.16];
 rho = interp1( Ti, rhoi, 32 );

All three of these methods provide exactly the same answer. However, using MATLAB's interp1 function allows you to use a vector for the points that you want to interpolate to. For example, if we wanted the density at 30, 50, 70, 10, and 110 °F, we could accomplish this by: 

  Ti = [-9.7 80 170 ];
  rhoi = [1.39 1.16 0.995];
  T = 30:20:110;
  rho = interp1( Ti, rhoi, T );


Polynomial Interpolation

Polynomial interpolation is a simple extension of linear interpolation. In general, an nth degree polynomial is given as

f(x)=\sum_{i=0}^n a_i x^i.

If n=1 then we recover a first-degree polynomial, which is linear. The formula gives f(x)=\sum_{i=0}^1 a_i x^i = a_0 + a_1 x.

In general, if we want to interpolate a set of data using an nth degree polynomial, then we must determine n+1 coefficients, a_0 \ldots a_n. Therefore, we require n+1 points to interpolate using an nth degree polynomial.

The equations that must be solved are given as

\left[\begin{array}{ccccc} 1 & x_{1} & x_{1}^2 & \cdots & x_{1}^{n}\\ 1 & x_{2} & x_{2}^2 & \cdots & x_{2}^{n}\\ \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & x_{n+1} & x_{n+1}^2 & \cdots & x_{n+1}^{n} \end{array}\right] \left(\begin{array}{c} a_{0}\\ a_{1}\\ \vdots\\ a_{n} \end{array}\right) = \left(\begin{array}{c} y_{1}\\ y_{2}\\ \vdots\\ y_{n+1} \end{array}\right)

Solving these equations provides the values for the coefficients, a_0 \ldots a_n. Once we know these coefficients, we can evaluate the polynomial at any point.


Example

Using a third degree (cubic) polynomial and the data in the table above, approximate the density of air at 50 °F and 300 °F.

A third degree polynomial (cubic polynomial) has four coefficients, a_0, \, a_1, \, a_2, \, a_3. To determine these coefficients we require four equations. Assuming that we have four observations for \rho(T), the equations to be solved are:

\begin{align} \rho_1 &= a_0 + a_1 T_1 + a_2 T_1^2 + a_3 T_1^3 \\ \rho_2 &= a_0 + a_1 T_2 + a_2 T_2^2 + a_3 T_2^3 \\ \rho_3 &= a_0 + a_1 T_3 + a_2 T_3^2 + a_3 T_3^3 \\ \rho_4 &= a_0 + a_1 T_4 + a_2 T_4^2 + a_3 T_4^3 \end{align}

Once we have solved these equations we have the polynomial coefficients and we can then evaluate the resulting polynomial at any temperature that we wish.

We still must choose what 4 points in the table to use in the above equations. We will choose the 4 points in the table that are closest to the temperature that we wish to interpolate to. Note that if we want to evaluate the density at different temperatures, we may want to use different points in the table to determine the coefficients!

Our procedure to perform the interpolation is:

  1. Determine which points we need to use.
  2. Set up and solve the linear system to obtain the polynomial coefficients a_0 \ldots a_3.
  3. Evaluate the polynomial to obtain the temperature.


Implementation in Matlab

The MATLAB code to obtain the density at 50 °F is 

% Set the data points from the table.
Ti = [ -9.7   80  170  260  350  440  530 ...
        620  710  800  890  980 1070 ]';
rhoi=[ 1.39 1.16 0.99 0.87 0.78 0.69 0.63 ...
       0.58 0.54 0.50 0.46 0.44 0.41 ]';

% Set the temperature we want to evaluate the density at.
T = 50;

% determine the set of points that we will use to construct
% the interpolant. We want to use points that are closest
% to the temperature that we are interested in.
points = 1:4;

% set up the linear system
A = [ ones(4,1), ones(4,1).*Ti(points), ones(4,1)*Ti(points).^2, ones(4,1)*Ti(points).^3 ];
b = rhoi(points);

% solve for the polynomial coefficients
a = A\b;

% evaluate the density.
rho = a(0) + a(1)*T + a(2)*T^2 + a(3)*T^3;

For T=300 we would probably want to choose a different set of points to determine the polynomial coefficients. Therefore, we would repeat the above with the only change being: 

     T=300;
 points=4:8;

Polynomial Interpolation in Matlab

Matlab has a command polyfit that can be used to fit a polynomial to data. Its syntax is 

  a = polyfit( x, y, order );
x - the independent variable vector.
y - the dependent variable vector.
order - the order of the polynomial to fit the data to.

Note: for n data points, we can fit an n-1^\mathrm{th} order polynomial. However, Matlab will allow you to use an order that is less than n-1. In the case where you use an order less than n-1, Matlab will apply regression to determine the polynomial coefficients. In this case, the resulting polynomial will not necessarily interpolate the data points!

Once you have the polynomial coefficients from polyfit, you can evaluate the polynomial easily using the polyval command: 

  yinterp = polyval( a, xinterp );
a - the polynomial coefficients obtained from polyfit
xinterp - the independent variable values we want to interpolate to.
yinterp - the interpolated values.


Example

    % Set the data points from the table.
Ti = [ -9.7   80  170  260  350  440  530 620  710  800  890  980 1070 ]';
rhoi=[ 1.39 1.16 0.99 0.87 0.78 0.69 0.63 0.58 0.54 0.50 0.46 0.44 0.41 ]';

% Set the temperature we want to evaluate the density at.
T = 50;

% fit the data to a polynomial.  We will use a third order polynomial.
% Therefore, we require four points.  We choose the points nearest T=50.
% If we used all of the data and a third order polynomial, we would end
% up with regression, not interpolation...
points = 1:4;
a = polyfit( x(points), y(points), 3 );

% interpolate to T=50:
rho = polyval( T, a );

% print out the result:
fprintf('The density at %.1f F is: %.2f kg/m^3\n',T,rho);


Cubic Spline Interpolation

Cubic spline interpolation uses cubic polynomials to interpolate datasets. The concept is illustrated in the following figure:

Spline1.png

The data points are connected with cubic functions, and on each interval the coefficients must be determined. If we have n data points then we can define n-1 cubic functions, one between each data point. For each spline we have four coefficients to determine. Therefore, we require 4(n-1) constraints to determine the 4(n-1) coefficients.

The constraints may be summarized in the following table:

Constraints for Cubic Splines
Description Mathematical Statement Number of Constraints
At each data point, two splines meet. Here the splines must be equal to the data values. f_{i}(x_{i+1}) = y_{i+1} \quad 2 \le i \le n-2

f_{i+1}(x_{i+1}) = y_{i+1} \quad 2 \le i \le n-2

2(n-2)
At the first and last point, the function must equal the data. f_1(x_1) = y_1

f_{n-1}(x_n)=y_n

2
At each interior point the first derivatives of adjacent splines must be equivalent. \left. \frac{\partial f_{i}}{\partial x} \right|_{i+1} = \left. \frac{\partial f_{i+1}}{\partial x} \right|_{i+1} n-2
At each interior point, the second derivatives of adjacent splines must be equivalent. \left. \frac{\partial^2 f_{i}}{\partial x^2} \right|_{i+1} = \left. \frac{\partial^2 f_{i+1}}{\partial x^2} \right|_{i+1} n-2
We specify the curvature at the end points. By choosing the curvature to be equal to zero, we obtain a natural spline. \left. \frac{\partial^2 f_1}{\partial x^2} \right|_{1} = 0

\left. \frac{\partial^2 f_{n-1}}{\partial x^2} \right|_{n} = 0

2

Adding up all of the constraints, we see that we have 4(n-1) constraints - precisely the number that we need. We could formulate equations for all of these constraints and solve them all to determine the coefficients for the cubic function on each interval.


Example 1

Using the data in the table above, determine the density of air at the following temperatures: 60, 120, 180, 240, 300, 360, 420, 480, 540 °F.

Implementation in Matlab

Matlab has built-in functions for cubic spline interpolation: 

  y = interp1( xi, yi, x, 'spline' );
(xi,yi) are the points at which we have data defined.
x is the point(s) where we want to interpolate.
'spline' tells Matlab to interpolate using cubic splines.

Using this built-in function, we may obtain the interpolated values as follows: 

    % Set the data points from the table.
Ti = [ -9.7   80  170  260  350  440  530 ...
        620  710  800  890  980 1070 ]';
rhoi=[ 1.39 1.16 0.99 0.87 0.78 0.69 0.63 ...
       0.58 0.54 0.50 0.46 0.44 0.41 ]';

% set the temperatures that we want to interpolate to.
T = 60:60:540;

% determine the density at the desired temperatures.
rho = interp1( Ti, rhoi, T, 'spline' );


Lagrange Polynomial Interpolation

Warn.jpg
This section is a stub and needs to be expanded.
If you can provide information or finish this section you're welcome to do so and then remove this message afterwards.

Given n_p points (x_k,y_k), the n^\mathrm{th} order Lagrange polynomial that interpolates these function values, f(x) are expressed as

f(x)=\sum_{k=0}^{n} y_{k} L_{k}(x),

where L_{k}(x) is the Lagrange polynomial given by

L_{k}(x)=\prod_{\stackrel{i=0}{i\ne k}}^{n}\frac{x-x_{i}}{x_{k}-x_{i}},

and n_p = n+1. In other words, for an n^\mathrm{th} order interpolation, we require n_p=n+1 points.

The \prod operator represents the continued product, and is analogous to the \sum operator for summations. For example, \prod_{i=1}^{3} (x+i) = (x+1)(x+2)(x+3).

Example

Retrieved from "https://sutherland.che.utah.edu/wiki/index.php?title=Interpolation&oldid=422"