Spline Toolbox™

Construction of the Chebyshev Spline

This section discusses these aspects of the Chebyshev spline construction:

• What Is a Chebyshev Spline?
• Choice of Spline Space
• Initial Guess
• Remez Iteration

What Is a Chebyshev Spline?

The Chebyshev spline of order for the knot sequence is the unique element of of max-norm 1 that maximally oscillates on the interval and is positive near . This means that there is a unique strictly increasing -sequence so that the function given by , all , has max-norm 1 on . This implies that , and that , all . In fact, , all . This brings up the point that the knot sequence is assumed to make such an inequality possible, i.e., the elements of are assumed to be continuous.

In short, the Chebyshev spline looks just like the Chebyshev polynomial. It performs similar functions. For example, its extreme sites are particularly good sites to interpolate at from since the norm of the resulting projector is about as small as can be; see the toolbox command chbpnt.

In this example, which can be run via the demo "Construction of the Chebyshev Spline", we try to construct for a particular knot sequence .

Choice of Spline Space

We 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))]);
  

Note the quadruple knot at each end. Since k = 4, this makes [0..8] = [breaks(1)..breaks(lp1)] the interval of interest, with n = length(t)-k the dimension of the resulting spline space . The same knot sequence would have been supplied by

• t=augknt(breaks,k);
  

Initial Guess

As our initial guess for the , we use the knot averages

recommended as good interpolation site choices. These are supplied by

• tau=aveknt(t,k); 
  

We plot the resulting first approximation to , i.e., the spline that satisfies , all :

• b = cumprod(repmat(-1,1,n)); b = b*b(end);
  c = spapi(t,tau,b); 
  fnplt(c,'-.') 
  grid
  

Here is the resulting picture.

Figure 9-3: First Approximation to a Chebyshev Spline

Remez Iteration

Starting from this approximation, we use the Remez algorithm to produce a sequence of splines converging to . This means that we construct new as the extrema of our current approximation to and try again. Here is the entire loop.

We find the new interior as the zeros of , i.e., the first derivative of , in several steps. First, we differentiate:

• Dc = fnder(c);
  

Next, we take the zeros of the control polygon of as our first guess for the zeros of . For this, we 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.

We refine this estimate for the zeros of by two steps of the secant method, taking tau and the resulting guess as our first approximations. First, we evaluate 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. We guard against division by zero by setting the function value difference to 1 in case it is zero. Since 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 we take these sites as our new tau,

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

and check the extrema values of our current approximation there:

• extremes = abs(fnval(c, tau));
  

The difference

• max(extremes)-min(extremes)
    ans = 0.6905
  

is an estimate of how far we are from total leveling.

We construct a new spline corresponding to our 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)
  

Here is the resulting picture.

Figure 9-4: A More Nearly Level Spline

If this is not close enough, one simply reiterates the loop. For this example, the next iteration already produces to graphic accuracy.

Provide feedback about this page

A Nonlinear ODE

Approximation by Tensor Product Splines

© 1984-2008- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS