Note: This page has been translated by MathWorks. Click here to see

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

To view all translated materials including this page, select Country from the country navigator on the bottom of this page.

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.'});