This example shows how to use commands from Curve Fitting Toolbox™ to construct a Chebyshev spline.

By definition, for given knot sequence `t`

of length `n+k`

, `C = C_{t,k}`

is the unique element of `S_{t,k}`

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 `tau`

of length `n`

so that the function `C`

in `S_{k,t}`

given by

` C(tau(i)) = (-1)^{n-i},`

for all `i`

, has max-norm 1 on `[t_k .. t_{n+1}]`

. This implies that

` tau(1) = t_k,`

` tau(n) = t_{n+1},`

and that

` t_i < tau(i) < t_{k+i},`

for all `i`

. In fact,

` t_{i+1} <= tau(i) <= t_{i+k-1},`

for all `i`

. This brings up the point that the knot sequence `t`

is assumed to make such an inequality possible, which turns out to be equivalent to having all the elements of `S_{k,t}`

continuous.

t = augknt([0 1 1.1 3 5 5.5 7 7.1 7.2 8], 4 ); [tau,C] = chbpnt(t,4); xx = sort([linspace(0,8,201),tau]); plot(xx,fnval(C,xx),'LineWidth',2); hold on breaks = knt2brk(t); bbb = repmat(breaks,3,1); sss = repmat([1;-1;NaN],1,length(breaks)); plot(bbb(:), sss(:),'r'); hold off ylim([-2 2]); title('The Chebyshev Spline for a Particular Knot Sequence'); legend({'Chebyshev Spline' 'Knots'});

In short, the Chebyshev spline `C`

looks just like the Chebyshev polynomial. It performs similar functions. For example, its extrema `tau`

are particularly good sites to interpolate at from `S_{k,t}`

since the norm of the resulting projector is about as small as can be.

hold on plot(tau,zeros(size(tau)),'k+'); hold off legend({'Chebyshev Spline' 'Knots' 'Extrema'});

In this example, we try to construct `C`

for a given spline space.

We deal with cubic splines with simple interior knots, specified by

k = 4; breaks = [0 1 1.1 3 5 5.5 7 7.1 7.2 8]; t = augknt(breaks, k)

`t = `*1×16*
0 0 0 0 1.0000 1.1000 3.0000 5.0000 5.5000 7.0000 7.1000 7.2000 8.0000 8.0000 8.0000 8.0000

thus getting a spline space of dimension

n = length(t)-k

n = 12

As our initial guess for the `tau`

, we use the knot averages

` tau(i) = (t_{i+1} + ... + t_{i+k-1})/(k-1)`

recommended as good interpolation site choices, and plot the resulting first approximation to `C`

.

tau = aveknt(t,k)

`tau = `*1×12*
0 0.3333 0.7000 1.7000 3.0333 4.5000 5.8333 6.5333 7.1000 7.4333 7.7333 8.0000

b = (-ones(1,n)).^(n-1:-1:0); c = spapi(t,tau,b); plot(breaks([1 end]),[1 1],'k', breaks([1 end]),[-1 -1],'k'); hold on fnplt(c,'r',1); hold off ylim([-2 2]); title('First Approximation to an Equioscillating Spline');

For the complete leveling, we use the Remez algorithm. This means that we construct a new `tau`

as the extrema of our current approximation, `c`

, to `C`

and try again.

Finding these extrema is itself an iterative process, namely, for finding the zeros of the derivative `Dc`

of our current approximation `c`

.

Dc = fnder(c);

We take the zeros of the control polygon of `Dc`

as our first guess for the zeros of `Dc`

. This control polygon has the vertices `(tstar(i),coefs(i))`

, where `coefs`

are the coefficients of `Dc`

, and `tstar`

the knot averages.

[knots,coefs,np,kp] = fnbrk(Dc, 'knots', 'coefs', 'n', 'order'); tstar = aveknt(knots,kp);

Since the control polygon is piecewise linear, its zeros are easy to compute. Here are those zeros.

npp = 1:np-1; guess = tstar(npp) - coefs(npp).*(diff(tstar)./diff(coefs)); fnplt(Dc,'r'); hold on plot(tstar,coefs,'k.:'); plot(guess,zeros(1,np-1),'o'); hold off title('First Derivative of the Approximation'); legend({'Dc' 'Control Polygon' 'Zeros of Control Polygon'});

This provides a very good first guess for the actual zeros of `Dc`

.

Now we evaluate `Dc`

at both these sets of sites.

sites = [guess; tau(2:n-1)]; values = fnval(Dc,sites);

Then we use two steps of the secant method, getting iterates `sites(3,:)`

and `sites(4,:)`

, with `values(3,:)`

and `values(4,:)`

the corresponding values of `Dc`

.

sites(3:4,:) = 0; values(3:4,:) = 0; for j = 2:3 rows = [j,j-1]; Dcd = diff(values(rows,:)); Dcd(Dcd==0) = 1; % guard against division by zero sites(j+1,:) = sites(j,:)-values(j,:).*(diff(sites(rows,:))./Dcd); values(j+1,:) = fnval(Dc,sites(j+1,:)); end

We take the last iterate as the computed zeros of `Dc`

, i.e., the extrema of our current approximation, `c`

. This is our new guess for `tau`

.

tau = [tau(1) sites(4,:) tau(n)]

`tau = `*1×12*
0 0.2759 0.9082 1.7437 3.0779 4.5532 5.5823 6.5843 7.0809 7.3448 7.7899 8.0000

plot(breaks([1 end]),[1 1],'k', breaks([1 end]),[-1 -1],'k'); hold on fnplt(c,'r',1); plot(guess,zeros(1,np-1),'o'); plot(tau(2:n-1),zeros(1,n-2),'x'); hold off title('First Approximation to an Equioscillating Spline'); ax = gca; h = ax.Children; legend(h([3 1 2]),{'Approximation' 'Extrema' ... 'Zeros of First Derivative''s Control Polygon'}); axis([0 8 -2 2]);

We compute the resulting new approximation to the Chebyshev spline using the new guess for `tau`

.

cnew = spapi(t,tau,b);

The new approximation is more nearly an equioscillating spline.

plot(breaks([1 end]),[1 1],'k', breaks([1 end]),[-1 -1],'k'); hold on fnplt(c,'r',1); fnplt(cnew, 'k', 1); hold off ax = gca; h = ax.Children; legend(h([2 1]),{'First Approximation' 'Updated Approximation'}); axis([0 8 -2 2]);

If this is not close enough, simply try again, starting from this new `tau`

. For this particular example, the next iteration already provides the Chebyshev spline to graphic accuracy.

The Chebyshev spline for a given spline space `S_{k,t}`

, along with its extrema, are available as optional outputs from the `chbpnt`

command in the toolbox. These extrema were proposed as good interpolation sites by Steven Demko, hence are now called the Chebyshev-Demko sites. This section shows an example of their use.

If you have decided to approximate the square-root function on the interval `[0 .. 1]`

by cubic splines with knot sequence

k = 4; n = 10; t = augknt(((0:n)/n).^8,k);

then a good approximation to the square-root function from that specific spline space is given by

tau = chbpnt(t,k); sp = spapi(t,tau,sqrt(tau));

as is evidenced by the near equioscillation of the error.

xx = linspace(0,1,301); plot(xx, fnval(sp,xx)-sqrt(xx)); title({'Error in Interpolant to Square Root','at Chebyshev-Demko Sites.'});