MA/CS 321 - Introduction to Numerical Methods - Programming Project 2: Interpolation
MA/CS 321 Programming Project 2
1 Goal
For ecological research purposes, a tracking device has been attached to some individuals on a migratory flock. The goal of this project is to construct the migration routes of these birds, by interpolating the signals obtained from the tracking device.
2 Detailed description
2.1 Method
The device attached to the birds sends information about, times as a matrix of the form: among other parameters, their positions at different
The pair (x(t), y(t)) describes a parametric curve corresponding to the path followed by the member of the flock, which we are interested in. To find it, we will use polynomial interpolation to get a description of the evolution of the variables x and y with respect to t, and evaluate it at different times in the interval [t0 , tn ].
2.1.1 First Approach: Polynomial Interpolation
One option that we have is to construct a single polynomial that interpolates all n nodes. More specifically, you have to code a routine that uses the (forward) divided difference method to compute the coefficients in the Newton form corresponding to a general set of interpolation nodes (ti , f (ti )), and another routine that evaluates the interpolating polynomial using the nested multiplication, i.e., a generalization of Horner’s method with various points (refer to Lecture 18): Pn (t) = (· · · (f [t0 , t1 , . . . , tn ](t − tn−1 ) + f [t0 , . . . , tn−1 ])(t − tn−2 ) · · · )(t − t0 ) + f [t0 ]. Then, a main function, which has the information of the tracking device, uses the previous routine with the first two columns (ti , xi ) to find the interpolation polynomial Pnx (t) that for each time t will give us the x coordinate of the position, and again with the first and third columns (ti , yi ) to find the interpolation polynomial Pny (t) that for each time t will give us the y coordinate. Once you have the coefficients of the polynomials, evaluate the both polynomials at the testing points τ0 , . . . , τm−1 to get the points in the path. Save the polynomial coefficients in the form
x[t0 , . . . , tn ] y[t0 , . . . , tn ]
2.1.2 Second Approach: Cubic Spline Interpolation
As you can see by plotting the results, sometimes the error is quite large near the boundaries of the interval [t0 , tn−1 ], and grows fast as we increase the number of points (see the second question in the implementation part). To handle this, we consider the cubic spline interpolant which is a piecewise function of the form t ∈ [t0 , t1 ]
sn−1 (t) t ∈ [tn−1 , tn ]
By the definition, the cubic spline S(t) satisfies the following conditions a. The function S(t) restricted on each subinterval [ti , ti+1 ], denoted by si (t), is a cubic polynomial of the form si (t) = yi + bi (t − ti ) + ci (t − ti )2 + di (t − ti )3
for i = 0, 1, . . . , n − 1.
b. (i) Natural boundary condition: S ′′ (t0 ) = S ′′ (tn ) = 0. (ii) Clamped boundary condition: S ′ (t0 ) = y0′
and S ′ (tn ) = yn′ .
To find the values of bi , ci and di , we denote mi = S ′′ (ti ) for i = 0, . . . n, and hi = ti+1 − ti for i = 0, 1, . . . , n − 1. After applying the interpolation conditions and continuity of the spline at the interior nodes, we have hi (mi+1 + 2mi yi+1 − yi The second derivatives m0 , . . . , mn at nodes can be obtained by solving the following tridiagonal linear system Am = r
for i = 1, . . . , n − 1.
In addition, if the natural boundary condition is applied, then µ0 = λn−1 = 0, r0 = rn = 0, m0 = mn = 0.
If a clamped boundary condition is applied, then Again, we want to find the corresponding cubic splines for the x and y coordinates respectively and save all coefficients in the form:
As long as the coefficients are available, we evaluate S x and S y at the test points τ0 . . . τm−1 to find the path.
2.2 Routines
You have to code:
- A routine divdif.m that finds the divided differences (refer to the Coef pseudocode in p.166 of the textbook): function [P, Q] = divdif(x,y)
Here P stores all coefficients of Newton interpolating polynomial, and Q is a matrix with entries Fi,j = f [xi−j , . . . , xi ].
- A routine evalP.m (refer to Eval pseudocode in p.166 of the textbook) function pt = evalP(P,t_int,t_test)
that evaluates the the polynomial with coefficients P associated with interpolation points t_int at the test point(s) t_test and returns that value in pt.
- A routine (refer to Spline3_Coef in p.271) function S = spline3(t,ft,varargin)
that computes the cubic splines corresponding to the interpolation points {(ti , f (ti ))} and stores all coefficients in S. Here varargin{1} specifies the boundary condition type, e.g., strings “natural” and “clamped”, and varargin{2} contains f ′ (t0 ) and f ′ (tn ) if the clamped boundary condition is used. By default, i.e., when varargin is empty, we use the natural boundary condition without any declarations. To test the clamped boundary condition using the given data, f ′ (t0 ) and f ′ (tn ) can be approximated by applying the five-point endpoint formula at t0 and tn .
- A routine (refer to Spline3_Eval pseudocode in p.271) function st = evalS(S,t_int,t_test)
that evaluates the spline with coefficients stored in S at the point(s) t_test and returns that value in st.
- A routine function [path, coeff] = findpath(ip,tp,method,varargin) that finds the path matrix of the form (1) and the coefficient matrix of the form (2) or (4). Here ip is a matrix of the form (1) storing all interpolation points, tp is a column vector storing τ0 , . . . , τm−1 , and method specifies the interpolation method, e.g., strings “lagrange” and “spline”.
- A MATLAB script file main.m that calls all aforementioned routines, and executes all the testings required in Summary.
3 Summary
You are expected to submit to Canvas the following files including:
- (50%) The source files “divdif.m”, “evalP.m”, “spline3.m”, “evalS.m”, “findpath.m” and “main.m”. To make your codes more readable, please add proper comments providing a brief description or intention of each step. Refer to Lectures 25-26 for the cubic spline1 .
- A .pdf file with a brief report of the project, including the following sections:
- (10%) User Guide: A concise description of all the routines and applications: what do they do, calling syntax, format and meaning of each input/output variable.
- (30%) Implementation:
– Provide an upper bound for the maximal error in your approximations, when using the Lagrange interpolation polynomial for the interpolation points “ip” in the “data.mat”(refer to Lectures 19-20 or Theorems 1-2 in Sec 4.2 of the textbook). Describe the tests that you have done to validate your code, the problems encountered, and their respective solutions. You must include at least the tests using the following test points as τj ’s: (a) 311 equidistant nodes in the interval [0, 3.1] generated by “(0:.01:3.1)’”; (b) 200 uniformly distributed random numbers in [0, 3.1] stored as “tp” in “data.mat”. (Hint: You can load the data by using the command “load(‘data.mat’)”.)
- (10%) Solutions: In main.m, you are expected to run the tests using the interpolation nodes “ip” and the two sets of test points described above, and plot the paths. The resulting figures should be included in the report.
Note that the notation and the position of coefficients may be different.
