#!/usr/bin/octave -q

%-----------------------------------------------------------------------------
% Jonathan R. Senning <senning@gordon.edu>
% Gordon College
% April 22, 1999
%
% $Id: crank_nicolson,v 1.7 2015/04/24 14:01:03 senning Exp senning $
%
% Use Crank-Nicolson scheme to solve the heat equation in a thin rod.
% The equation solved is
%
%	du     d  du
%	-- = K -- --
%	dt     dx dx
%
% along with boundary conditions
%
%	u(xmin,t) = a(t)
%	u(xmax,t) = b(t)
%
% and initial conditions
%
%	u(x,tmin) = f(x)
%
%-----------------------------------------------------------------------------

% Determine the desired plot types.  A value of 1 means show the plot and
% a value of 0 means don't show the plot.

individual_curves = 1;		% Show curves during computation
color_image = 1;		% Show pseudocolor plot of solution
surface_mesh = 1;		% Show surface plot of solution
contour_plot = 1;		% Show contour plot of solution

% Set the value of the thermal diffusivity.

K = 0.25;

% Set size of image.  Also, these values are combined with interval sizes to
% compute h (spatial stepsize) and k (temporal stepsize).  Define array
% to hold image.

n = 64;			% Number of spatial intervals
m = 64;			% Number of temporal intervals
u = zeros( n+1, m+1 );

% X position of left and right endpoints

xmin = 0;
xmax = 1;

% Interval of time: tmin should probably be left at zero

tmin = 0;
tmax = 2;

% Generate x and t values.  These aren't really needed to solve the PDE,
% but they are useful for computing boundary/initial conditions and graphing.

x = linspace( xmin, xmax, n+1 );
t = linspace( tmin, tmax, m+1 );

% Initial condition

u(:,1) = 100 * sin( pi * x )';

% Boundary conditions

u(1,:) = zeros( 1, m+1 );				% Left
u(n+1,:) = 60 * ( ( 1 .- cos( pi * t ) ) / 2.0 );	% Right

% We are using a Crank-Nicolson scheme, and can vary the weighting of the
% u(x,t+k) values with the u(x,t) values using the parameter c.
%
% If c == 0 then the iteration reduces to the fully explicit marching scheme
% which is not stable unless k <= (h^2)/2.
%
% If c == 1 then the iteration is fully implicit and is stable for any choice
% of k.
%
% The method should be unconditionally stable if c >= 0.5.  However, when
% c == 0.5 the solution can have undesirable oscillation.

c = 0.75;

%-----------------------------------------------------------------------------
%	Should not need to make changes below this point :)
%-----------------------------------------------------------------------------

% Compute step sizes.  Since we also need the square of the spatial step
% size, we go ahead and compute it now as well.

h = ( xmax - xmin ) / n;
k = ( tmax - tmin ) / m;
h2 = h * h;

% Find likely extremes for u

umin = min( min( u ) );
umax = max( max( u ) );

% Now we are basically done with the setup.  The discritation used leads to
% a tridiagonal linear system.  Here we construct the diagonals of tridiagonal
% matrix.

D = ( 2 * k * K * c + h2 ) * ones( n-1, 1 );
A = -k * K * c * ones( n-2, 1 );
C = -k * K * c * ones( n-2, 1 );

% Plot initial condition curve.  The "pause()" is used to allow time for the
% plot to appear on the screen before actually starting to solve the problem
% for t > 0.

if ( individual_curves != 0 )
    fmt = sprintf( 'o-;t = %5.2f;', 0 );
    plot( x, u(:,1), fmt );
    axis( [xmin, xmax, umin, umax] );
%    xlabel( 'x' );
%    ylabel( 'Temperature' );
    drawnow();
    pause( 2 );
end

% Main loop.  This consists of computing the appropriate right-hand-side
% vector and then solving the linear system.  We also plot the solution at
% each time step.  Note that this is not very efficient -- some of the work
% done as part of the tridiagonal solve processes need only be done once and
% could be done outside of this loop.

for j = 1 : m

    B = k * K * ( 1 - c ) * ( u(1:n-1,j) + u(3:n+1,j) ) ...
		- ( 2 * k * K * ( 1 - c ) - h2 ) * u(2:n,j);
    B(1)   = B(1)   + k * K * c * u(1,j+1);
    B(n-1) = B(n-1) + k * K * c * u(n+1,j+1);

    u(2:n,j+1) = tridiagonalSolver( A, D, C, B )';

    if ( individual_curves != 0 )
        fmt = sprintf( 'o-;t = %5.2f;', j * k );
        plot( x, u(:,j+1), fmt );
        axis( [xmin, xmax, umin, umax] );
%        xlabel( 'x' );
%        ylabel( 'Temperature' );
        drawnow();
    end

end

if ( individual_curves != 0 )
    input( 'Press Enter to continue... ' );
end
   
% All done computing solution, now show the desired plots...
%
% Alter the colormap to get a nice one and then display the image.  In the
% image, x is the horizontal axis and t is the vertical axis - however, t
% increases in the downward direction.  Thus, the first line of the image
% corresponds to the initial condition.

if ( color_image != 0 || surface_mesh != 0 || contour_plot != 0 )
    cm = colormap( 'default' );
    cm(:,3)=(1.0.-[0:63]/63)';
    cm(:,2)=4*cm(:,1).*(1.-cm(:,1));
    colormap( cm );
end

if ( color_image != 0 )
    image( x, t, u' );
    xlabel( 'x' );
    ylabel( 't' );
    drawnow();
    input( 'Press Enter to continue... ' );
end

% Display a mesh (surface) plot showing the solution.

if ( surface_mesh != 0 )
    meshc( x, t, u' );
    view( 37.5, 45 );
    xlabel( 'x' );
    ylabel( 't' );
    drawnow();
    input( 'Press Enter to continue... ' );
end

% We could also create a contour graph with the following command.  Note that
% correct contour levels depend on the solution itself.

if ( contour_plot != 0 )
    umin = min( min( u ) );
    umax = max( max( u ) );
    levels = linspace( umin, umax, 21 );
    contour( x, t, u', levels );
    xlabel( 'x' );
    ylabel( 't' );
    drawnow();
    input( 'Press Enter to continue... ' );
end

% End of file