Lesson 12: Numerical Techniques
Objectives
- To interpolate between data points, using either linear or cubic spline models.
- To solve differential equations numerically.
Background
12.1 Interpolation
Interpolation is a technique for adding new data points within a range of a set of known data points. You can use interpolation to fill-in missing data, smooth existing data, make predictions, and more. Interpolation in MATLAB is divided into techniques for data points on a grid and scattered data points.
| Syntax | Description |
|---|---|
| vq = interp1(x,y,xq) | Returns interpolated values of a 1-D function at specific query points using linear interpolation. Vector x contains the sample points, and v contains the corresponding values, v(x). Vector xq contains the coordinates of the query points. If you have multiple sets of data that are sampled at the same point coordinates, then you can pass v as an array. Each column of array v contains a different set of 1-D sample values. |
| vq = interp1(x,y,xq,method) | Specifies an alternative interpolation method: 'linear', 'nearest', 'next', 'previous', 'pchip', 'cubic', 'v5cubic', 'makima', or 'spline'. The default method is 'linear'. |
| vq = interp1(x,y,xq,method,extrapolation) | Specifies a strategy for evaluating points that lie outside the domain of x. Set extrapolation to 'extrap' when you want to use the method algorithm for extrapolation. Alternatively, you can specify a scalar value, in which case, interp1 returns that value for all points outside the domain of x. |
| vq = interp1(y,xq) | Returns interpolated values and assumes a default set of sample point coordinates. The default points are the sequence of numbers from 1 to n, where n depends on the shape of v:
|
| vq = interp1(y,xq,method) | Specifies any of the alternative interpolation methods and uses the default sample points. |
| vq = interp1(y,xq,method,extrapolation) | Specifies an extrapolation strategy and uses the default sample points. |
| pp = interp1(x,y,method,'pp') | Returns the piecewise polynomial form of v(x) using the method algorithm. |
Linear Interpolation
| Function | Description | Example |
|---|---|---|
| interp1(x,y,xq) interp1(x,y,xq,'linear') |
Linear interpolation, which is the default. | f1 = interp1(x,y,3.5); f1 = interp1(x,y,3.5,'linear'); |
In mathematics, linear interpolation is a method of curve fitting using linear polynomials to construct new data points within the range of a discrete set of known data points.
If the two known points are given by the coordinates (x0,y0) and (x1,y1), the linear interpolant is the straight line between these points. For a value x in the interval (x0,x1), the value y along the straight line is given from the equation of slopes: \(\frac{{y - {y_0}}}{{x - {x_0}}} = \frac{{{y_1} - {y_0}}}{{{x_1} - {x_0}}}\)
Solve this equation for y, which is the unknown value at x, gives: \(y = {y_0} + \left( {x - {x_0}} \right)\frac{{{y_1} - {y_0}}}{{{x_1} - {x_0}}}\)
Interpolation of a data set
Linear interpolation creates a piecewise continuous function by connecting a set of data points with straight lines. This is useful for interpolating between the data points or for generating a set of equally spaced data points. Equally spaced data points are sometimes needed for numerical integration or digital filtering.
The following figure is the graphic using linear interpolation for the data table.
| x | y = f(x) |
|---|---|
| 0 | 0 |
| 0.560 | 10 |
| 0.805 | 60 |
| 1.01 | 30 |
| 1.16 | 40 |
| 1.36 | 30 |
| 1.50 | 60 |
| 1.66 | 70 |
| 1.88 | 80 |
To calculate the value at x = 0.5: \(y = f(0.5) = 0 + (0.5 - 0)\frac{{10 - 0}}{{0.560 - 0}} = 8.9286\)
fy = interp1(x, y, 0.5, 'linear')
fy = 8.9286
Measured Linear Interpolated Data
Define the sample points, x, and corresponding sample values, y.
x = 0:5;
y = [15 10, 9, 6, 2, 0];
To perform a single interpolation, the input to interp1 is the x data, the y data, and the new x value for which you would like an estimate of y. For example, to estimate the value of y when x is equal to 3.5, type:
interp1(x, y, 3.5)
ans = 4
You can perform multiple interpolations all at the same time by putting a vector of x-value in the third field of the interp1 function.
xq = 0:0.2:5;
yq = interp1(x, y, xq)
Which returns
yq = 15.0000 14.0000 13.0000 12.0000 11.0000 10.0000 .... 0.4000 0
We can plot the results and the original data on the same graphic.
plot(x,y,':',xq,yq,'o')
Interpolation of Coarsely Sampled Sine Function
Define the sample points, x, and corresponding sample values, y.
x = 0:pi/4:2*pi;
y = sin(x);
Define the query points to be a finer sampling over the range of x.
xq = 0:pi/16:2*pi;
Interpolate the function at the query points and plot the result
figure
vq1 = interp1(x,y,xq);
plot(x,y,'o',xq,vq1,':.');
xlim([0 2*pi]);
title('(Default) Linear Interpolation');
xlabel('x'); ylabel('y = f(x)');
To find the value at \(x = 2\)
fx = interp1(x,y,2)
fx = 0.8399
Interpolation Without Specifying Points
Define a set of function values.
close all; clear; clc;
y = [0 1.41 2 1.41 0 -1.41 -2 -1.41 0];
Define a set of query points that fall between the default points, 1:9. In this case, the default points are 1:9 because y contains 9 values.
xq = 1.5:8.5;
Evaluate y at xq.
yq = interp1(y,xq);
Plot the result.
plot((1:9),y,'o',xq,yq,'*');
legend('y','yq');
Interpolation of Complex Values
Define a set of sample points.
x = 1:10;
Define the values of the function, \(y(x) = 5x + {x^2}i\) , at the sample points.
y = (5*x)+(x.^2*1i);
Define the query points to be a finer sampling over the range of x.
xq = 1:0.25:10;
Interpolate y at the query points.
yq = interp1(x,y,xq);
Plot the real part of the result in red and the imaginary part in blue.
figure
plot(x, real(y), '*r', xq, real(yq), '-r');
hold on
plot(x, imag(y), '*b', xq, imag(yq), '-b');
legend('real', 'imag');
Cubic Spline Interpolation
Cubic Spline Interpolation
| Function | Description | Example |
|---|---|---|
| interp1(x,y,xq,'split') | Piecewise cubic spline interpolation | f3 = interp1(x, y, 3.5 , 'spliine') |
In the mathematical field of numerical analysis, spline interpolation is a form of interpolation where the interpolant is a special type of piecewise polynomial called a spline. Spline interpolation is often preferred over polynomial interpolation because the interpolation error can be made small even when using low degree polynomials for the spline. Spline interpolation avoids the problem of Runge's phenomenon, in which oscillation can occur between points when interpolating using high degree polynomials.
Measured Spline Interpolated Data
Define the sample points, x, and corresponding sample values, y.
x = 0:5;
y = [15 10, 9, 6, 2, 0];
Connecting data points with straight lines probably isn’t the best way to estimate intermediate values, although it is surely the simplest. We can create a smoother curve by using the cubic spline interpolation technique. To call the cubic spline, we need to add a fourth field to interp1 :
interp1(x, y, 3.5,'spline')
This command will return an improved estimate of y when x is equal to 3.5:
ans = 3.9417
Of course, we could also use the cubic spline technique to create an array of new estimates for y for every member of an array of x-values:
xq = 0:0.2:5;
yq = interp1(x, y, xq,'spline');
A plot of these data on the same graph as the measured data using the command
plot(x, y, xq, yq, 'o')
Extrapolation Using Cubic Spline
Extrapolate a data set to predict population growth.
Create two vectors to represent the census years from 1900 to 1990 (t) and the corresponding United States population in millions of people (p).
t = 1900:10:1990;
p = [ 75.995 91.972 105.711 123.203 131.669 ...
150.697 179.323 203.212 226.505 249.633 ];
Extrapolate and predict the population in the year 2000 using a cubic spline.
p2000 = interp1(t,p,2000,'spline')
p2000 = 270.6060
Spline Interpolation of Sine Data
Use spline to interpolate a sine curve over unevenly-spaced sample points.
x = [0 1 2.5 3.6 5 7 8.1 10];
y = sin(x);
xx = 0:.25:10;
yy = interp1(x, y, xx, 'spline');
plot(x, y, 'o', xx, yy)
Multidimensional Interpolation
Multidimensional Interpolation
interp2 - Interpolation for 2-D gridded data in meshgrid format.
| Syntax | Description |
|---|---|
| Vq = interp2(X,Y,V,Xq,Yq) | Returns interpolated values of a function of two variables at specific query points using linear interpolation. The results always pass through the original sampling of the function. X and Y contain the coordinates of the sample points. V contains the corresponding function values at each sample point. Xq and Yq contain the coordinates of the query points. |
| Vq = interp2(V,Xq,Yq) | Assumes a default grid of sample points. The default grid points cover the rectangular region, X=1:n and Y=1:m, where [m,n] = size(V). Use this syntax when you want to conserve memory and are not concerned about the absolute distances between points. |
| Vq = interp2(V) | Returns the interpolated values on a refined grid formed by dividing the interval between sample values once in each dimension. |
| Vq = interp2(V,k) | Returns the interpolated values on a refined grid formed by repeatedly halving the intervals k times in each dimension. This results in 2k-1 interpolated points between sample values. |
| Vq = interp2(___,method) | Specifies an alternative interpolation method: 'linear', 'nearest', 'cubic', 'makima', or 'spline'. The default method is 'linear'. |
| Vq = interp2(___,method,extrapval) | Also specifies extrapval, a scalar value that is assigned to all queries that lie outside the domain of the sample points. If you omit the extrapval argument for queries outside the domain of the sample points, then based on the method argument interp2 returns one of the following:
|
There are a 3-D table of data v that depends on x and y variables:
To determine the value of v, for example, at \(y = 3\) and \(x = 1.5\), you need to perform two interpolations:
- Find all v-values at \(y = 3\) by using interp1 function.
- Do a second interpolation to find v at \(x = 1.5\) from the previous v-values results.
First, define x, y, and v in MATLAB.
x = 1:4;
y = 2:2:6;
v = [ 7, 15, 22, 30;
54, 109, 165, 218;
403, 807, 1201, 1614];
Find the value of v at \(y = 3\) for all the x-values using interp1 function:
vy = interp1(y, v, 3)
vy = 30.5000 62.0000 93.5000 124.0000
The vy is all v-values at y = 3. Use interp1 function again to find the v at \(x = 1.5\)
vq1 = interp1(x, vy, 1.5)
vq1 = 46.2500
Or you can use interp2 function to solve this problem in single step:
vq2 = interp2(x, y, v, 1.5, 3)
vq2 = 46.2500
The first field in the interp2 function must be a vector defining the value associated with each column x, and the second field must be a vector defining the values associated with each row y. The array v must have the same number of columns as the number of elements in x, and must have the same number of rows as the number of elements in y. The fourth and fifth fields correspond to the values of x and of y for witch you would like to determine new v-values.
MATLAB also includes a three-dimensional interpolation function, interp3. Using help feature to get details on how to use this function and interpn, which allows you to perform n-dimensional interpolation.
3-D Data Table
Practice Exercises 12.1
Create x and y vectors to represent the following data:
- Plot the data on an x–y plot.
- Use linear interpolation to approximate the value of y when \(x = 15\).
- Use cubic spline interpolation to approximate the value of y when \(x = 15\).
- Use linear interpolation to approximate the value of x when \(y = 80\).
- Use cubic spline interpolation to approximate the value of x when \(y = 80\).
- Use cubic spline interpolation to approximate y-values for x-values evenly spaced between 10 and 100 at intervals of 2.
- Plot the original data on an x–y plot as data points not connected by a line. Also, plot the values calculated in Exercise 6.
12.2 Curve Fitting
| Syntax | Description |
|---|---|
| p = polyfit(x,y,n) | Returns the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y. The coefficients in p are in descending powers, and the length of p is n+1: \(p(x) = {p_1}{x^n} + {p_2}{x^{n - 1}} + \cdots + {p_n}x + {p_{n + 1}}\) |
| [p,S] = polyfit(x,y,n) | Also returns a structure S that can be used as an input to polyval to obtain error estimates. |
| [p,S,mu] = polyfit(x,y,n) | Also returns mu, which is a two-element vector with centering and scaling values. mu(1) is mean(x), and mu(2) is std(x). Using these values, polyfit centers x at zero and scales it to have unit standard deviation, \(\hat x = \frac{{x - \bar x}}{{{\sigma _x}}}\) |
The polyfit function returns the coefficients of a polynomial that best fits the data.
Linear Regression
Linear Regression
The simplest way to model a set of data is as a straight line. Define the sample points, x, and corresponding sample values, y.
x = 0:5;
y = [15, 10, 9, 6, 2, 0];
Using linear interpolation to plot the data.
xq = 0:0.2:5;
yq = interp1(x, y, xq)
plot(x, y, ':', xq, yq, 'o')
Using linear regression to model the data. Since straight line is a first-order polynomial, so use polyfit to fit a first-degree polynomial to the points.
p = polyfit(x,y,1);
linear_y = p(1) * x + p(2);
plot(x, y, 'o', x, linear_y)
Linear Regression Through Zero
Linear Regression Through Zero
Sometime, based on the physics of a system, we know that a model should pass through zero. For example, a force applied to a kite by blowing wind. If the wind velocity is zero, the applied force should be zero too. However, the polyfit function does not have the option to forcing the fit through the zero point.
The left-division operator, \,can employs Gaussian elimination, which can be used to solve systems of linear equations. If a system of equation is over-specified, the left-division operator uses a least-squares fit strategy similar to that used bu polyfit.
Consider the following wind Velocity-vs-Force data:
| Wind Velocity, x | Resulting Force, y |
|---|---|
| 0 | 0 |
| 10 | 24 |
| 20 | 38 |
| 30 | 64 |
| 40 | 82 |
| 50 | 90 |
We can use polyfit function to model the data:
p = polyfit(x, y, 1)
p = 1.8571 3.2381
The corresponds to the linear equation: \(y = p(1) * x + p(2) = 1.8571x + 3.2381\)
We can force the model through zero with left-division:
pz = x'\y'
pz = 1.9455
which corresponds to the linear equation: \(y = pz(1) * x = 1.9455x\)
The results are plotted in the following figure.
Higher-Order Polynomial Regression
Higher-Order Polynomial Regression
The straight lines are not the only equation that could be analysed with the regression technique. Taylor's theorem tells us that any smooth function can be approximated as a polynomial. Therefore, a common approach is to fit the data with a higher-order polynomial of the form.
\(y = f(x) = p(x) = {p_1}{x^n} + {p_2}{x^{n - 1}} + \cdots + {p_n}x + {p_{n + 1}}\)
Second- and Thrid-order Models
Consider the following data:
| x | 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|---|
| y | 15 | 10 | 9 | 6 | 2 | 0 |
Second-Order Equation
Using polyfit function to fit the data to second-order equation, and find the sum of the the squares to determine whether this models fit the data better.
p2 = polyfit(x, y, 2);
y2 = p(1)*x.^2 + p(2)*x + p(3);
sum((y2-y).^2)
ans = 3.2643
Add more data points to plot the curves defined by the new equation.
xq = 0:0.2:5;
yq2 = p2(1)*(xq.^2)+p2(2)*xq+p2(3);
plot(x, y, 'o', xq, yq2)
title('Second-Order Model');
Third-Order Equation
Using polyfit function to fit the data to third-order equation, and calculate the sum of the the squares to compare with the second-order model.
p3 = polyfit(x, y, 3);
y3 = p3(1)*x.^3 + p3(2)*x.^2 + p3(3)*x + p3(4);
sum((y3-y).^2)
ans = 2.9921
Add more data points to plot the curves defined by the new equation.
xq = 0:0.2:5;
yq3 = p3(1)*xq.^3 + p3(2)*xq.^2 + p3(3)*xq + p3(4);
plot(x, y, 'o', xq, yq3)
title('Third-Order Model');
Since the result of the sum-of-the squares calculation for the third-order model is less than the value found for the second-order model, the third-order model is better fit to the data.
Notice the slight curvature in each model. Although mathematically these models fit the data better, they may not be as good a representation of reality as the straight line. As an engineer or scientist, you’ll need to evaluate any modeling you do. You’ll need to consider what you know about the physics of the process you’re modeling and how accurate and reproducible your measurements are.
The Polyval Function
The Polyval Function — Polynomial evaluation
| Syntax | Description |
|---|---|
| y = polyval(p,x) | Evaluates the polynomial p at each point in x. The argument p is a vector of length n+1 whose elements are the coefficients (in descending powers) of an nth-degree polynomial: \(p(x) = {p_1}{x^n} + {p_2}{x^{n - 1}} + \cdots + {p_n}x + {p_{n + 1}}\) The polynomial coefficients in p can be calculated for different purposes by functions like polyint, polyder, and polyfit, but you can specify any vector for the coefficients. To evaluate a polynomial in a matrix sense, use polyvalm instead. |
| [y,delta] = polyval(p,x,S) | Uses the optional output structure S produced by polyfit to generate error estimates. delta is an estimate of the standard error in predicting a future observation at x by p(x). |
| y = polyval(p,x,[],mu) [y,delta] = polyval(p,x,S,mu) |
Use the optional output mu produced by polyfit to center and scale the data. mu(1) is mean(x), and mu(2) is std(x). Using these values, polyval centers x at zero and scales it to have unit standard deviation, \(\hat x = \frac{{x - \bar x}}{{{\sigma _x}}}\) This centering and scaling transformation improves the numerical properties of the polynomial. |
In the previous section, we use polyfit function to find the coefficients. Those coefficients are used into a MATLAB expression for the corresponding polynomial, and used it to calculate new values of y . The polyval function can perform the same job without our having to reenter the coefficients.
The polyval function requires two inputs. The first is a coefficient array, such as that created by polyfit. The second is an array of x-values for which we would like to calculate new y-values.
For example, we might have
coef = polyfit(x, y, 1);
y_first_order_fit = polyval(coef, x);
The code can be shortened to one line by nesting functions:
y_first_order_fit = polyval( polyfit(x, y, 1), x);
Fourth- and Fifth-order Models
We can use our new understanding of the polyfit and polyval functions to write a program to calculate and plot the fourth- and fifth-order fits for the data from previous section.
Use the data from the previous section:
| x | 0 | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|---|
| y | 15 | 10 | 9 | 6 | 2 | 0 |
xq = 0:0.2:5;
y4 = polyval( polyfit(x, y, 4), xq);
y5 = polyval( polyfit(x, y, 5), xq);
subplot(1,2,1)
plot(x, y, 'o', xq, y4)
axis([0,6,-5,15])
title('Fourth-Order Model');
subplot(1,2,2)
plot(x, y, 'o', xq, y5)
axis([0,6,-5,15])
title('Fifth-Order Model');
Practice Exercises 12.2
Create x and y vectors to represent the following data:
- Use the polyfit function to fit the data for \(z = 15\) to a first-order polynomial.
- Create a vector of new x values from 10 to 100 in intervals of 2. Use your new vector in the polyval function together with the coefficient values found in Exercise 1 to create a new y vector.
- Plot the original data as circles without a connecting line and the calculated data as a solid line on the same graph. How well do you think your model fits the data?
- Repeat Exercises 1 through 3 for the x and y data corresponding to \(z = 30\).
12.3 Using The Interactive Fitting Tools
12.4 Differences and Numerical Differentiation
Questions
Consider a gas in a piston–cylinder device in which the temperature is held constant. As the volume of the device was changed, the pressure was measured. The volume and pressure values are reported in the following table:
- Answer the following questions:
- Use linear interpolation to estimate the pressure when the volume is 3.8 m3.
- Use cubic spline interpolation to estimate the pressure when the volume is 3.8 m3.
- Use linear interpolation to estimate the volume if the pressure is measured to be 1000 kPa.
- Use cubic spline interpolation to estimate the volume if the pressure is measured to be 1000 kPa.
- Use the interp2 function and the data to approximate a pressure value when the volume is 5.2 m3 and the temperature is 425 K.
- Fit the data with first-, second-, third-, and fourth-order polynomials, using the polyfit function:
- Plot your results on the same graph.
- Plot the actual data as a circle with no line.
- Calculate the values to plot from your polynomial regression results at intervals of 0.2 m3.
- Do not show the calculated values on the plot, but do connect the points with solid lines.
- Which model seems to do the best job? (Hit: compare the values of the sum-of-the squares calculation of each model.)
- Using a polynomial to model a function can be very useful, but it is always dangerous to extrapolate beyond your data. We can demonstrate this pitfall by modeling a sine wave as a third-order polynomial.
- Define x = -1:0.1:1
- Calculate y = sin(x)
- Use the polyfit function to determine the coefficients of a third-order polynomial to model these data.
- Use the polyval function to calculate new values of y (modeled_y) based on your polynomial, for your x vector from -1 to 1.
- Plot both sets of values on the same graph. How good is the fit?
- Create a new x vector, new_x = -4:0.1:4.
- Calculate new_y values by finding sin(new_x).
- Extrapolate new_modeled_y values by using polyval with the coefficient vector you found in part (c), and the new_x values.
- Plot both new sets of values on the same graph. How good is the fit outside of the region from -1 to 1?