%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  an example Matlab script to setup, visualize and solve a large
%  tridiagonal system as described in class i.e. the "fish-farm"
%  problem
%  	$Id: Tridiagonal_LU.m,v 1.2 2002/09/17 03:24:56 mspieg Exp $	
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all; 
clear all;
%
%  set the size of the problem (i.e. number of fish-tanks)
%

N=11;

%
%  set up sparse tridiagonal matrix A (which is mostly 2's on the diagonal and -1 
%  on the first super and subdiagonal
%
A=spalloc(N,N,3*N);  % allocate room for a NxN sparse matrix with
                               % max 3N entries

% 
% set the matrix for the Boundary tanks 1 and N
%
A(1,1)=1;
A(N,N)=1;

%
% set the interior rows of the matrix (p.s. this isn't the fastest way
% to do this but it might be the clearest
%

for k=2:N-1
  A(k,k-1:k+1)=[ -1 2 -1 ];  % just sets the 3 non-zero elements of row k
end

%
% Let's see what that looks like
%

figure
spy(A);
title('Tridiagonal matrix A');
disp(sprintf('\n Here are the first 5 rows of A\n'));
disp(full(A(1:5,1:5)))
disp('hit return to continue ');
pause;

%
%  Now let's do an LU decomposition on A  
%

[L,U]=lu(A);  % doesn't get much easier than that

%
%  and take a look
%

figure
subplot(1,2,1); spy(L); title('L');
subplot(1,2,2); spy(U); title('U');
disp(sprintf('\n Here are the first 5 rows of L and U\n'));
disp(full(L(1:5,1:5)))
disp(full(U(1:5,1:5)))
disp('hit return to continue ');
pause;

%
%  Now let's just solve the problem for a single input in tank i
%

i=4;
r=zeros(N,1); % start by make r a zero vector of length N
r(1)=1; % set the height of tank 1 to 1
r(N)=1; % set the height of tank N to 1
r(i)=1;       % and put a flux of 1 in component i

%
%  now solve by forward substitution for c
%

c=L\r;

%
%  and back substitution for the vector of heights h
%

h=U\c;


%
% now plot out the right hand side and the solution x
%
figure
bar(h,'b'); % plot water height as a bar plot
hold on;    
bar([0 ; r(2:N-1) ; 0],'r'); % plot the forcing function RHS (without the end points)
xlabel('Tank number');
ylabel('height or flux');
legend('water height','input flux');
disp('hit return to continue ');
pause;


%
%  lets do the same thing for a random RHS
%

r=zeros(N,1);
r([1 N])=1; % set the boundary tank heights again
r(2:N-1)=rand(N-2,1);  % generate random forcings between 0 and 1

%
%  and solve again, this time in one shot
%

h=U\(L\r);


figure
bar(h,'b'); % plot water height as a bar plot
hold on;    
bar([0 ; r(2:N-1) ; 0],'r'); % plot the forcing function RHS (without the end points)
xlabel('Tank number');
ylabel('height or flux');
legend('water height','input flux');
disp('hit return to continue ');
title('Solution of L*U*h=r');
pause;

%
%  Finally, let's do this the wrong way by using the matrix inverse
%

Ainverse=inv(A);  % calculate the inverse...also very easy

%
%  look at it...it's rather un-sparse
%

figure
subplot(1,2,1); spy(A); title('A');
subplot(1,2,2); spy(Ainverse); title('A^{-1}')
disp(sprintf('\n Here are the first 5 rows of inv(A)\n'));
disp(full(Ainverse(1:5,1:5)))
disp('hit return to continue ');
pause;

%
%  but it still should give a reasonable answer
%

h=inv(A)*r;

figure
bar(h,'b'); % plot water height as a bar plot
hold on;    
bar([0 ; r(2:N-1) ; 0],'r'); % plot the forcing function RHS (without the end points)
xlabel('Tank number');
ylabel('height or flux');
legend('water height','input flux');
title('Solution of h=inv(A)*r');