MA/CS 321
Fall 2019
1 Goal
Programming Project 2
Due: Nov 25th. NO late submissions are allowed.
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, among other parameters, their positions at different
times as a matrix of the form:
t0 x0 y0
t1 x1 y1
. . . (1) …
tn xn yn
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: Lagrange Interpolation
One option that we have is to construct a single polynomial that interpolates all n nodes, i.e., the Lagrange interpolation polynomial. More specifically, you have to code a routine that uses the (forward) divided difference method to compute the coefficients of the Lagrange polynomial corresponding to a general set of interpolation nodes (ti,f(ti)), and another routine that evaluates the Lagrange polynomial using the nested multiplication, i.e., a generalization of Horner’s method with various points:
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 Lagrange 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 Lagrange 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
t0 x[t0] y[t0]
. . . . (2)
tn 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
s0(t) t∈[t0,t1] S(t)= . . .
sn−1(t) t∈[tn−1,tn] 1
MA/CS 321 Fall 2019
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:
(ii) Clamped boundary condition:
S′′(t0) = S′′(tn) = 0. S′(t0) = y0′ and S′(tn) = yn′ .
Tofindthevaluesofbi,ci anddi,wedenotemi =S′′(ti)fori=0,…n,andhi =ti+1−ti fori= 0, 1, . . . , n − 1. After applying the interpolation conditions and continuity of the spline at the interior nodes, we
have
yi+1 − yi hi(mi+1 + 2mi) b= − ,
ih6 i
where
ci = mi , 2
m−m di=i+1 i.
6hi
The second derivatives m0, . . . , mn at nodes can be obtained by solving the following tridiagonal linear system
2 μ0 0 ··· ··· ··· 0 λ02μ10······0 m0
0 λ 2 μ ··· ··· 0
6
r0 y2−y1
y1−y0
Am⃗ =⃗r
(3)
12 m1h1+h0h1−h0 ……. ….0..
A= ,m⃗= . ,and⃗r= .
Here
. y−y y −y
. ………0 m 6 nn−1 n−1n−2
n−1 hn−1+hn−2 m
hn−1 − hn−2 rn
.
. λ2μn
n−2 n−1 0 ··· ··· ··· 0 λn−1 2
μi = hi , λi−1 =1−μi = hi−1 for i=1,…,n−1. hi−1 + hi hi−1 + hi
In addition, if the natural boundary condition is applied, then
μ0 =λn−1 =0, r0 =rn =0, m0 =mn =0.
h20 h0 hn−1 h2n−1
Again, we want to find the corresponding cubic splines for the x and y coordinates respectively and save all
coefficients in the form:
t0 x0 bx0 cx0 dx0 y0 by0 cy0 dy0 t1 x1 bx1 cx1 dx1 y1 by1 cy1 dy1
. . . . . . . . .. (4) ………
tn−1 xn−1 bxn−1 cxn−1 dxn−1 yn−1 byn−1 cyn−1 dyn−1
As long as the coefficients are available, we evaluate Sx and Sy at the test points τ0 . . . τm−1 to find the path.
If a clamped boundary condition is applied, then
μ0 =λn−1 =1, r0 = 6(y1−y0)− 6y0′, rn = 6 yn′ −6yn−yn−1.
2
MA/CS 321 Fall 2019
2.2 Routines
You have to code:
1. A routine divdif.m that implements Algorithm 3.2 in the textbook:
function [P, Q]=divdif(x,y)
Here P stores all coefficients of Lagrange polynomial, and Q is a matrix with entries Fi,j = f[xi−j,…,xi].
2. A routine
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.
3. A routine
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.
4. A routine
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.
5. 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”.
6. A MATLAB script file main.m that calls all aforementioned routines, and executes all the testings required in Summary.
3
MA/CS 321 Fall 2019
3 Summary
You are expected to submit to Canvas a .zip file containing:
1. (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 Lecture 19 for the cubic spline.
2. 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 and the clamped cubic splines for the interpolation points “ip” in the “data.mat”(refer to Theorem 3.13 in 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: Include the plots of the paths that you obtain with the interpolation nodes “ip”, and the two sets of test points described above. Save the path matrices and the coefficient matrices obtained from the four tests into a MATLAB data file “results.mat”. To name those output variables, follow the rule: “path#&” and “coeff#&” where “#” is 1 or 2 indicating the approach and “&” is “a” or “b” indicating the test point set, e.g., “path2a” represents the path obtained by using the natural cubic spline interpolant with the test point set (a).
4