www.gusucode.com > curvefit 案例源码程序 matlab代码 > curvefit/tspdem.m
%% Bivariate Tensor Product Splines % % This example shows how to use the spline commands in Curve Fitting Toolbox(TM) % to fit tensor product splines to bivariate gridded data. % Copyright 1987-2012 The MathWorks, Inc. %% Introduction % Since Curve Fitting Toolbox can handle splines with _vector_ 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. % % 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 construction not provided by the commands in the toolbox. %% Example: Least-Squares Approximation to Gridded Data % Consider, for example, least-squares approximation to given data % % z(i,j) = f(x(i),y(j)) for i = 1:I, j = 1:J. % % Here are some gridded data, taken from Franke's sample function. Note that % the grid is somewhat denser near the boundary, to help pin down the % approximation there. x = sort([(0:10)/10,.03 .07, .93 .97]); y = sort([(0:6)/6,.03 .07, .93 .97]); [xx,yy] = ndgrid(x,y); % note: ndgrid rather than meshgrid z = franke(xx,yy); mesh(x,y,z.'); xlabel('x'); ylabel('y'); view(150,50); title('Data from the Franke Function'); %% A note about NDGRID vs. MESHGRID % Note that the statements % % [xx,yy] = ndgrid(x,y); % z = franke(xx,yy); % % used above make certain that |z(i,j)| is the value of the function being % approximated at the grid point |(x(i),y(j))|. % % However, the MATLAB(R) command |mesh(x,y,z)| expects |z(j,i)| (note the reversed % order of |i| and |j|) as the value at the grid point |(x(i),y(j))|. For that % reason, the above plot was generated by the statement % % mesh(x,y,z.'); % % i.e., using the transpose of the matrix |z|. % % Such transposing would not have been necessary had we used |meshgrid| % instead of |ndgrid|. But the resulting |z| would not have followed % approximation theory standards. %% Choice of Spline Space in the Y-Direction % Next, we choose a spline order |ky| and a knot sequence |knotsy| for the % y-direction ky = 3; knotsy = augknt([0,.25,.5,.75,1],ky); %% % and then obtain sp = spap2(knotsy,ky,y,z); %% % a spline curve whose |i|-th component is an approximation to the curve % |(y,z(i,:))| for |i=1:I|. %% Evaluation % In particular, yy = -.1:.05:1.1; vals = fnval(sp,yy); %% % creates the matrix |vals| whose |(i,j)|-th element 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|. mesh(x,yy,vals.'); xlabel('x'); ylabel('y'); view(150,50); title('Simultaneous Approximation to All Curves in the Y-Direction'); %% % 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 |sp|. % % Also note the "ridges" that run along the y-direction, most noticeable near % the peaks of the surface. They confirm that we are plotting smooth curves in % one direction only. %% From Curves to a Surface; Choosing a Spline Space in the X-Direction % To get an actual surface, we now have to go one step further. Consider the % coefficients |coefsy| of the spline |sp|, as obtained by coefsy = fnbrk(sp,'c'); %% % Abstractly, you can think of the spline |sp| as the vector-valued function % % y |--> sum coefsy(:,r) B_{r,ky}(y) % r % % with the |i|-th element, |coefsy(i,r)|, of the vector coefficient |coefsy(:,r)| % corresponding to |x(i)| for |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|. kx = 4; knotsx = augknt(0:.2:1,kx); sp2 = spap2(knotsx,kx,x,coefsy.'); %% % The use of the |spap2| command here needs, perhaps, an explanation. % % Recall that |spap2(knots,k,x,fx)| treats |fx(:,j)| as the value at |x(j)|, % i.e., takes each _column_ of |fx| as a data value. Since we wanted to fit % the value |coefsy(i,:)| at |x(i)|, for all |i|, we have to provide |spap2| % with the _transpose_ of |coefsy|. %% % Now consider the _transpose_ of the coefficient matrix of the resulting spline % "curve" |sp2|, obtained as coefs = fnbrk(sp2,'c').'; %% % |coefs| provides the _bivariate_ spline approximation % % (x,y) |--> sum sum coefs(q,r) B_{q,kx}(x) B_{r,ky}(y) % q r % % to the original data % % (x(i),y(j)) |--> f(x(i),y(j)) = z(i,j). % % 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 plot the values. xv = 0:.025:1; yv = 0:.025:1; values = spcol(knotsx,kx,xv)*coefs*spcol(knotsy,ky,yv).'; mesh(xv,yv,values.'); xlabel('x'); ylabel('y'); view(150,50); title('The Spline Approximant'); %% Why Does This Evaluation Work? % The statement % % values = spcol(knotsx,kx,xv) * coefs * spcol(knotsy,ky,yv).' % % used above makes good sense since, for example, |spcol(knotsx,kx,xv)| is the matrix % whose |(i,q)|-th 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 % % sum sum coefs(q,r) B_{q,kx}(x) B_{r,ky}(y) % q r % = sum sum B_{q,kx}(x) coefs(q,r) B_{r,ky}(y) % q r %% % at |(x,y) = (xv(i),yv(j))|. %% More Efficient Alternatives % Since the matrices |spcol(knotsx,kx,xv)| and |spcol(knotsy,ky,yv)| are % banded, it may be more efficient for "large" |xv| and |yv| (though perhaps % more memory-consuming) to make use of |fnval|. value2 = fnval(spmak(knotsx,fnval(spmak(knotsy,coefs),yv).'),xv).'; %% % In fact, |fnval| and |spmak| can deal directly with multivariate splines, % hence the above statement can be replaced by value3 = fnval( spmak({knotsx,knotsy},coefs), {xv,yv} ); %% % Better yet, the construction of the approximation can be done by _one_ % call to |spap2|, therefore we can obtain these values directly from the % given data by the statement value4 = fnval( spap2({knotsx,knotsy},[kx ky],{x,y},z), {xv,yv} ); %% Check it Out % Here is a check, specifically, the _relative_ difference between the values % computed in these four different ways. diffs = abs(values-value2) + abs(values-value3) + abs(values-value4); max(max(diffs)) / max(max(abs(values))) %% % The four methods return the same values, up to round-off error. %% Error of the Approximation % Here is a plot of the error, i.e., the difference between the given data % value and the value of the spline approximation at those data sites. errors = z - spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y).'; mesh(x,y,errors.'); xlabel('x'); ylabel('y'); view(150,50); title('Error at the Given Data Sites'); %% % The _relative_ error is max(max(abs(errors))) / max(max(abs(z))) %% % This is perhaps not too impressive. On the other hand, the ratio % % (degrees of freedom used) / (number of data points) % % is only numel(coefs)/numel(z) %% Apparent Bias of This Approach % The approach followed here seems |biased|: 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|. % % 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|? % % Perhaps surprisingly, the final approximation is the same, up to roundoff. % The next section contains the numerical experiment confirming that. %% Doing It the Other Way Around: Start With a Spline Space in X % First, we fit a spline curve to the data, but this time with |x| as the % independent variable, hence it is the _rows_ of |z| which now become the % data values. Correspondingly, we must supply |z.'| (rather than |z|) to % |spap2|, and obtain spb = spap2(knotsx,kx,x,z.'); %% % a spline approximation to all the curves |(x,z(:,j))| for |j=1:J|. % In particular, valsb = fnval(spb,xv).'; %% % creates a matrix whose |(i,j)|-th element 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|. mesh(xv,y,valsb.'); xlabel('x'); ylabel('y'); view(150,50); title('Simultaneous Approximation to All Curves in the X-Direction'); %% % Again note the ridges, this time running along the x-direction. They confirm % that, once again, we are plotting smooth curves in one direction only. %% From Curves to a Surface: Using a Spline Space in the Y-Direction % Now comes the second step, to get the actual surface. % % Let |coefsx| be the coefficients for |spb|, i.e., coefsx = fnbrk(spb,'c'); %% % Abstractly, you can think of the spline |spb| as the vector-valued function % % x |--> sum coefsx(r,:) B_{r,kx}(x) % r % % with the |j|-th entry |coefsx(r,j)| of the vector coefficient |coefsx(r,:)| % corresponding to |y(j)|, for 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|. spb2 = spap2(knotsy,ky,y,coefsx.'); %% % In the construction of |spb2|, we again need to transpose the coefficient % matrix from |spb|, since |spap2| takes the columns of its last input % argument as the data values. % % For this reason, there is now no need to transpose the coefficient matrix % |coefsb| of the resulting "curve". coefsb = fnbrk(spb2,'c'); %% % Claim: |coefsb| equals the earlier coefficient array |coefs|, up to round-off. % For a proof of this, see the discussion of the tensor product construct in % Curve Fitting Toolbox documentation. Here, we simply make the following check. max(max(abs(coefs - coefsb))) %% % Thus, the _bivariate_ spline approximation % % (x,y) |--> sum sum coefsb(q,r) B_{q,kx}(x) B_{r,ky}(y) % q r % % to the original data % % (x(i),y(j)) |--> f(x(i),y(j)) = z(i,j) % % obtained coincides with the earlier one, which generated |coefs| rather than % |coefsb|. %% % As observed earlier, you can carry out the entire construction we % just went through (in two ways), using just two statements, one for the % construction of the least-squares approximant, the other for its evaluation % at a rectangular mesh. tsp = spap2({knotsx,knotsy},[kx,ky],{x,y},z); valuet = fnval(tsp,{xv,yv}); %% % Here, as another check, is the relative difference between the values % computed earlier and those computed now: max(max(abs(values-valuet))) / max(max(abs(values))) %% Another Example: Interpolation % Since the data come from a smooth function, we should be interpolating % it, i.e., using |spapi| instead of |spap2|, or, equivalently, use |spap2| % with the appropriate knot sequences. For illustration, here is the same % process done with |spapi|. % % To recall, the (univariate) data sites were x %% y %% % We use again quadratic splines in |y|, hence use knots midway between data % sites. knotsy = augknt( [0 1 (y(2:(end-2))+y(3:(end-1)))/2 ], ky); spi = spapi(knotsy,y,z); coefsy = fnbrk(spi,'c'); %% Interpolation of Resulting Coefficients % We use again cubic splines in |x|, and use the not-a-knot condition. We % therefore use all but the second and the second-to-last data points as knots. knotsx = augknt(x([1,3:(end-2),end]), kx); spi2 = spapi(knotsx,x,coefsy.'); icoefs = fnbrk(spi2,'c').'; %% Evaluation % We are now ready to evaluate the interpolant ivalues = spcol(knotsx,kx,xv)*icoefs*spcol(knotsy,ky,yv).'; %% % and plot the interpolant at a fine mesh. mesh(xv,yv,ivalues.'); xlabel('x'); ylabel('y'); view(150,50); title('The Spline Interpolant'); %% % Again, the steps above can be carried out using just two statements, one for % the construction of the interpolant, the other for its evaluation at a % rectangular mesh. tsp = spapi({knotsx,knotsy},{x,y},z); valuet = fnval(tsp,{xv,yv}); %% % For a check, we also compute the relative difference between the values % computed earlier and those computed now. max(max(abs(ivalues-valuet))) / max(max(abs(ivalues))) %% Error of the Approximation % Next, we compute the error of the interpolant as an approximation to the % Franke function. fvalues = franke(repmat(xv.',1,length(yv)),repmat(yv,length(xv),1)); error = fvalues - ivalues; mesh(xv,yv,error.'); xlabel('x'); ylabel('y'); view(150,50); title('Interpolation error'); %% % The _relative_ approximation error is max(max(abs(error))) / max(max(abs(fvalues))) displayEndOfDemoMessage(mfilename)