1. Numerical Computation Lecture 13: Cubic Splines United International College

2. Today • We will cover: • Cubic Splines • Readings: • Pav, section 6.2 (omit 6.2.1)

3. Cubic Spline Definition: A function S is a cubic spline on [a,b] if The domain of S is [a,b] S is continuous on [a,b] S’ is continuous on [a,b] S’’ is continuous on [a,b] There is a partition of points on [a,b] such that S is a polynomial of degree <= 3 on each sub-interval [ti , ti+1 ]

4. Degree k Spline Definition: A function S is a spline of degree k on [a,b] if The domain of S is [a,b] S is continuous on [a,b] S’, S’’, …, S(k-1) are all continuous on [a,b] There is a partition of points on [a,b] such that S is a polynomial of degree <= k on each sub-interval [ti , ti+1 ]

5. Computing Degree k Splines If there are n-1 knots then we have n (degree k) polynomials to find: S0 (t), S1 (t), …, Sn-1 (t) Each polynomial Si (t) is defined by (k+1) constants: ck tk + ck-1 tk-1 + ck-2 tk-2 + … + c1 t + c0 Thus, to find the spline we need to solve for n*(k+1) constants and so need n*(k+1) equations. … t0 t1 t2 t3 tn-2 tn-1 tn S0 S1 Sn-2 Sn-1

6. Computing Degree k Splines At the data points we have the following n+1 equations: S0 (t0) = y0, S1 (t1) = y1, … , Sn-1 (tn-1) = yn-1, Sn-1 (tn) = yn The continuity of S at the knots requires (n-1) equations S0 (t1) = y1 , S1 (t2) = y2, …, Sn-2 (tn-1) = yn-1 We have a total of 2n equations from the continuity of S. (Remember we need n*(k+1) equations) … t0 t1 t2 t3 tn-2 tn-1 tn S0 S1 Sn-2 Sn-1

7. Computing Degree k Splines The continuity of the (k-1) derivatives S’, S’’, …, S(k-1) at the (n-1) internal knots t1 ,t2 ,…, tn-1 produces (k-1)*(n-1) more equations. So, we have a total of 2n+(k-1)*(n-1) = 2n + kn – n –(k-1) = n*(k+1) – (k-1) equations. We have a total of n*(k+1) – (k-1) equations. We need n*(k+1) equations, so we are (k-1) short! Solution: Add (k-1) conditions on constants. (usually on derivatives) “Natural Spline”– set some derivatives=0 at t0 and/or t1 … t0 t1 t2 t3 tn-2 tn-1 tn S0’ S1’ Sn-2’ Sn-1’

8. Computing Natural Cubic Splines Example: A cubic spline has degree k=3. Suppose we have this data:Then, n+1=3, k+1=4. So, we need to determine n*(k+1) = 8 constants. Continuity of S: S0 3= -a+b–c+d, S1 h = -1, 3= 8e+4f+2g+h Continuity at knots: S0 d = -1

9. Computing Natural Cubic Splines Example: So far, have Continuity of S’ at internal knot: S0’ c= g Continuity of S’’ at internal knot: S0’’ 2b = 2f

10. Computing Natural Cubic Splines Example: Now, we have 3 = -a+b-c-1 3 = 8e+4b+2c-1 The “natural” spline has S’’(t0)=0=S’’(t2) So, -6a+2b=0 and 12e+2f=12e+2b=0. Thus, we have 4 equations in 4 unknowns: -a+b-c = 4 8e+4b+2c = 4 -6a+2b=0 12e+2b=0 Solve by Gaussian Elimination

11. Computing Natural Cubic Splines Example: Get

12. Computing Natural Cubic Splines General Case: Given we have to find n cubic polynomials: P0 (x) = a0 (x-t0 )3 + b0 (x-t0 )2 + c0 (x-t0 ) + d0 P1 (x) = a1 (x-t1 )3 + b1 (x-t1 )2 + c1 (x-t1 ) + d1 . . . Pn-1 (x) = an-1 (x-tn-1 )3 + bn-1 (x-tn-1 )2 + cn-1 (x-tn-1 ) + dn-1

13. Computing Natural Cubic Splines Solving for Constants: 1) Pk (tk) = yk -> dk = yk So, Pk(x) = ak(x-tk)3 + bk(x-tk)2 + ck(x-tk) + yk 2) Continuity at knots gives Pk(tk+1) = ak(tk+1-tk)3 + bk(tk+1-tk)2 + ck(tk+1-tk) + yk = yk+1 Thus, ak(tk+1-tk)3 + bk(tk+1-tk)2 + ck(tk+1-tk) = (yk+1 - yk) Divide by (tk+1-tk) to get ak(tk+1-tk)2 + bk (tk+1-tk) + ck = (yk+1 - yk)/ (tk+1-tk) Let hk = (tk+1-tk) and let δk = (yk+1 - yk)/ (tk+1-tk) Then, akhk2 + bk hk + ck = δk

14. Computing Natural Cubic Splines Solving for Constants: 3) Continuity at knots for P’ gives Pk’(tk+1) = 3ak(tk+1-tk)2 + 2bk(tk+1-tk) + ck = ck+1 Thus, 3akhk2 + 2bk hk + ck = ck+1 4) Continuity at knots for P’’ gives Pk’’(tk+1) = 6ak(tk+1-tk) + 2bk = 2bk+1 Thus, 6akhk + 2bk = 2bk+1 5) Finally, assume a “natural” condition: P0’(t0) = 0, so b0 = 0 and P0’’(t0) = 0, so c0 = 0 Using 2) above we get a0 = δ0 / h02

15. Computing Natural Cubic Splines Algorithm for Cubic Splines: Given Find Pk(x) = ak(x-tk)3 + bk(x-tk)2 + ck(x-tk) + dk (k=0,…,n-1) Setup: hk = (tk+1-tk)δk = (yk+1 - yk)/ (tk+1-tk) Initialization step: dk = yka0=δ0/h02 b0 = 0=c0 Iteration: ck+1=3akhk2 + 2bk hk + ck bk+1 = 3akhk + bk ak+1 = (δk+1- ck+1– bk+1 hk+1 )/(hk+12)

16. Matlab Cubic Splines function v = piececubic(x,y,u) % v = piececubic(x,y,u) finds the piecewise cubic value of S(x) n = length(x); a = zeros(n-1,1); b = zeros(n-1,1); c = zeros(n-1,1); h = diff(x); delta = diff(y)./diff(x); a(1) = delta(1)/(h(1)^2); b(1) = 0; c(1) = 0; for i = 1:n-2 c(i+1) = 3*a(i)*h(i)^2 + 2*b(i)*h(i) + c(i); b(i+1) = 3*a(i)*h(i) + b(i); a(i+1) = ( delta(i+1) - c(i+1) - b(i+1)*h(i+1) )/ ( h(i+1)^2 ); end % Find subinterval index k so that x(k) <= u < x(k+1) k = 1; for j = 2:n-1 if x(j) <= u; k = j; end end % Evaluate spline at u v = a(k)*((u-x(k))^3)+ b(k)*((u-x(k))^2) + c(k)*(u-x(k)) + y(k);

17. Plotting Cubic Splines % x,y data in vector form x = [0 1 4]; y = [1 -2 1]; % A series of values along the x-axis u = 0.0:.05:4.0; [m n] = size(u); % polyvals stores the linear spline values for the u-values splinevals = zeros(n); for i = 1:n % Compute each polynomial value splinevals(i) = piececubic(x,y,u(i)); end % Plot the x-y data as circles ('o') and the polynomial data as '-' plot(x,y,'o',u,splinevals,'-')

18. Plotting Cubic Splines