Documentation |
This section discusses these aspects of the Chebyshev spline construction:
The Chebyshev spline C=C_{t}=C_{k,t} of order k for the knot sequence t=(t_{i}: i=1:n+k) is the unique element of S_{k,t} of max-norm 1 that maximally oscillates on the interval [t_{k}..t_{n+1}] and is positive near t_{n+1}. This means that there is a unique strictly increasing n-sequence τ so that the function C=C_{t}∊S_{k,t} given by C(τ_{i})=(–1)^{n – 1}, all i, has max-norm 1 on [t_{k}..t_{n+1}]. This implies that τ_{1}=t_{k},τ_{n}=t_{n+1}, and that t_{i} < τ_{i }< t_{k+i}, for all i. In fact, t_{i+1} ≤ τ_{i }≤ t_{i+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 S_{k,t} are assumed to be continuous.
In short, the Chebyshev spline C looks just like the Chebyshev polynomial. It performs similar functions. For example, its extreme sites τ are particularly good sites to interpolate at from S_{k,t} because the norm of the resulting projector is about as small as can be; see the toolbox command chbpnt.
You can run the example Construction of a Chebyshev Spline to construct C for a particular knot sequence t.
You deal with cubic splines, i.e., with order
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] = [breaks(1)..breaks(lp1)] the interval [t_{k}..t_{n+1}] of interest, with n = length(t)-k the dimension of the resulting spline space S_{k,t}. The same knot sequence would have been supplied by
t=augknt(breaks,k);
As the initial guess for the τ, use the knot averages
$${t}_{i}=({t}_{i+1}+\mathrm{...}+{t}_{i+k-1})/(k-1)$$
recommended as good interpolation site choices. These are supplied by
tau=aveknt(t,k);
Plot the resulting first approximation to C, i.e., the spline c that satisfies c(τ_{i})=(–1)^{n-–i}, all i:
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.
You find the new interior τ_{i} as the zeros of Dc, i.e., the first derivative of c, in several steps. First, differentiate:
Dc = fnder(c);
Next, take the zeros of the control polygon of Dc as your first guess for the zeros of Dc. For this, you must take apart the spline Dc.
[knots,coefs,np,kp] = fnbrk(Dc,'knots','coefs','n','order');
The control polygon has the vertices (tstar(i),coefs(i)), with tstar the knot averages for the spline, provided by aveknt:
tstar = aveknt(knots,kp);
Here are the zeros of the resulting control polygon of Dc:
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.
Refine this estimate for the zeros of Dc by two steps of the secant method, taking tau and the resulting guess as your first approximations. First, evaluate Dc at both sets:
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
The check
max(abs(values.')) ans = 4.1176 5.7789 0.4644 0.1178
shows the improvement.
Now take these sites as your new tau,
tau = [tau(1) sites(4,:) tau(n)];
and check the extrema values of your current approximation there:
extremes = abs(fnval(c, tau));
The difference
max(extremes)-min(extremes) ans = 0.6905
is an estimate of how far you are from total leveling.
Construct a new spline corresponding to the new choice of tau and plot it on top of the old:
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.