% lininva.m
%
% Solve M604D HW07 Question B
%
% Philip D. Loewen
% 2 May 2001

saveplots = 0;      % 0-1 variable: Make *.eps files (y/n)?
format compact; format long;

% Initialize constants specified on problem sheet
a = 1/4; b = 1/2;   % Interval endpoints
%c = 1;              % Exponential decay rate will be set later
mmax = 15;          % Maximum problem dimension

% PART B:  Multiplots show effects of decay rate and more data.

% Set up basis function values in a big array useful for plotting
numt = 201;                           % Number of t-values to plot with
tt = linspace(0,1,numt);              % Vector of t-values
ttmtx = ones(mmax,1)*tt;              % Big matrix of t-values, all rows equal
kkmtx = (1:mmax)'*ones(1,numt);       % Big matrix of k-values, all cols equal

for c = [0.25,1,4],
  disp(' ');
  disp('---------------------------------------------------------');
  disp([' Question B:  Case c = ',sprintf('%3.2f',c),'.']);
  disp('---------------------------------------------------------');
  disp(' ');
  
  % Copy PART A(b):  Generate Data Vector d of maximum size.

  k = (1:mmax)';
  d = ( gamma(k) ./ (c*ones(size(k))).^k ) ...
      .* (gammainc(c*b,k)-gammainc(c*a,k))
  disp(' ');
  
  % From PART A(c):  Build matrix Gamma = G * G' of maximum size.

  ss = 1:(2*mmax-1);
  ggtvals = gamma(ss) .* gammainc(2*c,ss) ./ (2*c*ones(size(ss))).^ss;
  J = ones(mmax,1)*(1:mmax) + (1:mmax)'*ones(1,mmax) - 1;
  GGT = ggtvals(J);
  GT = ttmtx.^(kkmtx-1) .* exp(-c*ttmtx); % Row j holds vals of gj(t)
  
  ts = ['Min-norm function (exact data fit), c=',...
        sprintf('%3.2f',c),'.'];
  figure('Name',ts); hold on;
  title(ts);
  xlabel('t');
  ylabel('x');
  crvhndls = plot([0,a,a,b,b,1],[0,0,1,1,0,0],':'); % exact
  mvals = [9,11,13];
  for m=mvals,
    lambda = GGT(1:m,1:m)\d(1:m);
    xx = lambda'*GT(1:m,:);
    ph = plot(tt,xx,'LineWidth',1);
    crvhndls = [crvhndls,ph];
    
    N  =  0.5*lambda'*GGT(1:m,1:m)*lambda;
    mf =  GGT(1:m,1:m)*lambda - d(1:m);
    E  =  0.5*mf'*mf;   % Recall S=I (exact data) assumed here
    disp([' m = ',sprintf('%2.0f',m),' gives a curve x with ',...
        ' E(x) = ',sprintf('%12.4e',E),', ',...
        ' N(x) = ',sprintf('%12.4e',N),'.']);
  end;
  set(crvhndls(2),'LineStyle','--','Color','r');
  set(crvhndls(3),'LineStyle','-.','Color','c');
  legend(crvhndls,'Exact',['m=',num2str(mvals(1))],...
                          ['m=',num2str(mvals(2))],...
                          ['m=',num2str(mvals(3))]);
  set(gca,'YLim',[-0.4,1.2]);                      
                        
  if saveplots,
    set(gcf,'PaperPosition',[1,1,4,3]);
    cstr = num2str(100*c);
    while length(cstr)<3, cstr=['0',cstr]; end;
    fname = ['m',num2str(m),'c',cstr,'.eps'];
    disp(' ');
    disp([' Figure number ',num2str(gcf),' saved as ',fname,'.']);
    print(gcf,'-deps',fname);
  end;

end;

disp(' ');
return;


