Chapter 2 M-Files
Script File Index
ShowV | Illustrates InterpV and HornerV. |
ShowN | Illustrates InterpN and HornerN. |
InterpEff | Compares Vandermonde and Newton approaches. |
RungeEg | Examines accuracy of interpolating polynomial. |
TableLookUp2D | Illustrates SetUp and LinInterp2D. |
ShowPolyTools | Illustrates Matlab's polynomial tools. |
Pallas | Fits periodic data with CSInterp and CSEval. |
Function File Index
InterpV | Construction of Vandermonde interpolating polynomial. |
HornerV | Evaluates the Vandermonde interpolating polynomial. |
InterpNRecur | Recursive construction of the Newton interpolating polynomial. |
InterpN | Nonrecursive construction of the Newton interpolating polynomial. |
InterpN2 | Another nonrecursive construction of the Newton interpolating polynomial. |
HornerN | Evaluates the Newton interpolating polynomial. |
SetUp | Sets up matrix of f(x,y) evaluation. |
LinInterp2D | 2-Dimensional Linear Interpolation. |
Humps2D | A sample function of two variables. |
CSInterp | Fits periodic data with sines and cosines. |
CSEval | Evaluates sums of sines and cosines. |
% Script File: ShowV % Plots 4 random cubic interpolants of sin(x) on [0,2pi]. % Uses the Vandermonde method. close all x0 = linspace(0,2*pi,100)'; y0 = sin(x0); for eg=1:4 x = 2*pi*sort(rand(4,1)); y = sin(x); a = InterpV(x,y); pVal = HornerV(a,x0); subplot(2,2,eg) plot(x0,y0,x0,pVal,'--',x,y,'*') axis([0 2*pi -2 2]) end
% Script File: ShowN % Random cubic interpolants of sin(x) on [0,2pi]. % Uses the Newton representation. close all x0 = linspace(0,2*pi,100)'; y0 = sin(x0); for eg=1:20 x = 2*pi*sort(rand(4,1)); y = sin(x); c = InterpN(x,y); pvals = HornerN(c,x,x0); plot(x0,y0,x0,pvals,x,y,'*') axis([0 2*pi -2 2]) pause(1) end
% Script File: InterpEff % Compares the Vandermonde and Newton Approaches clc disp('Flop Counts:') disp(' ') disp(' n InterpV InterpN InterpNRecur') disp('--------------------------------------') for n = [4 8 16] x = linspace(0,1,n)'; y = sin(2*pi*x); flops(0); a = InterpV(x,y); f1 = flops; flops(0); c = InterpN(x,y); f2 = flops; flops(0); c = InterpNRecur(x,y); f3 = flops; disp(sprintf('%3.0f %7.0f %7.0f %7.0f',n,f1,f2,f3)); end
% Script File: RungeEg % For n = 7:2:15, the equal-spacing interpolants of f(x) = 1/(1+25x^2) on [-1,1]' % are of plotted. close all x = linspace(-1,1,100)'; y = ones(100,1)./(1 + 25*x.^2); for n=7:2:15 figure xEqual = linspace(-1,1,n)'; yEqual = ones(n,1)./(1+25*xEqual.^2); cEqual=InterpN(xEqual,yEqual); pValsEqual = HornerN(cEqual,xEqual,x); plot(x,y,x,pValsEqual,xEqual,yEqual,'*') title(sprintf('Equal Spacing (n = %2.0f)',n)) end
% Script File: TableLookUp2D % Illustrates SetUp and LinearInterp2D. close all clc a = 0; b = 5; n = 11; c = 0; d = 3; m = 7; fA = SetUp('Humps2D',a,b,n,c,d,m); contour(fA) x = input(sprintf('Enter x (%3.1f < = x < ="%3.1f" ):',a,b)); y = input(sprintf('Enter" y (%3.1f < ="y" < ="%3.1f" ):',c,d)); z= LinInterp2D(x,y,a,b,c,d,fA); disp(sprintf('f(x,y)="%20.16f\n',z))"
% Script File: ShowPolyTools close all x = linspace(0,4*pi); y = exp(-.2*x).*sin(x); x0 = [.5 4 8 12]; y0 = exp(-.2*x0).*sin(x0); p = polyfit(x0,y0,length(x0)-1); pvals = polyval(p,x); plot(x,y,x,pvals,x0,y0,'o') axis([0 4*pi -1 1]) set(gca,'XTick',[]) set(gca,'YTick',[]) title('Interpolating {\ite}^{-.2{\itx}}sin({\itx}) with a Cubic Polynomial')
% Script File: Pallas % Plots the trigonometric interpolant of the Gauss Pallas data. A = linspace(0,360,13)'; D = [ 408 89 -66 10 338 807 1238 1511 1583 1462 1183 804 408]'; Avals = linspace(0,360,200)'; F = CSInterp(D(1:12)); Fvals = CSeval(F,360,Avals); plot(Avals,Fvals,A,D,'o') axis([-10 370 -200 1700]) set(gca,'xTick',linspace(0,360,13)) xlabel('Ascension (Degrees)') ylabel('Declination (minutes)')
function a = InterpV(x,y) % a = InverpV(x,y) % This computes the Vandermonde polynomial interpolant where % x is a column n-vector with distinct components and y is a % column n-vector. % % a is a column n-vector with the property that if % % p(x) = a(1) + a(2)x + ... a(n)x^(n-1) % then % p(x(i)) = y(i), i=1:n n = length(x); V = ones(n,n); for j=2:n % Set up column j. V(:,j) = x.*V(:,j-1); end a = V\y;
function pVal = HornerV(a,z) % pVal = HornerV(a,z) % evaluates the Vandermonde interpolant on z where % a is an n-vector and z is an m-vector. % % pVal is a vector the same size as z with the property that if % % p(x) = a(1) + .. +a(n)x^(n-1) % then % pVal(i) = p(z(i)) , i=1:m. n = length(a); m = length(z); pVal = a(n)*ones(size(z)); for k=n-1:-1:1 pVal = z.*pVal + a(k); end
function c = InterpNRecur(x,y) % c = InterpNRecur(x,y) % The Newton polynomial interpolant. % x is a column n-vector with distinct components and y is % a column n-vector. c is a column n-vector with the property that if % % p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1)) % then % p(x(i)) = y(i), i=1:n. n = length(x); c = zeros(n,1); c(1) = y(1); if n > 1 c(2:n) = InterpNRecur(x(2:n),(y(2:n)-y(1))./(x(2:n)-x(1))); end
function c = InterpN(x,y) % c = InterpN(x,y) % The Newton polynomial interpolant. % x is a column n-vector with distinct components and y is % a column n-vector. c is a column n-vector with the property that if % % p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1)) % then % p(x(i)) = y(i), i=1:n. n = length(x); for k = 1:n-1 y(k+1:n) = (y(k+1:n)-y(k)) ./ (x(k+1:n) - x(k)); end c = y;
function c = InterpN2(x,y) % c = InterpN2(x,y) % The Newton polynomial interpolant. % x is a column n-vector with distinct components and y is % a column n-vector. c is a column n-vector with the property that if % % p(x) = c(1) + c(2)(x-x(1))+...+ c(n)(x-x(1))...(x-x(n-1)) % then % p(x(i)) = y(i), i=1:n. n = length(x); for k = 1:n-1 y(k+1:n) = (y(k+1:n)-y(k:n-1)) ./ (x(k+1:n) - x(1:n-k)); end c = y;
function pVal = HornerN(c,x,z) % pVal = HornerN(c,x,z) % Evaluates the Newton interpolant on z where % c and x are n-vectors and z is an m-vector. % % pVal is a vector the same size as z with the property that if % % p(x) = c(1) + c(2)(x-x(1))+ ... + c(n)(x-x(1))...(x-x(n-1)) % then % pval(i) = p(z(i)) , i=1:m. n = length(c); pVal = c(n)*ones(size(z)); for k=n-1:-1:1 pVal = (z-x(k)).*pVal + c(k); end
function fA = SetUp(fname,a,b,n,c,d,m) % fA = SetUp(fname,a,b,n,c,d,m) % Sets up a matrix of f(x,y) evaluations. % fname is a string that names a function of the form f(x,y). % a, b, c, and d are scalars that satisfy a<=b and c<="d." % n and m are integers>=2. % fA is an n-by-m matrix with the property that % % A(i,j) = f(a+(i-1)(b-a)/(n-1),c+(j-1)(d-c)/(m-1)) , i=1:n, j=1:m x = linspace(a,b,n); y = linspace(c,d,m); fA = zeros(n,m); for i=1:n for j=1:m fA(i,j) = feval(fname,x(i),y(j)); end end
function z = LinInterp2D(xc,yc,a,b,c,d,fA) % z = LinInterp2D(xc,yc,a,b,n,c,d,m,fA) % Linear interpolation on a grid of f(x,y) evaluations. % xc, yc, a, b, c, and d are scalars that satisfy a<=xc<=b and c<=yc<=d. % fA is an n-by-m matrix with the property that % % A(i,j)=f(a+(i-1)(b-a)/(n-1),c+(j-1)(d-c)/(m-1)) , i=1:n, j=1:m % % z is a linearly interpolated value of f(xc,yc).
[n,m]=size(fA); % xc=a+(i-1+dx)*hx 0<=dx<=1 hx=(b-a)/(n-1); i=max([1" ceil((xc-a)/hx)]); dx=(xc (a+(i-1)*hx))/hx;
% yc=c+(j-1+dy)*hy 0<=dy<=1 hy=(d-c)/(m-1); j=max([1 ceil((yc-c)/hy)]); dy=(yc (c+(j-1)*hy))/hy; z=(1-dy)*((1-dx)*fA(i,j)+dx*fA(i+1,j)) + dy*((1-dx)*fA(i,j+1)+dx*fA(i+1,j+1));
function z = Humps2D(x,y); z = 1/( (.2*(x-3)^2 +.1) + (.3*(y-1)^2 + .1) );
function F = CSInterp(f) % F = CSInterp(f) % f is a column n vector and n = 2m. % F.a is a column m+1 vector and F.b is a column m-1 vector so that if % tau = (0:n-1)'*pi/m, then % % f = F.a(1)*cos(0*tau) +...+ F.a(m+1)*cos(m*tau) + % F.b(1)*sin(tau) +...+ F.b(m-1)*sin((m-1)*tau) n = length(f); m = n/2; tau = (pi/m)*(0:n-1)'; P = []; for j=0:m, P = [P cos(j*tau)]; end for j=1:m-1, P = [P sin(j*tau)]; end y = P\f; F = struct('a',y(1:m+1),'b',y(m+2:n));
function Fvals = CSEval(F,T,tvals) % Fvals = CSEval(F,T,tvals) % F.a is a length m+1 column vector, F.b is a length m-1 column vector, % T is a positive scalar, and tvals is a column vector. % If % F(t) = F.a(1) + F.a(2)*cos((2*pi/T)*t) +...+ F.a(m+1)*cos((2*m*pi/T)*t) + % F.b(1)*sin((2*pi/T)*t) +...+ F.b(m-1)*sin((2*m*pi/T)*t) % % then Fvals = F(tvals). Fvals = zeros(length(tvals),1); tau = (2*pi/T)*tvals; for j=0:length(F.a)-1, Fvals = Fvals + F.a(j+1)*cos(j*tau); end for j=1:length(F.b), Fvals = Fvals + F.b(j)*sin(j*tau); end