% HW 15, # 2  (version 2010 of the HW;  previously this was Problem 3)
% Simple Explicit Method for the Heat equation in 2 dimensions with non-zero Boundary Conditions

% HERE WE SOLVE FOR U(x,t) BY "FILLING THE BOX"-METHOD.

format short e

h=0.1;
k=0.002;
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);         % M+1  in the notes 
Lp1=length(y);         % L+1  in the notes

tmin=0;
tmax=0.5;
t=tmin:k:tmax;
N=length(t);

% Open the movie file:
nFrames = N;         % upper bound for the # of frames in the movie 
                     % (this is more than needed since not every frame is
                     %  recorded)
mymovie = VideoWriter('hw15_p2.avi'); 
mymovie.Quality = 70;  % default is 75
mymovie.FrameRate = 2; % frames per second when movie is played
open(mymovie);


% The initial condition (the "bottom of the box"):
A1=10;
U0 = A1*sin(pi*X/xmax).*sin(pi*Y/ymax) + sin(2*pi*Y/ymax).*(1-X) + cos(2*pi*Y/ymax).*X; 
                   % Note that at the "bottom's edges", this agrees with the boundary conditions.

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


% Plot the initial condition and prepare the axes, labels, etc:
figure(15021);
subplot(2,1,1);
mesh(x,y,U);
axis([0 xmax 0 ymax min(min(U))*1.1 max(max(U))*1.1]);
title('Computed solution','fontsize',10)
subplot(2,1,2);
mesh( x,y,zeros(size(X)) );              % plot of the initial error
axis([0 xmax 0 ymax 0 10^(-3) ]);  
title('Numerical error','fontsize',10)

% Initialize the movie's frames:
writeVideo(mymovie, getframe(gcf));



% Begin the evolution in time:
for n=1:N-1;
    
    % Since the boundary conditions (the "sides of the box") are time-dependent, I create them inside the time loop:
	Unew(:,1) = sin(2*pi*y'/ymax) * exp( -(2*pi/ymax)^2*t(n+1) );     % along  x=0
	Unew(:,Mp1) = cos(2*pi*y'/ymax) * exp( -(2*pi/ymax)^2*t(n+1) );     % along  x=xmax
	Unew(1,:) = x' * exp( -(2*pi/ymax)^2*t(n+1) );     % along  y=0
	Unew(Lp1,:) = x' * exp( -(2*pi/ymax)^2*t(n+1) );     % along  y=ymax
    
	for m=2:Mp1-1
        for ell=2:Lp1-1
            Unew(ell,m) = U(ell,m)*(1-4*r) + r*( U(ell,m-1) + U(ell,m+1) + U(ell-1,m) + U(ell+1,m) );
        end
	end
    
    U = Unew;
    	
    rec_frame=4;      % 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 = A1*sin(pi*X/xmax).*sin(pi*Y/ymax)*exp(-(1/xmax^2+1/ymax^2)*pi^2*t(n+1)) + ...
                 ( sin(2*pi*Y/ymax).*(1-X) + cos(2*pi*Y/ymax).*X ) * exp( -(2*pi/ymax)^2*t(n+1) ) ;
        
		% Plot the computed solution and its error:
		subplot(2,1,1);
        mesh(x,y,U);
        axis([0 xmax 0 ymax  min(min(U))*1.1  max(max(U))*1.1 ]);
		title('Computed solution','fontsize',10)
		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',10)
        pause(0.3);

        writeVideo(mymovie, getframe(gcf));
        
    end 
end

pause(0.2)

Uexact = A1*sin(pi*X/xmax).*sin(pi*Y/ymax)*exp(-(1/xmax^2+1/ymax^2)*pi^2*t(n+1)) + ...
        ( sin(2*pi*Y/ymax).*(1-X) + cos(2*pi*Y/ymax).*X ) * exp( -(2*pi/ymax)^2*t(n+1) ) ;
% Plot the computed solution and its error:
subplot(2,1,1);
mesh(x,y,U);
axis([0 xmax 0 ymax  min(min(U))*1.1  max(max(U))*1.1 ]);
title(['Computed solution;  Umax = ' num2str(max(max(U))) ],   'fontsize',10)
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',10)
pause(0.5);

writeVideo(mymovie, getframe(gcf));
        
close(mymovie);       

% figure(15022);
%         mesh(x,y,Uexact - U);
%         vert_extent=max(max(abs(U-Uexact)))*1.1 + 10^(-8);   % see below
% 		axis([0 xmax 0 ymax -0.1*vert_extent  vert_extent ]);  
% 		title(['Numerical error of magnitude   ' num2str(max(max(abs(U-Uexact))))],'fontsize',14)
%         xlabel('X', 'fontsize',12)
%         ylabel('Y', 'fontsize',12)
%         zlabel('U_{exact} - U', 'fontsize',12)

        