📄 chebdem.html
字号:
<html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Example: Construction of a Chebyshev Spline</title> <meta name="generator" content="MATLAB 7.1"> <meta name="date" content="2005-07-27"> <meta name="m-file" content="chebdem"> <link rel="stylesheet" type="text/css" href="../../matlab/demos/private/style.css"> </head> <body> <div class="header"> <div class="left"><a href="matlab:edit chebdem">Open chebdem.m in the Editor</a></div> <div class="right"><a href="matlab:echodemo chebdem">Run in the Command Window</a></div> </div> <div class="content"> <h1>Example: Construction of a Chebyshev Spline</h1> <introduction> <p>Illustration of toolbox use in a nontrivial problem.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">Chebyshev (aka equioscillating) spline defined</a></li> <li><a href="#3">Choice of spline space</a></li> <li><a href="#5">Initial guess</a></li> <li><a href="#6">Iteration</a></li> <li><a href="#11">End of first iteration step</a></li> <li><a href="#13">Use of Chebyshev-Demko points</a></li> </ul> </div> <h2>Chebyshev (aka equioscillating) spline defined<a name="1"></a></h2> <p>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))=(-)^{n-i} , 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} , all i . In fact, t_{i+1} <= TAU(i) <= t_{i+k-1} , all i . This brings up the point that the knot sequence t is assumed to make such an inequality possible, i.e., the elements of S_{k,t} are assumed to be continuous. </p><pre class="codeinput">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),<span class="string">'linew'</span>,2)hold <span class="string">on</span>breaks = knt2brk(t); bbb = repmat(breaks,3,1);sss = repmat([1;-1;NaN],1,length(breaks));plot(bbb(:), sss(:),<span class="string">'r'</span>)title(<span class="string">'The Chebyshev spline for a particular knot sequence (|) .'</span>)</pre><img vspace="5" hspace="5" src="chebdem_01.png"> <p>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. </p><pre class="codeinput">plot(tau,zeros(size(tau)),<span class="string">'k+'</span>)xlabel(<span class="string">' Also its extrema (+).'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="chebdem_02.png"> <h2>Choice of spline space<a name="3"></a></h2> <p>In this example, we try to construct C for a given spline space.</p> <p>We deal with cubic splines with simple interior knots, specified by</p><pre class="codeinput">k = 4;breaks = [0 1 1.1 3 5 5.5 7 7.1 7.2 8];t = augknt( breaks, k )</pre><pre class="codeoutput">t = Columns 1 through 13 0 0 0 0 1.0000 1.1000 3.0000 5.0000 5.5000 7.0000 7.1000 7.2000 8.0000 Columns 14 through 16 8.0000 8.0000 8.0000</pre><p>... thus getting a spline space of dimension</p><pre class="codeinput">n = length(t)-k</pre><pre class="codeoutput">n = 12</pre><h2>Initial guess<a name="5"></a></h2> <p>As our initial guess for the TAU, we use the knot averages</p><pre> TAU(i) = (t_{i+1} + ... + t_{i+k-1})/(k-1)</pre><p>recommended as good interpolation site choices, and plot the resulting first approximation to C :</p><pre class="codeinput">tau = aveknt(t,k)b = (-ones(1,n)).^(n-1:-1:0);c = spapi(t,tau,b);plot(breaks([1 end]),[1 1],<span class="string">'k'</span>, breaks([1 end]),[-1 -1],<span class="string">'k'</span>), hold <span class="string">on</span>grid <span class="string">off</span>, lw = 1; fnplt(c,<span class="string">'r'</span>,lw)title(<span class="string">'First approximation to an equioscillating spline'</span>)</pre><pre class="codeoutput">tau = 0 0.3333 0.7000 1.7000 3.0333 4.5000 5.8333 6.5333 7.1000 7.4333 7.7333 8.0000</pre><img vspace="5" hspace="5" src="chebdem_03.png"> <h2>Iteration<a name="6"></a></h2> <p>For the complete levelling, 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. </p> <p>Finding these extrema is itself an iterative process, namely for finding the zeros of the derivative Dc of our current approximation c . </p> <p>We take the zeros of the control polygon of Dc as our first guess for the zeros of Dc .</p> <p>The control polygon has the vertices (TSTAR(i),COEFS(i)) , with TSTAR the knot averages for the spline, as supplied by AVEKNT, and COEFS supplied by SPBRK(Dc). </p><pre class="codeinput">Dc = fnder(c);[knots,coefs,np,kp] = spbrk(Dc);tstar = aveknt(knots,kp);</pre><p>Here are the zeros of the control polygon of Dc :</p><pre class="codeinput">npp = 1:np-1;guess = tstar(npp) - coefs(npp).*(diff(tstar)./diff(coefs));plot(guess,zeros(1,np-1),<span class="string">'o'</span>)xlabel(<span class="string">'...and the zeros of the control polygon of its first derivative'</span>)</pre><img vspace="5" hspace="5" src="chebdem_04.png"> <p>This provides already a very good first guess for the actual extrema of Dc .</p> <p>Now we evaluate Dc at both these sets of sites:</p><pre class="codeinput">sites = repmat( tau(2:n-1), 4,1 );sites(1,:) = guess;values = zeros(4,n-2);values(1:2,:) = fnval(Dc,sites(1:2,:));</pre><p>... and 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 (but guard against division by zero): </p><pre class="codeinput"><span class="keyword">for</span> j = 2:3 rows = [j,j-1]; Dcd = diff(values(rows,:)); Dcd(Dcd==0) = 1; sites(j+1,:) = sites(j,:)-values(j,:).*(diff(sites(rows,:))./Dcd); values(j+1,:) = fnval(Dc,sites(j+1,:));<span class="keyword">end</span></pre><p>We take the last iterate as our new guess for TAU</p><pre class="codeinput">tau = [tau(1) sites(4,:) tau(n)]plot(tau(2:n-1),zeros(1,n-2),<span class="string">'x'</span>)xlabel(<span class="string">' ... and the computed extrema (x) of the current Dc.'</span>)</pre><pre class="codeoutput">tau = 0 0.2759 0.9082 1.7437 3.0779 4.5532 5.5823 6.5843 7.0809 7.3448 7.7899 8.0000</pre><img vspace="5" hspace="5" src="chebdem_05.png"> <h2>End of first iteration step<a name="11"></a></h2> <p>We plot the resulting new approximation c = spapi(t,tau,b) to the Chebyshev spline.</p><pre class="codeinput">c = spapi(t,tau,b);fnplt( c, <span class="string">'k'</span>, lw )xlabel(<span class="string">'... and a more nearly equioscillating spline'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="chebdem_06.png"> <p>If this is not close enough, simply try again, starting from this new TAU. For this particular example, already the next iteration provides the Chebyshev spline to graphic accuracy. </p> <p>The Chebyshev spline for a given spline space S_{k,t} is available as optional output from the command CHBPNT in this toolbox, along with its extrema. These extrema were proposed as good interpolation sites by Steven Demko, hence are now called the Chebyshev-Demko sites. The next slide shows an example of their use. </p> <h2>Use of Chebyshev-Demko points<a name="13"></a></h2> <p>If you have decided to approximate the square-root function on the interval [0 .. 1] by cubic splines with knot sequence</p><pre> k = 4; n = 10; t = augknt(((0:n)/n).^8,k);</pre><p>then a good approximation to the square-root function from that specific spline space is given by</p><pre> tau = chbpnt(t,k); sp = spapi(t,tau,sqrt(tau));</pre><p>as is evidenced by the near equi-oscillation of the error.</p><pre class="codeinput">k = 4; n = 10; t = augknt(((0:n)/n).^8,k);tau = chbpnt(t,k); sp = spapi(t,tau,sqrt(tau));xx = linspace(0,1,301); plot(xx, fnval(sp,xx)-sqrt(xx))title(<span class="string">'The error in interpolant to sqrt at Chebyshev-Demko sites.'</span>)</pre><img vspace="5" hspace="5" src="chebdem_07.png"> <p class="footer">Copyright 1987-2005 C. de Boor and The MathWorks, Inc.<br> Published with MATLAB® 7.1<br></p> </div> <!--##### SOURCE BEGIN #####%% Example: Construction of a Chebyshev Spline%% Illustration of toolbox use in a nontrivial problem.% Copyright 1987-2005 C. de Boor and The MathWorks, Inc.% $Revision: 1.18.4.2 $%% Chebyshev (aka 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))=(-)^{n-i} , 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} , all i . In fact,% t_{i+1} <= TAU(i) <= t_{i+k-1} , all i . This brings up the point% that the knot sequence t is assumed to make such an inequality possible,% i.e., the elements of S_{k,t} are assumed to be 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),'linew',2)hold onbreaks = knt2brk(t); bbb = repmat(breaks,3,1);sss = repmat([1;-1;NaN],1,length(breaks));plot(bbb(:), sss(:),'r')title('The Chebyshev spline for a particular knot sequence (|) .')%%% 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.plot(tau,zeros(size(tau)),'k+')xlabel(' Also its extrema (+).')hold off%% 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 byk = 4;breaks = [0 1 1.1 3 5 5.5 7 7.1 7.2 8];t = augknt( breaks, k )%%% ... thus getting a spline space of dimensionn = length(t)-k%% 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)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 ongrid off, lw = 1; fnplt(c,'r',lw)title('First approximation to an equioscillating spline')%% Iteration% For the complete levelling, 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 .%% We take the zeros of the control polygon of Dc as our first guess for the% zeros of Dc .%% The control polygon has the vertices (TSTAR(i),COEFS(i)) , with% TSTAR the knot averages for the spline, as supplied by AVEKNT, and COEFS% supplied by SPBRK(Dc).Dc = fnder(c);[knots,coefs,np,kp] = spbrk(Dc);tstar = aveknt(knots,kp);%%% Here are the zeros of the control polygon of Dc :npp = 1:np-1;guess = tstar(npp) - coefs(npp).*(diff(tstar)./diff(coefs));plot(guess,zeros(1,np-1),'o')xlabel('...and the zeros of the control polygon of its first derivative')%%% This provides already a very good first guess for the actual extrema of Dc .%% Now we evaluate Dc at both these sets of sites:sites = repmat( tau(2:n-1), 4,1 );sites(1,:) = guess;values = zeros(4,n-2);values(1:2,:) = fnval(Dc,sites(1:2,:));%%% ... and 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 (but guard against division by zero):for j = 2:3 rows = [j,j-1]; Dcd = diff(values(rows,:)); Dcd(Dcd==0) = 1; 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 our new guess for TAUtau = [tau(1) sites(4,:) tau(n)]plot(tau(2:n-1),zeros(1,n-2),'x')xlabel(' ... and the computed extrema (x) of the current Dc.')%% End of first iteration step% We plot the resulting new approximation% c = spapi(t,tau,b)% to the Chebyshev spline.c = spapi(t,tau,b);fnplt( c, 'k', lw )xlabel('... and a more nearly equioscillating spline')hold off%%% If this is not close enough, simply try again, starting from this new TAU.% For this particular example, already the next iteration provides the% Chebyshev spline to graphic accuracy.%% The Chebyshev spline for a given spline space S_{k,t} is available% as optional output from the command CHBPNT in this toolbox, along with its% extrema.% These extrema were proposed as good interpolation sites by Steven Demko,% hence are now called the Chebyshev-Demko sites.% The next slide shows an example of their use.%% Use of Chebyshev-Demko points% 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 equi-oscillation of the error.k = 4; n = 10; t = augknt(((0:n)/n).^8,k);tau = chbpnt(t,k); sp = spapi(t,tau,sqrt(tau));xx = linspace(0,1,301); plot(xx, fnval(sp,xx)-sqrt(xx))title('The error in interpolant to sqrt at Chebyshev-Demko sites.')displayEndOfDemoMessage(mfilename)##### SOURCE END #####--> </body></html>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -