📄 splexmpl.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>Spline Toolbox Examples</title> <meta name="generator" content="MATLAB 7.1"> <meta name="date" content="2005-07-27"> <meta name="m-file" content="splexmpl"> <link rel="stylesheet" type="text/css" href="../../matlab/demos/private/style.css"> </head> <body> <div class="header"> <div class="left"><a href="matlab:edit splexmpl">Open splexmpl.m in the Editor</a></div> <div class="right"><a href="matlab:echodemo splexmpl">Run in the Command Window</a></div> </div> <div class="content"> <h1>Spline Toolbox Examples</h1> <introduction> <p>Here are some simple examples that illustrate the use of the spline toolbox.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">Interpolation</a></li> <li><a href="#2">Estimating the error in the interpolation</a></li> <li><a href="#8">Smoothing</a></li> <li><a href="#14">Least-squares approximation</a></li> <li><a href="#15">Knot selection</a></li> <li><a href="#18">Gridded data</a></li> <li><a href="#20">Curves</a></li> <li><a href="#22">Surfaces</a></li> <li><a href="#25">Scattered data</a></li> </ul> </div> <h2>Interpolation<a name="1"></a></h2> <p>One can construct a cubic spline that matches cosine at the following points:</p> <p>(Note that one can view the interpolating spline by using FNPLT)</p><pre class="codeinput">x = 2*pi*[0 1 .1:.2:.9];y = cos(x);cs = csapi(x,y);fnplt(cs,2);xxyy = [-1 7 -1.2 1.2];axis(xxyy)hold <span class="string">on</span>, plot(x,y,<span class="string">'o'</span>), hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_01.png"> <h2>Estimating the error in the interpolation<a name="2"></a></h2> <p>The cosine is 2*pi-periodic. How well does our interpolant in CS do in that regard?</p> <p>We compute the difference in the first derivative at the two endpoints:</p><pre class="codeinput">diff( fnval( fnder(cs), [0 2*pi] ) )</pre><pre class="codeoutput">ans = -0.1375</pre><p>To enforce periodicity, use CSAPE instead of CSAPI:</p><pre class="codeinput">csp = csape( x, y, <span class="string">'periodic'</span> );hold <span class="string">on</span>, fnplt(csp,<span class="string">'g'</span>), hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_02.png"> <pre>Now, the check gives</pre><pre class="codeinput">diff( fnval( fnder(csp), [0 2*pi] ) )</pre><pre class="codeoutput">ans = -4.5612e-017</pre><p>Even the second derivative now matches at the endpoints:</p><pre class="codeinput">diff( fnval( fnder( csp, 2 ), [0 2*pi] ) )</pre><pre class="codeoutput">ans = -4.4409e-016</pre><p>The piecewise linear interpolant to the same data is available via</p><pre class="codeinput">pl = spapi (2, x, y );</pre><p>Here we add it to the previous plot, in red:</p><pre class="codeinput">hold <span class="string">on</span>, fnplt(pl, <span class="string">'r'</span>, 2 ), hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_03.png"> <h2>Smoothing<a name="8"></a></h2> <p>If the data are noisy, one would approximate rather than interpolate. For example, with</p><pre class="codeinput">x = linspace(0,2*pi,51);noisy_y = cos(x) + .2*(rand(size(x))-.5);plot(x,noisy_y,<span class="string">'x'</span>)</pre><img vspace="5" hspace="5" src="splexmpl_04.png"> <p>... interpolation would give the wiggly interpolant marked in blue</p><pre class="codeinput">hold <span class="string">on</span>fnplt( csapi( x, noisy_y ) )hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_05.png"> <p>... while smoothing with a proper tolerance</p><pre class="codeinput">tol = (.05)^2*(2*pi);hold <span class="string">on</span>fnplt( spaps( x, noisy_y, tol ), <span class="string">'r'</span>, 2 )hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_06.png"> <p>... gives the smoothed approximation, shown here in red.</p> <p>The approximation is much worse near the ends of the interval, and is far from periodic. To enforce periodicity, approximate to periodically extended data, then restrict approximation to the original interval: </p><pre class="codeinput">noisy_y([1 end]) = mean( noisy_y([1 end]) );lx = length(x);lx2 = round(lx/2);range = [lx2:lx 2:lx 2:lx2];sps = spaps([x(lx2:lx)-2*pi x(2:lx) x(2:lx2)+2*pi],noisy_y(range),2*tol);hold <span class="string">on</span>fnplt( sps, [0 2*pi], <span class="string">'k'</span>, 2)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_07.png"> <p>giving the more nearly periodic approximation, shown in black.</p> <h2>Least-squares approximation<a name="14"></a></h2> <p>Alternatively, one could use least-squares approximation to the noisy data by a spline with few degrees of freedom.</p> <p>For example, one might try a cubic spline with just four pieces:</p><pre class="codeinput">spl2 = spap2(4, 4, x, noisy_y);fnplt(spl2,<span class="string">'b'</span>,2);axis(xxyy)hold <span class="string">on</span>plot(x,noisy_y,<span class="string">'x'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_08.png"> <h2>Knot selection<a name="15"></a></h2> <p>When using an SP... command to construct a spline, one usually has to specify a particular spline space. This is done by specifying a k n o t s e q u e n c e and an o r d e r , and this may be a bit of a problem. When doing spline interpolation, to data X , Y by splines of order K , then OPTKNT will supply a good knot sequence, as in the following example: </p><pre class="codeinput">k = 5; <span class="comment">% i.e., we are working with quartic splines</span>x = 2*pi*sort([0 1 rand(1,10)]);y = cos(x);sp = spapi( optknt(x,k), x, y );fnplt(sp,2,<span class="string">'g'</span>);hold <span class="string">on</span>, plot(x,y,<span class="string">'o'</span>), hold <span class="string">off</span>axis([-1 7 -1.1 1.1])</pre><img vspace="5" hspace="5" src="splexmpl_09.png"> <p>When doing least-squares approximation, one can use the current approximation to determine a possibly better knot selection with the aid of NEWKNT. For example, here is an approximation </p><pre class="codeinput">x = linspace(0,10,101);y = exp(x);sp0 = spap2( augknt(0:2:10,4), 4, x, y );plot(x,y-fnval(sp0,x),<span class="string">'r'</span>,<span class="string">'linew'</span>,2)</pre><img vspace="5" hspace="5" src="splexmpl_10.png"> <p>... whose error is plotted above in red, and which isn't all that good, compared to the following approximation with the s a m e number of knots, but better distributed (whose error is plotted in black): </p><pre class="codeinput">sp1 = spap2( newknt(sp0), 4,x,y);hold <span class="string">on</span>plot(x,y-fnval(sp1,x),<span class="string">'k'</span>,<span class="string">'linew'</span>,2)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_11.png"> <h2>Gridded data<a name="18"></a></h2> <p>All the interpolation and approximation commands in the toolbox can also handle gridded data (in any number of variables).</p> <p>For example, here is a bicubic spline interpolant to the Mexican Hat function:</p><pre class="codeinput">x =.0001+(-4:.2:4);y = -3:.2:3;[yy,xx] = meshgrid(y,x);r = pi*sqrt(xx.^2+yy.^2);z = sin(r)./r;bcs = csapi( {x,y}, z );fnplt( bcs )axis([-5 5 -5 5 -.5 1])</pre><img vspace="5" hspace="5" src="splexmpl_12.png"> <p>... and here is the l e a s t - s q u a r e s approximation to noisy values of that function on the same grid:</p><pre class="codeinput">knotsx = augknt(linspace(x(1), x(end), 21), 4);knotsy = augknt(linspace(y(1), y(end), 15), 4);bsp2 = spap2({knotsx,knotsy},[4 4], {x,y},z+.02*(rand(size(z))-.5));fnplt(bsp2)axis([-5 5 -5 5 -.5 1])</pre><img vspace="5" hspace="5" src="splexmpl_13.png"> <h2>Curves<a name="20"></a></h2> <p>Gridded data can be handled easily because the toolbox can deal with v e c t o r - v a l u e d splines. This also makes it easy to work with curves. </p> <p>Here, for example, is an approximation to infinity (obtained by putting a cubic spline curve through the marked points):</p><pre class="codeinput">t = 0:8;xy = [0 0;1 1; 1.7 0;1 -1;0 0; -1 1; -1.7 0; -1 -1; 0 0].';infty = csape( t , xy, <span class="string">'periodic'</span>);fnplt( infty , 2 )axis([-2 2 -1.1 1.1])hold <span class="string">on</span>plot(xy(1,:),xy(2,:),<span class="string">'o'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_14.png"> <p>... and here is the same curve, but with motion in a third dimension:</p><pre class="codeinput">roller = csape( t , [ xy ;0 1/2 1 1/2 0 1/2 1 1/2 0], <span class="string">'periodic'</span>);fnplt( roller , 2, [0 4],<span class="string">'b'</span> )hold <span class="string">on</span>fnplt( roller, 2, [4 8], <span class="string">'r'</span>)plot3(0,0,0,<span class="string">'o'</span>)hold <span class="string">off</span><span class="comment">% I have plotted the two halves of the curve in different colors and have</span><span class="comment">% marked the origin, as an aid to visualizing this two-winged space curve.</span></pre><img vspace="5" hspace="5" src="splexmpl_15.png"> <h2>Surfaces<a name="22"></a></h2><pre>Bivariate tensor-product splines with values in R^3 give surfaces.</pre><p>For example, here is a good approximation to a doughnought:</p><pre class="codeinput">x = 0:4; y=-2:2; R = 4; r = 2; clear <span class="string">v</span>v(3,:,:) = [0 (R-r)/2 0 (r-R)/2 0].'*[1 1 1 1 1];v(2,:,:) = [R (r+R)/2 r (r+R)/2 R].'*[0 1 0 -1 0];v(1,:,:) = [R (r+R)/2 r (r+R)/2 R].'*[1 0 -1 0 1];dough0 = csape({x,y},v,<span class="string">'periodic'</span>);fnplt(dough0)axis <span class="string">equal</span>axis <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_16.png"> <p>And here is a crown of normals to that surface.</p><pre class="codeinput">nx = 43;xy = [ones(1,nx); linspace(2,-2,nx)];points = fnval(dough0,xy)';ders = fnval(fndir(dough0,eye(2)),xy);normals = cross(ders(4:6,:),ders(1:3,:));normals = (normals./repmat(sqrt(sum(normals.*normals)),3,1))';pn = [points;points+normals];hold <span class="string">on</span><span class="keyword">for</span> j=1:nx plot3(pn([j,j+nx],1),pn([j,j+nx],2),pn([j,j+nx],3))<span class="keyword">end</span>hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_17.png"> <p>And here is its projection onto the (x,y)-plane.</p><pre class="codeinput">fnplt(fncmb(dough0, [1 0 0; 0 1 0]))axis <span class="string">equal</span>axis <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_18.png"> <h2>Scattered data<a name="25"></a></h2> <p>It is also possible to interpolate to values given at ungridded data sites in the plane. Consider, for example, the task of mapping the unit square smoothly to the unit disk. We construct the data values </p><pre class="codeinput">n = 64; t = linspace(0,2*pi,n+1); t(end) = [];values = [cos(t); sin(t)];plot(values(1,:),values(2,:),<span class="string">'or'</span>), axis <span class="string">equal</span>, axis <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_19.png"> <p>... and corresponding data sites (marked here as x, and each connected to its associated value by an arrow)</p><pre class="codeinput">sites = values./repmat(max(abs(values)),2,1);hold <span class="string">on</span>, plot(sites(1,:),sites(2,:),<span class="string">'xk'</span>)quiver( sites(1,:),sites(2,:), <span class="keyword">...</span> values(1,:)-sites(1,:), values(2,:)-sites(2,:) )hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_20.png"> <p>... and then use TPAPS to construct a bivariate interpolating vector-valued thin-plate spline</p><pre class="codeinput">st = tpaps(sites, values, 1);</pre><p>... that does indeed map the unit square smoothly (approximately) to the unit disk, as its plot via FNPLT indicates. That plot shows the image of a uniformly spaced square grid under this spline map ST. </p><pre class="codeinput">hold <span class="string">on</span>, fnplt(st), hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="splexmpl_21.png"> <p class="footer">Copyright 1987-2005 The MathWorks, Inc.<br> Published with MATLAB® 7.1<br></p> </div> <!--##### SOURCE BEGIN #####%% Spline Toolbox Examples% Here are some simple examples that illustrate the use of the spline toolbox.% Copyright 1987-2005 The MathWorks, Inc.% $Revision: 1.16.4.2 $ $Date: 2005/06/21 19:46:13 $%% Interpolation% One can construct a cubic spline that matches cosine at the following% points:%% (Note that one can view the interpolating spline by using FNPLT)x = 2*pi*[0 1 .1:.2:.9];y = cos(x);cs = csapi(x,y);fnplt(cs,2);xxyy = [-1 7 -1.2 1.2];axis(xxyy)hold on, plot(x,y,'o'), hold off%% Estimating the error in the interpolation% The cosine is 2*pi-periodic. How well does our interpolant in CS do in that% regard?%% We compute the difference in the first derivative at the two endpoints:diff( fnval( fnder(cs), [0 2*pi] ) )%%% To enforce periodicity, use CSAPE instead of CSAPI:csp = csape( x, y, 'periodic' );hold on, fnplt(csp,'g'), hold off%%% Now, the check gives⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -