Cubic Splines and Matlab

Cubic Splines and Matlab In this section, we introduce the concept of the cubic spline, and how they are implemented in Matlab. Of particular importance are the new Matlab data structures that we will see. Cubic Splines Defined Definition: Given n data points, (x1 , y1 ), . . . , (xn , yn ), a cubic spline is a piecewise-defined function of the form: S1 (x) = y1 + b1 (x − x1 ) + c1 (x − x1 )2 + d1 (x − x1 )3 S2 (x) = y2 + b2 (x − x2 ) + c2 (x − x2 )2 + d2 (x − x2 )3 .. .. . . Sn−1 (x) = yn−1 + bn−1 (x − xn−1 ) + cn−1 (x − xn−1 )2 + dn−1 (x − xn−1 )3 for x ∈ [x1 , x2 ] for x ∈ [x2 , x3 ] for x ∈ [xn−1 , xn ] Note the difference between this and in class- We can go ahead and substitute the y’s in directly. The 3n − 3 unknowns (b’s, c’s, d’s) are chosen to satisfy the interpolation constraints and also some smoothness constraints: • (Interpolation Conditions) We already have Si (xi ) = yi for i = 1 to n − 1. In addition, we want the overall function to be continuous. At the interior points, i = 1, 2, . . . , n−2, Si (xi+1 ) = Si+1 (xi+1 ) = yi+1 and at the far right value of x, Sn−1 (xn ) = yn Altogether, there are n − 1 constraints here. • (Smoothness Condition 1) At the interior points, Si 0 (xi+1 ) = Si+1 0 (xi+1 ) There are n − 2 constraints here. • (Smoothness Condition 2) At the interior points, Si 00 (xi+1 ) = Si+1 00 (xi+1 ) There are n − 2 constraints here. So far, we have (n − 1) + (n − 2) + (n − 2) = 3n − 5 constraints, but we have 3(n − 1) = 3n − 3 coefficients. We need two more constraints to define the cubics uniquely. 1 Using the Two Extra Constraints We can have different cubic splines depending on how we want to use our two extra constraints. Here are some common ones: • (The Natural Spline) S1 00 (x1 ) = 0 = Sn−1 00 (xn ) • (The Clamped Spline) The user defines S1 0 (x1 ) and Sn−1 0 (xn ) • “Not-A-Knot”: (Matlab Default) Third derivatives match at x2 and xn−1 ; That is, d1 = d2 and dn−2 = dn−1 , or S 000 1 (x2 ) = S 000 2 (x2 ) S 000 n−2 (xn−2 ) = S 000 n−1 (xn−2 ) Derivatives and Integrals We just note here the derivative and antiderivative of each piece of the cubic spline, Si (x) = yi + bi (x − xi ) + ci (x − xi )2 + di (x − xi )3 Si 0 (x) = bi + 2ci (x − xi ) + 3di (x − xi )2 and Si 00 (x) = 2ci + 6di (x − xi ) The antiderivative (note the shift): Z 1 1 1 Si (x) dx = yi (x − xi ) + bi (x − xi )2 + ci (x − xi )3 + di (x − xi )4 2 3 4 In particular, if we integrate over the interval [xi , xi+1 ] and the data is evenly spaced with width h, then: Z xi+1 1 1 1 Si (x) dx = yi h + bi h2 + ci h3 + di h4 2 3 4 xi Matlab Data Structures In Matlab, we can collect different types of data together into a single “structure”. For example, suppose that you have a vector x and a matrix A and several constants, a, b, c that you are going to be using over and over again. To wrap them up into a structure, you could write: Examp.x=[0;1;2]; Examp.A=[1 2;3 4;5 6]; Examp.constants=[1,2,3]; To access, for example, the (3, 2) element you’ve stored in the matrix A, you would type: Examp.A(3,2) 2 Piecewise Polynomials Matlab has built-in commands for dealing with piecewise-defined polynomials, like cubic splines. A piecewise-defined polynomial is defined in Matlab by a vector containing the breaks and a matrix defining the polynomial coefficients. Caution: A vector of coefficients, like [3,2,1], over an interval like [2,3] is interpreted by the piecewise polynomial routines as: 3(x − 2)2 + 2(x − 2) + 1 NOT 3x2 + 2x + 1 Example: Constructing a piecewise polynomial by hand: Suppose we want to define the following: ( P (x) = x2 − 3x + 1 3x − 5 if x ∈ [1, 2) if x ∈ [2, 3] (NOTE: In the Matlab documentation, it is unclear which function is used to evaluate the interior breaks. This is the way it is done- The intervals are all assumed to be of the form [a, b) except the last interval.) First we need to rewrite these1 : x2 − 3x + 1 = (x − 1)2 − (x − 1) − 1 and 3x − 5 = 3(x − 2) + 1 In Matlab, make the piecewise polynomial and plot the results: breaks=[1,2,3]; coefs=[1,-1,-1;0,3,1]; pp=mkpp(breaks,coefs); x=linspace(1,3); y=ppval(pp,x); plot(x,y); In Matlab, there is a command polyval that can be used to evaluate a polynomial. For example, the command: y=polyval([1,2,3,4],x,[],[3,5]); will evaluate the polynomial: (x − 3)3 (x − 3)2 (x − 3) +2 +3 +4 5 5 5 1 Is there an easy way to do this? Use a Taylor expansion based at the left endpoint. 3 1 Implementing Splines In Matlab The relevent commands are: spline ppval mkpp (Also see the help page for any of these commands). The relevent command is: pp = spline(x,y) Builds the cubic spline using the data in x, y, and outputs the result as a piecewise polynomial, pp. The structure pp will contain several pieces of data. Before continuing, let’s try one: x=0:10; y=sin(x); A=spline(x,y); If you look at A, you will see: A = form: breaks: coefs: pieces: order: dim: ’pp’ [0 1 2 3 4 5 6 7 8 9 10] [10x4 double] 10 4 1 The coefficients for the spline (highest power first) are in the array A.coefs Let’s see how we might pull this structure apart and look at functions separately. Here, we will take out the third cubic function, and plot it together with the original spline on the interval [−1, 10]. c = A.coefs(3,:); xx= linspace(-1,10); %Coefficients for cubic %Domain for plot %We could evaluate the polynomial manually: yy=c(1)*(xx-2).^3+c(2)*(xx-2).^2+c(3)*(xx-2)+c(4); %We could use Matlab’s polyval command to evaluate. y2=polyval(c,xx,[],[2,1]); %Here we evaluate the spline: 4 y3=ppval(A,xx); %Plot the spline with the third cubic: plot(xx,yy,xx,y3,’k-’); %Show that y2 and yy are the same thing: max(abs(y2-y)) Matlab Exercises 1. Use Maple to construct the spline for the data x 0 1 2 3 4 y 1 −1 2 1 0 2. How would we manipulate the coefficients to get the first, and second derivatives of a spline? Below are two sample M −files. The first will return the second derivative, and the second will return the third derivative. Your homework: Write a similar M −file named dspline that will return the first derivative. function pp=dspline(pp) % function p2 = dspline(pp) % Returns the piecewise polynomial that represents the second % derivative of the cubic spline in the data structure pp. % Example: x=0:10; y=sin(x); pp=spline(x,y); % p2=dspline(pp); % xx=linspace(0,10); % yy=ppval(p2,xx); % plot(xx,yy); A=zeros(size(pp.coefs)); A(:,3)=6*pp.coefs(:,1); A(:,4)=2*pp.coefs(:,2); pp.coefs=A; Here’s an M −file for the third derivative. function pp=d3spline(pp) % function p2 = d3spline(pp) % Returns the piecewise polynomial that represents the third % derivative of the cubic spline in the data structure pp. 5 % % % % % Example: x=0:10; y=sin(x); pp=spline(x,y); p2=d3spline(pp); xx=linspace(0,10); yy=ppval(p2,xx); plot(xx,yy); A=zeros(size(pp.coefs)); A(:,4)=6*pp.coefs(:,1); pp.coefs=A; %Initializes a matrix of zeros 3. Do the same thing as before, but implement the antiderivative of the spline. Call the function antispline, and assume the constant of integration is zero. Be sure you also change pp.order. Use your code to estimate Z 1 −1 1 dx 1 + x2 Extra: If you plot the antiderivative, you will get a very choppy curve (because each piece starts at zero). Is it possible to change this so that you get a smooth curve starting at zero? How would this code then be used to numerically estimate a definite integral? 6