Construction of the Chebyshev Spline
This section discusses these aspects of the Chebyshev spline construction:
What Is a Chebyshev Spline?
The Chebyshev spline C=Ct=Ck,t of order k for the knot sequence t=(ti: i=1:n+k) is the unique element of Sk,t of max-norm 1 that maximally oscillates on the interval [tk..tn+1] and is positive near tn+1. This means that there is a unique strictly increasing n-sequence τ so that the function C=Ct∊Sk,t given by C(τi)=(–1)n – 1, all i, has max-norm 1 on [tk..tn+1]. This implies that τ1=tk,τn=tn+1, and that ti < τi < tk+i, for all i. In fact, ti+1 ≤ τi ≤ ti+k–1, all i. This brings up the point that the knot sequence is assumed to make such an inequality possible, i.e., the elements of Sk,t are assumed to be continuous.
In short, the Chebyshev spline C looks just like
polynomial. It performs similar functions. For example, its extreme
sites τ are particularly good
sites to interpolate at from Sk,t because
the norm of the resulting projector is about as small as can be; see
the toolbox command
You can run the example Construction of a Chebyshev Spline to construct C for a particular knot sequence t.
Choice of Spline Space
k = 4;
and use the break sequence
breaks = [0 1 1.1 3 5 5.5 7 7.1 7.2 8]; lp1 = length(breaks);
and use simple interior knots, i.e., use the knot sequence
t = breaks([ones(1,k) 2:(lp1-1) lp1(:,ones(1,k))]); n = length(t)-k;
Note the quadruple knot at each end. Because
k = 4,
this makes [0..8] = [
the interval [tk..tn+1]
of interest, with
n = length(t)-
dimension of the resulting spline space Sk,t.
The same knot sequence would have been supplied by
recommended as good interpolation site choices. These are supplied by
b = cumprod(repmat(-1,1,n)); b = b*b(end); c = spapi(t,tau,b); fnplt(c,'-.') grid
Here is the resulting plot.
First Approximation to a Chebyshev Spline
Starting from this approximation, you use the Remez algorithm to produce a sequence of splines converging to C. This means that you construct new τ as the extrema of your current approximation c to C and try again. Here is the entire loop.
Dc = fnder(c);
[knots,coefs,np,kp] = fnbrk(Dc,'knots','coefs','n','order');
tstar = aveknt(knots,kp);
Here are the zeros of the resulting control polygon of
npp = (1:np-1); guess = tstar(npp) -coefs(npp).*(diff(tstar)./diff(coefs));
This provides already a very good first guess for the actual zeros.
sites = tau(ones(4,1),2:n-1); sites(1,:) = guess; values = zeros(4,n-2); values(1:2,:) = reshape(fnval(Dc,sites(1:2,:)),2,n-2);
Now come two steps of the secant method. You guard against division by zero by setting the function value difference to 1 in case it is zero. Because Dc is strictly monotone near the sites sought, this is harmless:
for j=2:3 rows = [j,j-1];Dcd=diff(values(rows,:)); Dcd(find(Dcd==0)) = 1; sites(j+1,:) = sites(j,:) ... -values(j,:).*(diff(sites(rows,:))./Dcd); values(j+1,:) = fnval(Dc,sites(j+1,:)); end
max(abs(values.')) ans = 4.1176 5.7789 0.4644 0.1178
shows the improvement.
Now take these sites as your new
tau = [tau(1) sites(4,:) tau(n)];
and check the extrema values of your current approximation there:
extremes = abs(fnval(c, tau));
max(extremes)-min(extremes) ans = 0.6905
is an estimate of how far you are from total leveling.
c = spapi(t,tau,b); sites = sort([tau (0:100)*(t(n+1)-t(k))/100]); values = fnval(c,sites); hold on, plot(sites,values)
The following code turns on the grid and plots the locations of the extrema.
grid on plot( tau(2:end-1), zeros( 1, np-1 ), 'o' ) hold off legend( 'Initial Guess', 'Current Guess', 'Extreme Locations',... 'location', 'NorthEastOutside' );
Following is the resulting figure (legend not shown).
A More Nearly Level Spline
If this is not close enough, one simply reiterates the loop. For this example, the next iteration already produces C to graphic accuracy.