% HW 15, # 9
% Peaceman-Rachford method for the Heat equation in 2 dimensions 
% with non-zero, time-dependent boundary conditions.
%
% Here I did NOT include movie-making commands in order to focus one's
%  attention on the details of the P-R method.  

clear all

format short e

h = 0.1;
k = h/4; 
r = k/h^2;
xmax = 1;
x = 0:h:xmax;
ymax = 2.5;
y = 0:h:ymax;
[X,Y] = meshgrid(x,y);   
Mp1=length(x);   % Mp1 = M+1 (in accordance with the notations in Lecture 15)
Lp1=length(y);   % Lp1 = L+1
M = Mp1-1;
L = Lp1-1;

tmin = 0;
tmax = 0.5;
t = tmin+k :k: tmax;
N = length(t);  % N (in Lecture 15's notations)


scale=1;         % This parameter allows one to scale the nonzero part 
                 %  of the boundary conditions
                 %  and, correspondingly, a piece of the initial condition.
                 %  scale =1 correspond to the original Problem 2 (or 12), while
                 %  scale=0 corresponds to Problem 1 (or 11) 
				 %  where the solution is 10 times bigger than there.
                 %
                 % The reason I introduced this factor is to facilitate debugging the problem.
                 %  Namely, if I set scale=0 and see a problem, then it does NOT come from 
                 %  the boundary conditions. On the other hand, if the code works for scale=0
                 %  but does not for scale=1, then the problem IS with the boundary conditions.
                 
% The initial condition:
Uinit = 10*sin(pi*X/xmax).*sin(pi*Y/ymax) + ...
        scale*( (1-X/xmax).*sin(2*pi*Y/ymax) +(X/xmax).*cos(2*pi*Y/ymax) );
                      % Note that at the "bottom's edges", this agrees with the boundary conditions.

U=Uinit;              % U at time level n=1


% Prepare vectors that will be used in the Thomas algorithm:
ax = -r/2*ones(1,M-2);    % vector  a  in  delta_x^2
bx=(1+r)*ones(1,M-1);   % vector  b  in  delta_x^2
cx=-r/2*ones(1,M-2);    % vector  c  in  delta_x^2
ay=-r/2*ones(1,L-2);    % vector  a  in  delta_y^2
by=(1+r)*ones(1,L-1);   % vector  b  in  delta_y^2
cy=-r/2*ones(1,L-2);    % vector  c  in  delta_y^2


% Plot the initial condition etc:
figure(15121);
subplot(2,1,1);
mesh(x,y,U);
% axis([0 xmax 0 ymax 0 max(max(U))*1.1]);
title('Numerical solution','fontsize',12)
subplot(2,1,2);
mesh(x,y,zeros(Lp1,Mp1));              % plot of the initial error
% axis([0 xmax 0 ymax 0 10^(-3) ]);  
title('Numerical error','fontsize',12)
% pause   % for monitoring/debugging purposes only

% Begin the evolution in time:
for n=1:N;

    % Step 1:
    % Set the boundary conditions at levels  n+1 and "star".
    %  The b.c. for the level  n  are already set since we know the solution (denoted as U) at level  n.
    time_factor = exp(-(2*pi/ymax)^2*t(n));
	Unew(:,1) = scale* sin(2*pi*y'/ymax)*time_factor;     % along  x=0
	Unew(:,Mp1) = scale* cos(2*pi*y'/ymax)*time_factor;   % along  x=xmax
	Unew(1,:) = scale* x*time_factor;                     % along  y=0
	Unew(Lp1,:) = scale* x*time_factor;                   % along  y=ymax
    %
    Ustar(2:L,1) = 1/2*( U(2:L,1) + Unew(2:L,1) ) + ...   % at x=0 for y=h...ymax-h
                   r/4*( (U(3:Lp1,1) - 2*U(2:L,1) + U(1:L-1,1)) - ...
                         (Unew(3:Lp1,1) - 2*Unew(2:L,1) + Unew(1:L-1,1)) );
    Ustar(2:L,Mp1) = 1/2*( U(2:L,Mp1) + Unew(2:L,Mp1) ) + ...  % at x=xmax for y=h...ymax-h
                     r/4*( (U(3:Lp1,Mp1) - 2*U(2:L,Mp1) + U(1:L-1,Mp1)) - ...
                           (Unew(3:Lp1,Mp1) - 2*Unew(2:L,Mp1) + Unew(1:L-1,Mp1)) );
      
    % Step 2:
    % Solve tridiagonal systems along the x-direction.
    for ell=2:L
        rhsx = (1-r)*U(ell,2:M) + r/2*( U(ell+1,2:M) + U(ell-1,2:M) ) + ...
               r/2* [Ustar(ell,1) zeros(1,M-3) Ustar(ell,Mp1)];
        Ustar(ell,2:M) = thomas(ax,bx,cx,rhsx);
    end
    
    % Step 3:
    % Solve tridiagonal systems along the y-direction.
    for m=2:M
        rhsy = ( (1-r)*Ustar(2:L,m) + r/2*( Ustar(2:L,m+1) + Ustar(2:L,m-1) ) )' + ...
               r/2* [Unew(1,m) zeros(1,L-3) Unew(Lp1,m)];
        Unew(2:L,m) = ( thomas(ay,by,cy,rhsy) )';
    end
    
    U = Unew;
    Uexact = 10*sin(pi*X/xmax).*sin(pi*Y/ymax) * exp(-(1/xmax^2+1/ymax^2)*pi^2*t(n)) + ...
             scale*( (1-X/xmax).*sin(2*pi*Y/ymax) +(X/xmax).*cos(2*pi*Y/ymax) ) *exp(-(2*pi/ymax)^2*t(n));
    
   
    % Plot the solution and record frames for the movie:    
    rec_frame=4;      % plot and/or record only every "rec_frame"-th frame (to confine the size of the .avi file)
    if rem(n,rec_frame) == 0
    
        % The exact solution:
%         Uexact = scale*(1-Y/ymax) + sin(pi*X/xmax).*sin(pi*Y/ymax)*exp(-(1/xmax^2+1/ymax^2)*pi^2*t(n+1));
        
		% Plot the computed solution and its error:
		subplot(2,1,1);
        mesh(x,y,U);
        vert_extent=max(max(abs(U)))*1.1;
        axis([0 xmax 0 ymax -vert_extent vert_extent]);
		title('Numerical solution','fontsize',12)
		subplot(2,1,2);
        mesh(x,y,U-Uexact);
		vert_extent=max(max(abs(U-Uexact)))*1.1 + 10^(-8);   % see below
		axis([0 xmax 0 ymax -vert_extent  vert_extent ]);   
		title(['Numerical error of magnitude   ' num2str(max(max(abs(U-Uexact))))],'fontsize',12)
        % pause       % for monitoring/debugging purposes only

    end 
end

    
Uexact = 10*sin(pi*X/xmax).*sin(pi*Y/ymax) * exp(-(1/xmax^2+1/ymax^2)*pi^2*t(n)) + ...
         scale*( (1-X/xmax).*sin(2*pi*Y/ymax) +(X/xmax).*cos(2*pi*Y/ymax) ) *exp(-(2*pi/ymax)^2*t(n));
subplot(2,1,1);
mesh(x,y,U);
vert_extent=max(max(abs(U)))*1.1;
axis([0 xmax 0 ymax -vert_extent vert_extent]);
title(['Numerical solution;  Umax =   ' num2str(max(max(U)))],'fontsize',12)
subplot(2,1,2);
mesh(x,y,U-Uexact);
vert_extent=max(max(abs(U-Uexact)))*1.1 + 10^(-8);   % see below
axis([0 xmax 0 ymax -vert_extent  vert_extent ]);   
title(['Numerical error of magnitude   ' num2str(max(max(abs(U-Uexact))))],'fontsize',12)
