# Construction of a Chebyshev Spline

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

### Chebyshev (a.k.a. Equioscillating) Spline Defined

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'});``` ### Choice of Spline Space

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

### Initial Guess

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');``` ### Iteration

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]);``` ### End of First Iteration Step

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.

### Use of Chebyshev-Demko Points

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