Spline Interpolation for Polynomial Approximation

Sec:18.6 Spline Interpolation

Sec:18.6 Spline Interpolation Polynomial with degree 20 Sections 18.1 and 18.2 concerned the approximation of arbitrary functions on closed intervals using a single polynomial. However, high-degree polynomials can oscillate erratically, that is, a minor fluctuation over a small portion of the interval can induce large fluctuations over the entire range. natural cubic spline An alternative approach is to divide the approximation interval into a collection of subintervals and construct a different approximating polynomial on each subinterval. This is called piecewise- polynomial approximation.

Sec:18.6 Spline Interpolation Linear Splines We start with n+1 nodes (knots): ??< ??< < ?? We have n subintervals: [??,??] ,[??,??] , ,[?? ?,??] ??(? ??) + ?(?0) ??(? ??) + ?(?1) [??,??] [??,??] [?? ?,??] ? ? = ??=? ??+? ?(??) ??+? ?? ?? ?(? ?? ?) + ?(?? 1)

Sec:18.6 Spline Interpolation Cubic Splines The most common piecewise-polynomial approximation uses cubic polynomials between each successive pair of nodes and is called cubic spline. Example Construct a cubic spline that passes through the points (1, 2), (2, 3), and (3, 5). ?0= 1 ?1= 2 ?2= 3 ?0(?) = ??+ ??? 1 + ??(? 1)2+??(? 1)3 ?1(?) = ??+ ??? 2 + ??(? 2)2+??(? 2)3 There are 8 constants to be determined, which requires 8 conditions. Four conditions come from the fact that the splines must agree with the data at the nodes. Two more come from The final two come from the natural boundary conditions ?0 (2) = ?1 (2) ?0 (1) = 0 ?0 (2) = ?1 (2) ?0(1) = 2 ?1(2) = 3 ?1 (3) = 0 ?0(2) = 3 ?1(3) = 5

Sec:18.6 Spline Interpolation Example Construct a cubic spline that passes through the points (1, 2), (2, 3), and (3, 5). ?0(?) = ??+ ??? 1 + ??(? 1)2+??(? 1)3 ?1(?) = ??+ ??? 2 + ??(? 2)2+??(? 2)3 ?0=2, ?1=3 ?0=0 1 1 0 0 0 0 0 1 1 1 1 3 -1 0 0 0 6 0 -2 0 0 0 0 2 6 Four conditions come from the fact that the splines must agree with the data at the nodes. ?0 ?0 ?1 ?1 ?1 1 2 0 0 0 = ?1(3) = 5 ?0(1) = 2 ?0(2) = 3 ?1(2) = 3 Two more come from ?? ?? ?? ?? ?0 (2) = ?1 (2) 0.2500 0 0.7500 2.0000 ?0 (2) = ?1 (2) The final two come from the natural boundary conditions ?? ?? ?? ?? ?0 (1) = 0 -0.2500 0.7500 1.5000 3.0000 ?1 (3) = 0

Sec:18.6 Spline Interpolation Example Construct a cubic spline that passes through the points (1, 2), (2, 3), and (3, 5). Matlab built-in command: spline , interp1 , csape clear ; clc; close all x=[1 2 3]; y=[2 3 5]; cs = csape(x,y,'second'); cs.coefs xq=1:0.1:3; plot(x,y,'ro',xq,ppval(cs,xq),'-'); grid on ?? ?? ?? ?? 0.2500 0 0.7500 2.0000 ?0(?) = ??+ ??? 1 + ??(? 1)2+??(? 1)3 ?? ?? ?? ?? ?1(?) = ??+ ??? 2 + ??(? 2)2+??(? 2)3 -0.2500 0.7500 1.5000 3.0000

Sec:18.6 Spline Interpolation ezplot and plot ezplot('sin(x)',[0,6.4]) grid on; x=[0:0.1:8.4]; y=sin(x); figure plot(x,y,'-r') grid on

Sec:18.6 Spline Interpolation Cubic Splines The most common piecewise-polynomial approximation uses cubic polynomials between each successive pair of nodes and is called cubic spline. ??? ? = 0,1,2, ,? 1 ??(?) = ??+ ??? ?? + ??(? ??)2+??(? ??)3 ? [??,??+1] There are 4n constants to be determined, which requires 4n conditions.

Sec:18.3 Spline Interpolation Given (n+1) data points Cubic Splines: There are 4n constants to be determined, which requires 4n conditions. n+1 1) interpolating data n-1 2) continuity at interior points n-1 3) continuous slope at interior points n-1 4) continuous curvature at interior points 2 5) natural spline

Sec:18.3 Spline Interpolation Cubic Splines: 1 2 3 4 5

Sec:18.3 Spline Interpolation ? ? = ?

Sec:18.3 Spline Interpolation ? ? = ?

Sec:18.3 Spline Interpolation Cubic Splines: Algorithm: (Natural Cubic Spline) 1 2 3 Calculate: Assemble the matrix ? Assemble the vector ? 4 solve the linear system ? ? = ? 5

Sec:18.3 Spline Interpolation The matrix ? is tridiagonal ??, ??, ?? ?, ?? ?] Above the Main diagonal : [?, Main diagonal : ? ??+ ??, ? ??+ ??, ? ?? ?+ ?? ?, ?] [?, ??, ??, ?? ?, ?] [??, below the Main diagonal :

Sec:18.3 Spline Interpolation Example Algorithm: (Natural Cubic Spline) Construct a cubic spline that passes through the points (1, 2), (2, 3), and (3, 5). 1 2 3 Calculate: ? is 3X3 1= 1 0= 1 ? = 2 Assemble the matrix ? Above: (0 1) ? ? ? ? ? ? ? ? ? (1 4 1) Main diagonal : Assemble the vector ? ? = below : (1 0) 4 solve the linear system ? ? = ? ? ? ? ? = 5 0 ? ? ? ? ? ? ? ? ? ?? ?? ?? 0 3 0 = ? = 3/4 0 1/4 1/4 2 3 5 3/4 3/2 ? = ? = ? = Above the Main diagonal : ?? ?, ?? ?] [?, ??, ??, Main diagonal : [?, ? ??+ ??, ? ??+ ??, ? ?? ?+ ?? ? ?] below the Main diagonal : [??, ??, ??, ?? ?, ?]

Sec:18.3 Spline Interpolation Example 18.14 page 522 ? = # ?? ???????? = 5 ? is 6X6 Develop cubic splines for the data 3= 2 4= 1 2= 2 1= 1 1 1 0= 1 0 4 1 0 30 81 168 171 0 x f (x) 3 6 19 99 291 444 1 2 3 5 7 8 1 6 2 2 8 2 ? = 2 6 0 1 1 Above the Main diagonal : ?? ?, ?? ?] [?, ??, ??, [ 0 1 2 2 1 ] Main diagonal : [?, ? ??+ ??, ? ??+ ??, [ 1 4 6 8 6 1 ] ? ?? ?+ ?? ? ?] below the Main diagonal : [??, ??, ??, ?? ?, ?] 0 [ 1 1 2 2 0 ] 5.4236 8.3057 12.8712 24.2096 0

Sec:18.3 Spline Interpolation x=[1 2 3 5 7 8]; y=[3 6 19 99 291 444]; n=length(y)-1; for i=1:n h(i) = x(i+1)-x(i); end v=zeros(n+1,1); for i=2:n v(i) = 3*( (y(i+1)-y(i))/h(i) - (y(i)-y(i-1))/h(i-1) ); end A = zeros(n+1,n+1); diagA = [1 4 6 8 6 1]; a_diag = [0 1 2 2 1]; b_diag = [1 1 2 2 0]; A = diag(diagA) + diag(a_diag,1) + diag(b_diag,-1) c = A\v cs = csape(x,y,'second'); cs.coefs

Sec:18.3 Spline Interpolation function [a,b,c,d,coeff] = myspline(x,y) Here we calculate spline coefficients using myspline.m and store them in coeff matrix. Also we use MATLAB built-in command csape to calculate the coefficients and store them in cof. n1=length(y); n = n1-1; for i=1:n; h(i) = x(i+1)-x(i); end v=zeros(n+1,1); for i=2:n v(i) = 3*( (y(i+1)-y(i))/h(i) - (y(i)-y(i-1))/h(i-1) ); end md(1) = 1; ud(1)=0; for i=2:n md(i) = 2*(h(i-1)+h(i)); ud(i) = h(i); end md = [md , 1]; for i=1:n-1 dd(i) = h(i); end dd = [ dd 0]; clear, clc x=[1 2 3 5 7 8]; y=[3 6 19 99 291 444]; [a,b,c,d,coeff] = myspline(x,y); cs = csape(x,y,'second'); cof=cs.coefs coeff A = zeros(n+1,n+1); A = diag(md) + diag(ud,1) + diag(dd,-1); c = A\v; for i=1:n+1; a(i) = y(i); end for i=1:n d(i) = ( c(i+1) - c(i) )/(3*h(i)); b(i) = (a(i+1)-a(i))/h(i) - h(i)*( 2*c(i) + c(i+1))/3; end coeff =[d' c(1:n) b' a(1:n)']; end

Sec:18.6 Spline Interpolation Quadratic Splines The objective in quadratic splines is to derive a second-order polynomial for each interval between data points. Example Construct a cubic spline that passes through the points (1, 2), (2, 3), and (3, 5). ?0= 1 ?1= 2 ?2= 3 ?0(?) = ??+ ??? + ???2 ?1(?) = ??+ ??? + ???2 There are 6 constants to be determined, which requires 6 conditions. Four conditions come from the fact that the splines must agree with the data at the nodes. The first derivatives at the interior knots must be equal the second derivative is zero at the first point. ?0(1) = 2 ?0 (1) = 0 ?1(2) = 3 ?0 (2) = ?1 (2) ?0(2) = 3 ?1(3) = 5

Sec:18.6 Spline Interpolation Parametric Curves None of the techniques developed in this section can be used to generate curves of the form shown in the adjacent Figure. because this curve cannot be expressed as a function of one coordinate variable in terms of the other. A straightforward parametric technique for determining a polynomial or piecewise polynomial to connect the points in the order given is to use a parameter t xi 1 0 1 0 1 yi 0 1 0.5 0 1 idx=[ 0 1 2 3 4]; t=[0 0.25 0.5 0.75 1]; x=[-1 0 1 0 1]; y=[ 0 1 0.5 0 -1]; ? = ? ? , ? = ?(?) i ti xi 1 0 1 0 1 yi 0 1 0.5 0 1 0 1 2 3 4 0 0.25 0.5 0.75 1 csx = csape(t,x,'second'); csy = csape(t,y,'second'); tq=0:0.01:1; xq=ppval(csx,tq); yq=ppval(csy,tq); plot(x,y,'ro',xq,yq,'-'); grid on

Sec:18.6 Spline Interpolation Parametric Curves Try this - Place your hand on the computer screen. Use the mouse to select a few dozen data points outlining your hand. You might find it easier to trace your hand on a piece of paper and then put the paper on the computer screen. You should be able to see the ginput cursor through the paper. - Find parametric splines for the x and y data you have collected and then plot these splines to get a smooth representation of your hand. http://homepages.gac.edu/~hvidsten/cours es/MC355/Labs-Projects/Lab3/Lab3.html

Sec:18.6 Spline Interpolation MULTIDIMENSIONAL INTERPOLATION Data points X = -3.0000 -1.5000 0 1.5000 3.0000 -3.0000 -1.5000 0 1.5000 3.0000 -3.0000 -1.5000 0 1.5000 3.0000 -3.0000 -1.5000 0 1.5000 3.0000 -3.0000 -1.5000 0 1.5000 3.0000 Y = -3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 0 0 0 0 0 1.5000 1.5000 1.5000 1.5000 1.5000 3.0000 3.0000 3.0000 3.0000 3.0000 [X,Y,V] = peaks(5); [Xq,Yq] = meshgrid(-3:.1:3,-3:.1:3); Vq = interp2(X,Y,V,Xq,Yq,'spline'); mesh(X,Y,V); figure mesh(Xq,Yq,Vq) V = 0.0001 0.0042 -0.2450 -0.0298 -0.0000 -0.0005 0.3265 -5.6803 -0.4405 0.0036 -0.0365 -2.7736 0.9810 3.2695 0.0331 -0.0031 0.4784 7.9966 1.1853 0.0044 0.0000 0.0312 0.2999 0.0320 0.0000

Sec:18.6 Spline Interpolation MULTIDIMENSIONAL INTERPOLATION Data points X = -3.0000 -1.5000 0 1.5000 3.0000 -3.0000 -1.5000 0 1.5000 3.0000 -3.0000 -1.5000 0 1.5000 3.0000 -3.0000 -1.5000 0 1.5000 3.0000 -3.0000 -1.5000 0 1.5000 3.0000 Y = -3.0000 -3.0000 -3.0000 -3.0000 -3.0000 -1.5000 -1.5000 -1.5000 -1.5000 -1.5000 0 0 0 0 0 1.5000 1.5000 1.5000 1.5000 1.5000 3.0000 3.0000 3.0000 3.0000 3.0000 V = 0.0001 0.0042 -0.2450 -0.0298 -0.0000 -0.0005 0.3265 -5.6803 -0.4405 0.0036 -0.0365 -2.7736 0.9810 3.2695 0.0331 -0.0031 0.4784 7.9966 1.1853 0.0044 0.0000 0.0312 0.2999 0.0320 0.0000 [X,Y,V] = peaks(5); [Xq,Yq] = meshgrid(-3:.1:3,-3:.1:3); Vq = interp2(X,Y,V,Xq,Yq,'spline'); mesh(X,Y,V); figure mesh(Xq,Yq,Vq)