📄 tspdem.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>The Tensor Product Construct</title> <meta name="generator" content="MATLAB 7.1"> <meta name="date" content="2005-07-27"> <meta name="m-file" content="tspdem"> <link rel="stylesheet" type="text/css" href="../../matlab/demos/private/style.css"> </head> <body> <div class="header"> <div class="left"><a href="matlab:edit tspdem">Open tspdem.m in the Editor</a></div> <div class="right"><a href="matlab:echodemo tspdem">Run in the Command Window</a></div> </div> <div class="content"> <h1>The Tensor Product Construct</h1> <introduction> <p>Illustrate approximation to data on rectangular grid.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#2">Example: least-squares approximation to gridded data</a></li> <li><a href="#3">NDGRID vs MESHGRID</a></li> <li><a href="#4">Choice of spline space in y-direction</a></li> <li><a href="#5">Evaluation</a></li> <li><a href="#7">From curves to a surface; choosing a spline space in the x-direction</a></li> <li><a href="#10">Why does this evaluation work?</a></li> <li><a href="#11">More efficient alternatives</a></li> <li><a href="#12">Check this out</a></li> <li><a href="#13">Error</a></li> <li><a href="#14">Relative error</a></li> <li><a href="#16">Apparent bias of this approach</a></li> <li><a href="#17">Doing it the other way around</a></li> <li><a href="#19">Comparison with earlier plot of curves</a></li> <li><a href="#20">Approximating to the coefficients</a></li> <li><a href="#26">Interpolation</a></li> <li><a href="#27">Interpolation of resulting coefficients</a></li> <li><a href="#28">Evaluation</a></li> <li><a href="#29">Error</a></li> <li><a href="#30">Relative error</a></li> </ul> </div> <p>Since the toolbox can handle splines with v e c t o r coefficients, it is easy to implement interpolation or approximation to gridded data by tensor product splines. Most spline construction commands in the toolbox take advantage of this. </p> <p>However, you might be interested in seeing a detailed description of how approximation to gridded data by tensor products is actually done for bivariate data. This will also come in handy when you need some tensor product constructon not provided by the M-files in this toolbox. </p> <h2>Example: least-squares approximation to gridded data<a name="2"></a></h2> <p>Consider, for example, least-squares approximation to given data z(i,j) = f(x(i),y(j)), i = 1,...,I, j = 1,...,J . Here are some gridded data, taken from Franke's sample function. Note that I have chosen the grid somewhat denser near the boundary, to help pin down the approximation there. </p><pre class="codeinput">x = sort([(0:10)/10,.03 .07, .93 .97]);y = sort([(0:6)/6,.03 .07, .93 .97]);[xx,yy] = ndgrid(x,y); z = franke(xx,yy);<span class="comment">% NOTE the use of NDGRID in preference to MESHGRID.</span>mesh(x,y,z.'), view(150,50)title(<span class="string">'Data from Franke function'</span>)</pre><img vspace="5" hspace="5" src="tspdem_01.png"> <h2>NDGRID vs MESHGRID<a name="3"></a></h2> <p>Note that the statement [xx,yy] = ndgrid(x,y); z = franke(xx,yy); used here makes certain that Z(i,j) is the value (of the function being approximated) at the grid point (X(i),Y(j)). However, the MATLAB command MESH(X,Y,Z) takes Z(j,i) as the value at the grid point (X(i),Y(j)). For that reason, the above plot is generated by mesh(x,y,z.') i.e., using the transpose of the matrix Z . Such transposing would not have been necessary, had I used MESHGRID instead of NDGRID. But the resulting Z would not have followed approximation theory standards. </p> <h2>Choice of spline space in y-direction<a name="4"></a></h2> <p>Next, I choose some knot sequence KNOTSY and spline order KY for the y-direction, and then use the command</p><pre> sp = spap2(knotsy,ky,y,z);</pre><p>to obtain, in SP, a spline curve whose i-th component is an approximation to the curve (y,z(i,:)), i=1:I.</p><pre class="codeinput">ky = 3; knotsy = augknt([0,.25,.5,.75,1],ky);sp = spap2(knotsy,ky,y,z);</pre><h2>Evaluation<a name="5"></a></h2> <p>In particular,</p><pre> yy = -.1:.05:1.1; vals = fnval(sp,yy);</pre><p>provides the array VALS whose entry VALS(i,j) can be taken as an approximation to the value f(X(i),YY(j)) of the underlying function f at the grid point (X(i),YY(j)) . This is evident when we plot VALS using MESH: </p><pre class="codeinput">yy = -.1:.05:1.1; vals = fnval(sp,yy);mesh(x,yy,vals.'), view(150,50)title(<span class="string">'Simultaneous approximation to all curves in y-direction'</span>)</pre><img vspace="5" hspace="5" src="tspdem_02.png"> <p>Note that, for each X(i), both the first two and the last two values are zero since both the first two and the last two sites in YY are outside the basic interval for the spline in SP. </p> <p>Note also the ridges. They confirm that we are plotting here smooth curves in one direction only.</p> <h2>From curves to a surface; choosing a spline space in the x-direction<a name="7"></a></h2> <p>To get an actual surface, we now have to go one step further. Consider the coefficients COEFSY of the spline in SP, as obtained by coefsy = fnbrk(sp,'c'); </p> <p>Abstractly, you can think of the spline in SP as the vector-valued function y |--> sum_r COEFSY(:,r) B_{r,KY}(y)</p> <p>with the i-th entry COEFSY(i,r) of the vector coefficient COEFSY(:,r) corresponding to X(i) , i=1:I. This suggests approximating each curve (X, COEFSY(:,r)) by a spline, using the same order KX and the same appropriate knot sequence KNOTSX for every r : </p><pre class="codeinput">coefsy = fnbrk(sp,<span class="string">'c'</span>);kx = 4; knotsx = augknt(0:.2:1,kx);sp2 = spap2(knotsx,kx,x,coefsy.');</pre><p>The command sp2 = spap2(knotsx,kx,x,coefsy.'); used here needs, perhaps, an explanation.</p> <p>Recall that spap2(knots,k,x,fx) treats fx(:,j) as the value at x(j) , i.e., takes each c o l u m n of fx as a data value. Since we wanted to fit the value COEFSY(i,:) at X(i), all i , we have to present SPAP2 with the t r a n s p o s e of COEFSY. </p> <p>Now consider the transpose COEFS = CXY.' of the coefficient matrix CXY of the resulting spline `curve', as obtained by the command coefs = fnbrk(sp2,'c').'; </p> <p>COEFS provides the b i v a r i a t e spline approximation</p><pre> (x,y) |--> sum_q sum_r COEFS(q,r) B_{q,KX}(x) B_{r,KY}(y)</pre><p>to the original data</p><pre> (X(i),Y(j)) |--> f(X(i),Y(j)) = Z(i,j) .</pre><p>We use SPCOL to provide the values B_{q,KX}(XV(i)) and B_{r,KY}(YV(j)) needed to evaluate this spline surface at some grid points (XV(i),YV(j)) and then use MESH to plot the values. </p><pre class="codeinput">coefs = fnbrk(sp2,<span class="string">'c'</span>).';xv = 0:.025:1; yv = 0:.025:1;values = spcol(knotsx,kx,xv)*coefs*spcol(knotsy,ky,yv).';mesh(xv,yv,values.'), view(150,50)title(<span class="string">'The spline approximant'</span>)</pre><img vspace="5" hspace="5" src="tspdem_03.png"> <h2>Why does this evaluation work?<a name="10"></a></h2> <p>The command values = spcol(knotsx,kx,xv)*coefs*spcol(knotsy,ky,yv).' used here makes good sense since, e.g., SPCOL(KNOTSX,KX,XV) is the matrix whose (i,q)-entry equals the value B_{q,KX}(XV(i)) at XV(i) of the q-th B-spline of order KX for the knot sequence KNOTSX, while we want to evaluate the expression </p><pre> sum_q sum_r COEFS(q,r) B_{q,KX}(x) B_{r,KY}(y) =</pre><pre> = sum_q sum_r B_{q,KX}(x) COEFS(q,r) B_{r,KY}(y)</pre><p>at (x,y) = (XV(i),YV(j)).</p> <h2>More efficient alternatives<a name="11"></a></h2> <p>Since the matrices SPCOL(KNOTSX,KX,XV) and SPCOL(KNOTSY,KY,YV) are banded, it may be more efficient (though perhaps more memory-consuming) for `large' XV and YV to make use of FNVAL, as follows: </p><pre> value2 = fnval(spmak(knotsx,fnval(spmak(knotsy,coefs),yv).'),xv).';</pre><p>In fact, FNVAL and SPMAK can deal directly with multivariate splines, hence this statement can be replaced by</p><pre> value3 = fnval( spmak({knotsx,knotsy},coefs), {xv,yv} );</pre><p>Better yet, the construction of the approximation can be done by o n e call to SPAP2, therefore we can obtain these values directly from the given data by the statement </p><pre> value4 = fnval( spap2({knotsx,knotsy},[kx ky],{x,y},z), {xv,yv} );</pre><pre class="codeinput">value2 = fnval(spmak(knotsx,fnval(spmak(knotsy,coefs),yv).'),xv).';value3 = fnval( spmak({knotsx,knotsy},coefs), {xv,yv} );value4 = fnval( spap2({knotsx,knotsy},[kx ky],{x,y},z), {xv,yv} );</pre><h2>Check this out<a name="12"></a></h2> <p>Here is a check, viz. the r e l a t i v e difference between the values computed in these four different ways:</p><pre class="codeinput">disp( <span class="keyword">...</span> max(max(abs(values-value2)+abs(values-value3)+abs(values-value4)))/ <span class="keyword">...</span> max(max(abs(values))) <span class="keyword">...</span> )</pre><pre class="codeoutput"> 1.1206e-015</pre><h2>Error<a name="13"></a></h2> <p>Here is also a plot of the error, i.e., the difference between the given data value and the value of the approximation at those data sites: </p><pre class="codeinput">errors = z - spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y).';mesh(x,y,errors.'), view(150,50)title(<span class="string">'Error at the given data sites'</span>)</pre><img vspace="5" hspace="5" src="tspdem_04.png"> <h2>Relative error<a name="14"></a></h2> <p>The r e l a t i v e error is</p><pre class="codeinput">disp( max(max(abs(errors)))/max(max(abs(z))) )</pre><pre class="codeoutput"> 0.0539</pre><p>This is perhaps not too impressive. On the other hand, the ratio of</p><pre> degrees of freedom used ----------------------- number of data points</pre><p>is only</p><pre class="codeinput">disp( numel(coefs)/numel(z) )</pre><pre class="codeoutput"> 0.2909</pre><h2>Apparent bias of this approach<a name="16"></a></h2> <p>The approach followed here seems b i a s e d : We first think of the given data values Z as describing a vector-valued function of y , and then we treat the matrix formed by the vector coefficients of the approximating curve as describing a vector-valued function of x . </p> <p>What happens when we take things in the opposite order, i.e., think of Z as describing a vector-valued function of x , and then treat the matrix made up from the vector coefficients of the approximating curve as describing a vector-valued function of y ? </p> <p>Perhaps surprisingly, the final approximation is the same (up to roundoff). Here is the numerical experiment.</p> <h2>Doing it the other way around<a name="17"></a></h2> <p>First, we fit a spline curve to the data, but this time with x as the independent variable, hence it is the r o w s of Z which now become the data values. Correspondingly, we must supply Z.' (rather than Z) to SPAP2, and obtain via the statement spb = spap2(knotsx,kx,x,z.'); a spline approximation to all the curves (X,Z(:,j)), j=1:J. In particular, the statement valsb = (fnval(spb,xv)).'; provides the array VALSB, whose entry VALSB(i,j) can be taken as an approximation to the value f(XV(i),Y(j)) of the underlying function f at the grid point (XV(i),Y(j)). This is evident when we plot VALSB using MESH : </p><pre class="codeinput">spb = spap2(knotsx,kx,x,z.');valsb = (fnval(spb,xv)).';mesh(xv,y,valsb.'), view(150,50)title(<span class="string">'Simultaneous approximation to all curves in the x-direction'</span>)</pre><img vspace="5" hspace="5" src="tspdem_05.png"> <p>Note the ridges. They confirm that we are plotting here smooth curves in one direction only.</p> <h2>Comparison with earlier plot of curves<a name="19"></a></h2> <p>Here, for a quick comparison is the earlier mesh of curves, with the smooth curves running in the other direction:</p><pre class="codeinput">mesh(x,yy,fnval(sp,yy).'), view(150,50)title(<span class="string">'Simultaneous approximation to all curves in y-direction'</span>)</pre><img vspace="5" hspace="5" src="tspdem_06.png"> <h2>Approximating to the coefficients<a name="20"></a></h2> <p>Now comes the second step, to get the actual surface. Let COEFSX be the coefficients for SPB, i.e., coefsx = fnbrk(spb,'c');</p> <p>Abstractly, you can think of the spline in SPB as the vector-valued function</p><pre> x |--> sum_r coefsx(r,:) B_{r,kx}(x)</pre><p>with the j-th entry COEFSX(r,j) of the vector coefficient COEFSX(r,:) corresponding to Y(j) , all j . Thus we now fit each curve (Y, COEFSX(r,:)) by a spline, using the same order KY and the same (appropriate) knot sequence KNOTSY for each r : </p><pre class="codeinput">mesh(xv,y,valsb.'), view(150,50)title(<span class="string">'Simultaneous approximation to all curves in the x-direction'</span>)coefsx = fnbrk(spb,<span class="string">'c'</span>);spb2 = spap2(knotsy,ky,y,coefsx.');</pre><img vspace="5" hspace="5" src="tspdem_07.png"> <p>In our construction of spb2 = spap2(knotsy,ky,y,coefsx.') we need again to transpose the coefficient array from SPB, since SPAP2 takes the columns of its last input argument as the data values. </p> <p>For this reason, there is now no need to transpose the coefficient array COEFSB of the resulting `curve':</p><pre class="codeinput">coefsb = fnbrk(spb2,<span class="string">'c'</span>);</pre><p>I claim that COEFSB equals the earlier coefficient array COEFS (up to round-off); see the discussion of the tensor product construct in the spline toolbox Tutorial for a proof of this. Here, I simply try the following test: </p><pre class="codeinput">disp( max(max(abs(coefs - coefsb))) )</pre><pre class="codeoutput"> 1.6653e-015</pre><p>Thus, the b i v a r i a t e spline approximation</p><pre> (x,y) |--> sum_q sum_r coefsb(q,r) B_{q,kx}(x) B_{r,ky}(y)</pre><p>to the original data</p><pre> (X(i),Y(j)) |--> f(X(i),Y(j)) = Z(i,j)</pre><p>obtained coincides with the earlier one (which generated COEFS rather than COEFSB).</p> <p>As already observed earlier, you can carry out the entire construction we just went through (even two ways) using just two commands, one for the construction of the least-squares approximant, the other for its evaluation at a rectangular mesh. </p><pre class="codeinput">tsp = spap2({knotsx,knotsy},[kx,ky],{x,y},z);valuet = fnval(tsp,{xv,yv});</pre><p>Here, for another check, is the relative difference between the values computed earlier and those computed now:</p><pre class="codeinput">disp( max(max(abs(values-valuet)))/max(max(abs(values))) )</pre><pre class="codeoutput"> 3.7353e-016</pre><h2>Interpolation<a name="26"></a></h2> <p>Since the data come from a smooth function, we should be interpolating it, i.e., use SPAPI instead of SPAP2, or, equivalently, use SPAP2 with the appropriate knot sequences. For illustration, here is the same process done with SPAPI. </p>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -