echo on %%%%% ch07_exa.m %%%%% % Run script and study the effects of each line. echo off format compact; format short; disp(sprintf('**** Examples from Chapter 07 with 5 data points ****')); t=[-2 -1 0 1 2]'; % equidistant independent variable y=[-27 -10 0 -6 10]'; % measured data values n=size(t); % number of points xx=t(1):0.01:t(n); % points for ploting disp(sprintf('**** Example 7.1 Monomial Basis *')); echo on %%%% Use polynomial interpolation with monomial basis % p(t)=x_1+x_2*t+x_3*t^2+x_4*t^3+x_5*t^4. % Find coefficients x_i. A=[1 t(1) t(1)^2 t(1)^3 t(1)^4;... 1 t(2) t(2)^2 t(2)^3 t(2)^4;... 1 t(3) t(3)^2 t(3)^3 t(3)^4;... 1 t(4) t(4)^2 t(4)^3 t(4)^4;... 1 t(5) t(5)^2 t(5)^3 t(5)^4;... ] % Vandermonde matrix x=A\y % Solve dense system to get coefficients. tt=sym('tt'); % symbolic variable tt p5=simplify (x(1)+x(2)*tt+x(3)*tt^2+x(4)*tt^3+x(5)*tt^4) % Interpolating polynomial of second order % It can be evaluated for t efficiently by Horner's algorithm. echo off % Plot data and interpolant. yy=zeros(401,1); for i=1:401 tt=xx(i); yy(i)=eval(p5); end echo on figure(1); plot(xx,yy,'b', t,y,'ro') title('Data points (-2 -27),(-1 5),(0 0),(1 -6),(2 10) and their Interpolants.') text(0.3,-10,'Monomial Interpolation for n=5','Color','b') disp(sprintf('**** Example 7.6 Cubic Spline Interpolation *')); echo on %%%% Use a piecewise interpolation with polynomials of degree 3, one for % each pair of points. %(n-1) polynomials with 4 paramethers have to be found so A has dimension %of 4*(n-1)=16 %knote values equal to the measured values A=[1 t(1) t(1)^2 t(1)^3 0 0 0 0 0 0 0 0 0 0 0 0;... 1 t(2) t(2)^2 t(2)^3 0 0 0 0 0 0 0 0 0 0 0 0;... 0 0 0 0 1 t(2) t(2)^2 t(2)^3 0 0 0 0 0 0 0 0;... 0 0 0 0 1 t(3) t(3)^2 t(3)^3 0 0 0 0 0 0 0 0;... 0 0 0 0 0 0 0 0 1 t(3) t(3)^2 t(3)^3 0 0 0 0;... 0 0 0 0 0 0 0 0 1 t(4) t(4)^2 t(4)^3 0 0 0 0;... 0 0 0 0 0 0 0 0 0 0 0 0 1 t(4) t(4)^2 t(4)^3;... 0 0 0 0 0 0 0 0 0 0 0 0 1 t(5) t(5)^2 t(5)^3;... %first derivatives in inner knots are equal 0 1 2*t(2) 3*t(2)^2 0 -1 -2*t(2) -3*t(2)^2 0 0 0 0 0 0 0 0;... 0 0 0 0 0 1 2*t(3) 3*t(3)^2 0 -1 -2*t(3) -3*t(3)^2 0 0 0 0;... 0 0 0 0 0 0 0 0 0 1 2*t(4) 3*t(4)^2 0 -1 -2*t(4) -3*t(4)^2;... %second derivatives in inner knots are equal 0 0 2 6*t(2) 0 0 -2 -6*t(2) 0 0 0 0 0 0 0 0;... 0 0 0 0 0 0 2 6*t(3) 0 0 -2 -6*t(3) 0 0 0 0;... 0 0 0 0 0 0 0 0 0 0 2 6*t(4) 0 0 -2 -6*t(4);... %second derivatives in boundary knots = 0 (normal condition) 0 0 2 6*t(1) 0 0 0 0 0 0 0 0 0 0 0 0;... 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 6*t(5);... ] b=[y(1) y(2) y(2) y(3) y(3) y(4) y(4) y(5) 0 0 0 0 0 0 0 0 ]' x=A\b %x=[4 coeficients for 4 polynomials] tt=sym('tt'); % symbolic variable tt C1=simplify(x(1)+x(2)*tt+x(3)*tt^2+x(4)*tt^3) C2=simplify(x(5)+x(6)*tt+x(7)*tt^2+x(8)*tt^3) C3=simplify(x(9)+x(10)*tt+x(11)*tt^2+x(12)*tt^3) C4=simplify(x(13)+x(14)*tt+x(15)*tt^2+x(16)*tt^3) echo off % Plot Cubic Spline interpolants C=zeros(400,1); % Space for function values for i=1:101 % First cubic polynomial tt=xx(i); C(i,1)=eval(C1); end for i=101:201 % Second cubic polynomial tt=xx(i); C(i,1)=eval(C2); end for i=201:301 % Third cubic polynomial tt=xx(i); C(i,1)=eval(C3); end for i=301:401 % Fourth cubic polynomial tt=xx(i); C(i,1)=eval(C4); end figure(1) hold on; plot(xx,C,'g') % Plot monomial again - some other solution? text(-1.8,0,'Cubic monomial Spline','Color','g') echo on disp(sprintf('**** Test - linear B-splines *')); echo on %%%% Use polynomial interpolation with cubic B-splines % p(t)=x_1*B0(t)+x_2*B1(t). % Find coefficients x_i. %Matrix A is banded, main diagonal is 1 and subdiagonals 1/4 A=diag(ones(1,n)) x=A\y % Solve dense system to get coefficients x, which are equal to y in this case. %Tabulate cubic B-spline (only ones) echo off B=zeros(201,1); % Space for linear B-spline h=t(2)-t(1); hh=-h:0.01:h; % interval for evaluation +-h for i=1:101 % left interval tt=hh(i); B(i,1)= (1/h)*(tt-t(2)); end for i=101:201 % right interval (mirror of left part if h constant) tt=hh(i); B(i,1)= (1/h)*(t(4)-tt); end % Evaluate linear B-spline by segments S=zeros(401,1); % Space for interpolated function values by linear B-spline for i=1:101 % calculate by intervals S(i,1)= x(1)*B(i+100,1) + x(2)*B(i,1); S(i+100,1)= x(2)*B(i+100,1) + x(3)*B(i,1); S(i+200,1)= x(3)*B(i+100,1) + x(4)*B(i,1); S(i+300,1)= x(4)*B(i+100,1) + x(5)*B(i,1); end figure(1) hold on; plot(xx,S,'r') % Plot monomial again - some other solution? text(-0.5,-7,'Linear B-Spline','Color','r') echo on disp(sprintf('**** Test - cubic B-splines *')); echo on %%%% Use polynomial interpolation with cubic B-splines % p(t)=x_1*B0(t)+x_2*B1(t)+x_3*B2(t)+x_4*B3(t). % Find coefficients x_i. %Matrix A is banded, main diagonal is 1 and subdiagonals 1/4 A=diag(ones(1,5))+diag(1/4*ones(1,4),1)+diag(1/4*ones(1,4),-1) x=A\y % Solve dense system to get coefficients x. %Tabulate cubic B-spline (only ones) echo off B=zeros(401,1); % Space for cubic B-spline h=t(2)-t(1); hh=-2*h:0.01:2*h; % interval for evaluation +-2*h for i=1:101 % second left interval tt=hh(i); B(i,1)= 1/(4*h^3)*(tt-t(1))^3; end for i=101:201 % left central interval tt=hh(i); B(i,1)=1/4*(1 + 3/h*(tt-t(2)) + 3/h^2*(tt-t(2))^2 - 3/h^3*(tt-t(2))^3); end for i=201:301 % right central interval tt=hh(i); B(i,1)=1/4*(1 + 3/h*(t(4)-tt) + 3/h^2*(t(4)-tt)^2 - 3/h^3*(t(4)-tt)^3); end for i=301:401 % second right interval (mirror of left part if h constant) tt=hh(i); B(i,1)= 1/(4*h^3)*(t(5)-tt)^3; end % Evaluate cubic B-spline by segments S=zeros(401,1); % Space for interpolated function values by cubic B-spline for i=1:101 % calculate by intervals S(i,1)= x(1)*B(i+200,1) + x(2)*B(i+100,1) + x(3)*B(i,1); S(i+100,1)= x(1)*B(i+300,1) + x(2)*B(i+200,1) + x(3)*B(i+100,1) + x(4)*B(i,1); S(i+200,1)= x(2)*B(i+300,1) + x(3)*B(i+200,1) + x(4)*B(i+100,1) + x(5)*B(i,1); S(i+300,1)= x(3)*B(i+300,1) + x(4)*B(i+200,1) + x(5)*B(i+100,1) ; end figure(1) hold on; plot(xx,S,'k') % Plot monomial again - some other solution? text(-1.8,-5,'Cubic B-Spline','Color','k')