% quick discussion of inverse linear interpolation % adaptive piecewise linear interpolation % here is a version of the book's pwLAdapt that treats points in the plane % as one item, rather than as two. >> type pwladapt function xy = pwLAdapt(fname,xyL,xyR,delta,hmin) % xy = pwLAdapt(fname,xyL,xyR,delta,hmin) % Adaptively determines interpolation points for a piecewise linear % approximation to the function specified. % % fname is a string that specifies a function of the form y = f(u). % xyL =: [xL;f(xL)] and xyR =: [xR;f(xR)] specify the interval [xL .. xR] % and also provide the corresponding function values. % delta is a (positive) uppper bound on the error. % hmin/2 is a positive lower bound on the interval length. % % xy(1,:) =: x is an vector with the property that % xL = x(1) < ... < x(end) = xR, % and xy(2,:) is the corresponding vector of values of f. % % Each subinterval is either <= hmin in length or has the property % that at its midpoint m, |f(m) - L(m)| <= delta, where % L(x) is the line that connects (xR,fR) with (xL,fL) % % if (xyR(1)-xyL(1)) <= hmin % Subinterval is acceptable xy = [xyL, xyR]; else mid = (xyL(1)+xyR(1))/2; fmid = feval(fname,mid); if (abs(((xyL(2)+xyR(2))/2) - fmid) <= delta ) % Subinterval acceptable xy = [xyL, xyR]; else % Produce left and right partitions, then synthesize xyLeft = pwLAdapt(fname,xyL,[mid;fmid],delta,hmin); xyRight = pwLAdapt(fname,[mid;fmid],xyR,delta,hmin); xy = [xyLeft, xyRight(:,2:end)]; end end % Then discussed cubic Hermite interpolation, as cubic interpolation to data % at four points, as the left two coalesce and the right two coalesce. % The book's script ShowHermite gives a nice illustration: >> showhermite % Discussed the construction of the cubic Hermite polynomial interpolant using % divided difference table and the Newton form, and using the definition % f[t_0, ...,t_k] = f^{(k)}(t_0))/k! in case t_0 = t_1 = ... t_k . % Here is the use of the resulting formulas for the construction of the % piecewise cubic Hermite interpolant: >> type pwch function pc = pwch(x,y,s) % % Pre: x = strictly increasing sequence of length n % y, s = sequence of the same orientation and length as x % % Post: pc = a description, ready for use with ppval, of the piecewise % cubic function q that satisfies % q(x_i) = y_i, q'(x_i) = s_i, i=1:n . h = diff(x); dy = diff(y)./h; dzzh = (dy-s(1:end-1))./h; dzhh = (s(2:end)-dy)./h; pc = mkpp(x,[(dzhh-dzzh)./h 2*dzzh-dzhh s(1:end-1) y(1:end-1)]); %% Note how it differs from the book's pwch command, in that it outputs just one %% item, for direct use in matlab's ppval (rather than four vectors, two which %% one then has to remember to adjoin also the sequence x of breakpoints %% when one wants to evaluate the interpolant using the book's evaluator). %% Finally, started on spline interpolation. diary off