% HW 13 - Problem 5 - part (q2).

clear all

% Part 1:
% Investigate amplification factor for a fixed h.
h = 0.05;     
r2 = 1/(2*h);     
r4 = 1/(4*h);
theta1 = 1;      % backward Euler
theta2 = 1/2;    % Crank-Nicolson
betah = [0 : 0.01 : pi];
zz = 2*(cos(betah) - 1);
rho_theta1_r2 = (1 + r2*(1-theta1)*zz)./(1 - r2*theta1*zz);
rho_theta1_r4 = (1 + r4*(1-theta1)*zz)./(1 - r4*theta1*zz);
rho_theta2_r2 = (1 + r2*(1-theta2)*zz)./(1 - r2*theta2*zz);
rho_theta2_r4 = (1 + r4*(1-theta2)*zz)./(1 - r4*theta2*zz);
rho_true_r2 = exp( -r2*betah.^2);
rho_true_r4 = exp( -r4*betah.^2);

figure(130504);
plot(betah,rho_theta1_r2, ':', betah,rho_theta1_r4, '-.', ...
     betah,rho_theta2_r2, '--', betah,rho_theta2_r4, '-', 'linewidth',3)
hold on
plot(betah,rho_true_r2, 'k:', betah,rho_true_r4, 'k-.', 'linewidth',2)
hold off
xlabel('\beta h', 'fontsize',14);
ylabel('\rho', 'fontsize',14);
legend('\theta=1, r=1/2h', '\theta=1, r=1/4h', ...
      '\theta=0.5, r=1/2h', '\theta=0.5, r=1/4h', ...
      'exact, r=1/2h', 'exact, r=1/4h', 'fontsize',14)

  
% Part 2:
% Investigate amplitucation factor of the fastest harmonic
% as a function of  h. Do this for theta=1/2 (C-N) only.

hvec = [0.01 : 0.005 : 0.05];
kappa2 = hvec/2;
kappa4 = hvec/4;
rvec2 = kappa2./hvec.^2;
rvec4 = kappa4./hvec.^2;
tmax = 0.5;            
N2 = round(tmax./kappa2);
N4 = round(tmax./kappa4);

% Define the amplification factor |rho|:
log_totalrho2 = N2.*log10(abs( (1 - 2*rvec2)./(1 + 2*rvec2) ));
log_totalrho4 = N4.*log10(abs( (1 - 2*rvec4)./(1 + 2*rvec4) ));

%Plot |rho|^(t/k)  vs  r:
figure(130505);
plot(hvec, log_totalrho2, '--', hvec,  log_totalrho4, '-', 'linewidth',2)
xlabel('h', 'fontsize',14)
ylabel('log_{10} ( |\rho|^{tmax / \kappa} )', 'fontsize',13)
legend('r=1/2h', 'r=1/4h')

