%
%   *******************   Collocation Project   ********************
%
% colloc.m
%
order = input('Maximum polynomial degree? ');
order = order + 1;                                  % Order = max deg + 1.
n = input('How many basis functions? ' )
knots = augknt( linspace(0,6,n-2), order )
m = input('How many sampling (design) points m? ' )
t = linspace(2,4,m)';                               % Column vector.
%
y = Ax(t);                                          % Get system output data.
G = Ggen('ksamp',t,knots,order);                    % Generate the G matrix.
c = G\y;                                            % Solve Gc=y in the
%                                                       least squares sense.
sp = spmak(knots,c');                               % Generate spline data
%                                                       stucture.  c'=row vec.
npts = 100;                                         % Number of points to plot.
s = linspace(0,6,npts);
v = zeros(size(s));
v = fnval(sp,s);
plot(s,v)


Ax.m
----
function y=Ax(t)
%Ax  Formula for ECE 735 integral operator output.
y = ( sin(6*(t-1))./(t-1) + sin(6*(t+1))./(t+1) )/2;


ksamp.m
-------
function y = ksamp(s)
global TOUT
y = cos(TOUT .* s);


Ggen.m
------
function y=Ggen(krnl,t,knots,order)
%Ggen  Generate the matrix
%
%             G(i,j) = int_knots(1)^knots(n) krnl(t(i),s) B_j(s) ds
%
%      for i = 1, ... , length(t), j = 1, ... , n,
%      where n = length(knots) - order
%      krnl is a string variable containing the name of the kernel
%           function.  Usage of krnl is  krnl(s),  where the first
%           argument is the global variable TOUT
%
global GKERNEL TOUT SPDATA
GKERNEL = krnl;
lt = length(t);
n = length(knots) - order;
y = zeros(lt,n);
coef = 1;
for j = 1:n
   fprintf( 'Working on column (basis function) %g\n', j )
   fprintf( '  Design point' )
   SPDATA = spmak( knots(j:j+order), coef );
   a = knots(j);
   b = knots(j+order);
   for i = 1:lt
      fprintf( ' %g', i )
      TOUT = t(i);
      y(i,j) = quad8('Gint',a,b);
   end
   fprintf( '\n' )
end


Gint.m
------
function y = Gint(s)
%Gint  called by Ggen
global GKERNEL SPDATA
y = feval(GKERNEL,s) .* fnval(SPDATA,s);


doreg.m
-------
%
% Regularization script to be run after colloc script if necessary.
%
lambda = input('Regularization parameter? ')
H=G'*G;                                             % G is put into the work-
                                                    % space by colloc
alpha = lambda*ones(1,size(H,2))
A = diag(alpha)+H;
c = A\(G'*y)
pause
sp = spmak(knots,c');                               % Generate new spline data
%                                                       stucture.  c'=row vec.
npts = 100;                                         % Number of points to plot.
s = linspace(0,6,npts);
v = zeros(size(s));
v = fnval(sp,s);
plot(s,v)
