% quick discussion of polyomial interpolation, showing slightly edited % versions of the book's m-files for it. >> type interpn2 function p = InterpN2(x,y) % p = InterpN2(x,y) % The Newton polynomial interpolant. % x is a column n-vector with distinct components and y is % a column n-vector. c = p.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 p.c = y; p.x = x; >> type hornern function pVal = HornerN(p,z) % pVal = HornerN(p,z) % Evaluates the Newton interpolant on z where % c = p.c and x = p.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. pVal = repmat(p.c(end),size(z)); for k=length(p.c)-1:-1:1 pVal = (z-p.x(k)).*pVal + p.c(k); end % Note that I made InterpN2 return a structure, with the two fields c and x, % and that I made, correspondingly, HornerN, to have just two input arguments % rather than 3, with the first one the structure p, expected to contain both % the coefficients and the corresponding points x . % (To be honest, in class I only talked about editing these two m-files in this % fashion. Here, I actually made those changes, for the record.) % Here is an example >> x = rand(1,10)*(2*pi); y = sin(x); >> p = interpn2(x,y) p = c: [1x10 double] x: [1x10 double] >> xx = linspace(0,2*pi); >> plot(x,y,'o',xx,hornern(p,xx)) % note that this looks very good, even though it uses a polynomial of degree % <= 9. Compare it to the function being interpolated: >> hold on, plot(xx,sin(xx),'k'), hold off % .. it's the same graphically! % On to piecewise polynomials, i.e., chapter 3. Will do it slightly % differently; the notes for this are on the web already. % matlab command diff % What is diff(x) in case x has only one entry? >> diff(1) ans = Empty matrix: 1-by-0 % Now a version of the book's piecewise linear interpolation example, but % using matlab's commands mkpp and ppval , instead of the book's homemade % ones. >> type myshowpl % script: myshowpl n = 10; x=linspace(0,2*pi,n); y = sin(x); pl = mkpp(x,[diff(y)./diff(x) y(1:n-1)]); plot(x,y,'o') xx = linspace(-2,2*pi+2,100); hold on, plot(xx,ppval(pl,xx)), hold off title('broken line interpolant to sin at linspace(0,2*pi,10)') %% Now run it: >> myshowpl %% discussed at length the use of matlab's ppval ... %% specifically how to vectorize the locating of a given z within a given %% increasing sequence x % stressed the matlab command [yo,I] = sort(y) >> [yo,I] = sort([3 1 2]) yo = 1 2 3 I = 2 3 1 %% ... i.e., the second output tells you how to reorder the input vector in %% order to achieve the ordered output vector % Looked at I in case [ignored, I] = sort([x z]) in case x is increasing %% ran the book's script for comparing static with adaptive piecewise linear % interpolation >> showpwl2 diary off