% Matlab batch file to plot amplification factors for Fully Implicit
% (BTCS) finite-difference scheme for model parabolic problem
% 
%     u_t = u_xx 
% 
% vs true exponential-decay factors for continuous problem:
% 
%                           1
%     lambda(k) = ----------------------
%                 1 + 4*nu*sin(k*dx/2)^2
% vs
% 
%     lambda(k) = exp(-nu*(k*dx)^2)
% 
% Used in Math/CS 6/72262 "Numer PDE" course.

% File: ~/matlab/courses/lambda_k.m
% 
% C.G. (2-00)

n = 128 ;                             % # grid cells for plotting resolution

kdx = linspace( 0, pi, n ) ;          % uniform grid

clf ;  hold on ;                      % initialize figure (plotting window)

y = sin(kdx/2).^2 ;                   % vector of values used in several plots

for nu = [ 1/4 1/2 1 4 ]              % range of nu := dt / dx^2

  yh = 1 ./ ( 1 + 4 * nu * y ) ;      % discrete lambda(k) --- solid line

  ytrue = exp( - nu * kdx.^2 ) ;      % continuous/exact lambda(k) --- dashed

  plot( kdx, yh, 'b-', kdx, ytrue, 'r--' ) ;          % plot both on same plot

end

set( gca, ...
    'FontSize'  , 12                            , ...     % .. bigger fonts ..
    'XTick'     , [ 0 pi/4 pi/2 3*pi/4 pi ]     , ...     % 'x' axis ticks and
    'XTickLabel', '0|pi / 4|pi / 2|3 pi / 4| pi'  ...     % labels as in class
   ) ;

axis( [ 0 pi 0 1 ] ) ;

xlabel( '\fontsize{14}k\Delta{}x' ) ;                         % 'x' axis label

ylabel( '\fontsize{14}\lambda(k)' ) ;                         % 'y'  ""    "

title( '\fontsize{14}Amplification Factors:  \nu = 1/4, 1/2, 1, 4' ) ;

legend( 'discrete', 'continuous' ) ;  % legend in upper right corner (default)

box on ;