tspdem.html

演示matlab曲线拟和与插直的基本方法

%% Check this out
% Here is a check, viz. the  r e l a t i v e  difference between the values
% computed in these four different ways:

disp( ...
     max(max(abs(values-value2)+abs(values-value3)+abs(values-value4)))/ ...
     max(max(abs(values))) ...
    )

%% Error
% 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:

errors = z -  spcol(knotsx,kx,x)*coefs*spcol(knotsy,ky,y).';
mesh(x,y,errors.'), view(150,50)
title('Error at the given data sites')

%% Relative error
% The  r e l a t i v e  error is
disp( max(max(abs(errors)))/max(max(abs(z))) )

%%
% This is perhaps not too impressive. On the other hand, the ratio of
%
%                degrees of freedom used
%                REPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASHREPLACE_WITH_DASH_DASH-
%                 number of data points
%
% is only
disp( numel(coefs)/numel(z) )

%% Apparent bias of this approach
% 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 .
%
% 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).
% Here is the numerical experiment.

%% Doing it the other way around
% 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 :

spb = spap2(knotsx,kx,x,z.');
valsb = (fnval(spb,xv)).';
mesh(xv,y,valsb.'), view(150,50)
title('Simultaneous approximation to all curves in the x-direction')

%%
% Note the ridges. They confirm that we are plotting here smooth curves in
% one direction only.
%
%% Comparison with earlier plot of curves
% Here, for a quick comparison is the earlier mesh of curves, with the smooth
% curves running in the other direction:

mesh(x,yy,fnval(sp,yy).'), view(150,50)
title('Simultaneous approximation to all curves in y-direction')

%% Approximating to the coefficients
% 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 in SPB as the vector-valued function
%
%       x |REPLACE_WITH_DASH_DASH> sum_r coefsx(r,:) B_{r,kx}(x)
%
% 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 :

mesh(xv,y,valsb.'), view(150,50)
title('Simultaneous approximation to all curves in the x-direction')
coefsx = fnbrk(spb,'c');
spb2 = spap2(knotsy,ky,y,coefsx.');

%%
% 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.
%
% For this reason, there is now no need to transpose the coefficient array
% COEFSB of the resulting `curve':

coefsb = fnbrk(spb2,'c');

%%
% 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:

disp( max(max(abs(coefs - coefsb))) )

%%
% Thus, the  b i v a r i a t e  spline approximation
%
%     (x,y) |REPLACE_WITH_DASH_DASH> sum_q sum_r coefsb(q,r) B_{q,kx}(x) B_{r,ky}(y)
%
% to the original data
%
%    (X(i),Y(j)) |REPLACE_WITH_DASH_DASH>  f(X(i),Y(j)) = Z(i,j)
%
% obtained coincides with the earlier one (which generated COEFS rather than
% COEFSB).

%%
% 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.

tsp = spap2({knotsx,knotsy},[kx,ky],{x,y},z);
valuet = fnval(tsp,{xv,yv});
%%
% Here, for another check, is the relative difference between the values
% computed earlier and those computed now:
disp( max(max(abs(values-valuet)))/max(max(abs(values))) )


%% Interpolation
% 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.
%
% To recall, the data points were
% x = sort([(0:10)/10,.03 .07, .93 .97]);
% y = sort([(0:6)/6,.03 .07, .93 .97]);
%
% 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 cubics in  x , and use the not-a-knot condition, therefore use
% all but the second and the second-last data point as knots:

knotsx = augknt( x([1,3:(end-2),end]), kx );
spi2 = spapi(knotsx,x,coefsy.');

icoefs = fnbrk(spi2,'c').';

%% Evaluation
% We are ready to evaluate and plot the interpolant at a fine mesh:

ivalues = spcol(knotsx,kx,xv)*icoefs*spcol(knotsy,ky,yv).';

mesh(xv,yv,ivalues.'), view(150,50)
title('The spline interpolant')

%% Error
% Its error, as an approximation to the Franke function, is computed next:

fvalues = franke(repmat(xv.',1,length(yv)),repmat(yv,length(xv),1));
error =  fvalues - ivalues;
mesh(xv,yv,error.'), view(150,50)
title('Interpolation error')

%% Relative error
% The error is of  r e l a t i v e  size ...
disp( max(max(abs(error)))/max(max(abs(fvalues))) )

%%
% The above steps can be carried out by just two commands, one for the
% construction of the interpolant, the other for its evaluation
% at a rectangular mesh, as shown below.
% For a check, we also compute the relative difference between the values
% computed earlier and those computed now:

tsp = spapi({knotsx,knotsy},{x,y},z);
valuet = fnval(tsp,{xv,yv});
disp( max(max(abs(ivalues-valuet)))/max(max(abs(ivalues))) )


displayEndOfDemoMessage(mfilename)

##### SOURCE END #####
-->
   </body>
</html>