clear

t_0 = 0;
t_f = .1;
M = 16;
N = 8;

tic
[u1,err1,x1,t1] = heat_eul_dir(t_0,t_f,M,N);
rtime(1)=toc;

M = 4*M;
N = 2*N;
tic
[u2,err2,x2,t2] = heat_eul_dir(t_0,t_f,M,N);
rtime(2) = toc;

M = 4*M;
N = 2*N;
tic
[u3,err3,x3,t3] = heat_eul_dir(t_0,t_f,M,N);
rtime(3)=toc;


M = 4*M;
N = 2*N;
tic
[u4,err4,x4,t4] = heat_eul_dir(t_0,t_f,M,N);
rtime(4)=toc;

M = 4*M;
N = 2*N;
tic
[u5,err5,x5,t5] = heat_eul_dir(t_0,t_f,M,N);
rtime(5)=toc;


M = 4*M;
N = 2*N;
tic
[u6,err6,x6,t6] = heat_eul_dir(t_0,t_f,M,N);
rtime(6)=toc;

M = 4*M;
N = 2*N;
tic
[u7,err7,x7,t7] = heat_eul_dir(t_0,t_f,M,N);
rtime(7)=toc;

M = 4*M;
N = 2*N;
tic
[u8,err8,x8,t8] = heat_eul_dir(t_0,t_f,M,N);
rtime(8)=toc;


% at each moment in time, compute the L^infty norm of the error
E1 = max(abs(err1));
E2 = max(abs(err2));
E3 = max(abs(err3));
E4 = max(abs(err4));
E5 = max(abs(err5));
E6 = max(abs(err6));
E7 = max(abs(err7));
E8 = max(abs(err8));
% plot these errors versus time for the 5 computations
figure(1)
clf
plot(t1,E1,':','LineWidth',2);
hold on
plot(t2,E2,'--','LineWidth',2);
plot(t3,E3,'-.','LineWidth',2);
plot(t4,E4,'-','LineWidth',2);
plot(t5,E5,'r:','LineWidth',2);
plot(t6,E6,'r--','LineWidth',2);
plot(t7,E7,'r-.','LineWidth',2);
plot(t8,E8,'r-','LineWidth',2);
title('errors versus time')

% for a more informative plot, plot the log10 of the errors
figure(2)
clf
plot(t1,log10(E1),':','LineWidth',2);
hold on
plot(t2,log10(E2),'--','LineWidth',2);
plot(t3,log10(E3),'-.','LineWidth',2);
plot(t4,log10(E4),'-','LineWidth',2);
plot(t5,log10(E5),'r:','LineWidth',2);
plot(t6,log10(E6),'r--','LineWidth',2);
plot(t7,log10(E7),'r-.','LineWidth',2);
plot(t8,log10(E8),'r-','LineWidth',2);
title('log10(error) versus time')

% It appears that the errors are decreasing --- how to find the rate?  I
% have to be careful because twice as many time-intervals were used in
% computing E2 as were used in computing E1.  So If I'm going to do ratios
% of E1/E2, I need to sample E2 at the right times.  I do thise by defining
% the index I and then computing E2(I).  To check that this works, ask
%
% matlab for t1(1:5) and t2(I(1:5))
I = 1:4:length(t2);
rat1 = E1./E2(I);

I = 1:4:length(t3);
rat2 = E2./E3(I);

I = 1:4:length(t4);
rat3 = E3./E4(I);

I = 1:4:length(t5);
rat4 = E4./E5(I);

I = 1:4:length(t6);
rat5 = E5./E6(I);

I = 1:4:length(t7);
rat6 = E6./E7(I);

I = 1:4:length(t8);
rat7 = E7./E8(I);

% now plot the ratios as a function of time.  They are going to 4, as
% expected seeing as forward Euler has error O(k) if I'm on a fixed
% refinement path and I'm doing k -> k/4 at each refinement...
figure(3)
clf
plot(t1,rat1,':','LineWidth',2)
hold on
plot(t2,rat2,'--','LineWidth',2)
plot(t3,rat3,'-.','LineWidth',2);
plot(t4,rat4,'-','LineWidth',2)
plot(t5,rat5,'r:','LineWidth',2)
plot(t6,rat6,'r--','LineWidth',2)
plot(t7,rat7,'r-.','LineWidth',2)
title('error ratios versus time')

% to get slope more quantitatively, I can take the norm of the error at the
% final time and see how that depends on h...
dx(1) = x1(2)-x1(1);
dx(2) = x2(2)-x2(1);
dx(3)=dx(2)/2;
dx(4)=dx(3)/2;
dx(5)=dx(4)/2;
dx(6) = dx(5)/2;
dx(7) = dx(6)/2;
dx(8) = dx(7)/2;
% take the errors at the final times
E(1) = E1(end);
E(2)=E2(end);
E(3)=E3(end);
E(4)=E4(end);
E(5)=E5(end);
E(6) = E6(end);
E(7) = E7(end);
E(8) = E8(end);
figure(4)
clf
plot(log10(dx),log10(E),'o-')
xlabel('log10(h)')
ylabel('log10(||error(final time)||)')
hold on
% compute the least squares fit
X = log10(dx(1:8));
Y = log10(E(1:8));
clear A
A(:,1)=X;
A(:,2)= ones(size(X));
B = A'*Y';
A = A'*A;
mm = inv(A)*B;
hold on
% plot the least squares fit
plot(X,mm(1)*X+mm(2),'r')
title('log10(error) versus log10(h)')
% if mm(1) is close to 2, we're happy:
mm(1)


figure(5)
clf
plot(log10(dx),log10(rtime),'o-')  
X = log10(dx(4:8));
Y = log10(rtime(4:8));
clear A
A(:,1)=X;
A(:,2)= ones(size(X));
B = A'*Y';
A = A'*A;
mm = A\B;
hold on
plot(X,mm(1)*X+mm(2),'r')
xlabel('log10(dx)')
ylabel('log10(runtime)')
title('log10(runtime) versus log10(h) along with least-squares fit')




% Now I can test what would happen if I didn't have the errors but I just
% had u_coarse, u_medium, and u_fine.  Now I have to 
I1 = 1:2:length(x2);
J1 = 1:4:length(t2);
I2 = 1:4:length(x3);
J2 = 1:16:length(t3);
diff1 = u1-u2(I1,J1);
diff2 = u2(I1,J1)-u3(I2,J2);
Rat1 = diff1./diff2;
r1 = max(abs(diff1))./max(abs(diff2));

I1 = 1:2:length(x3);
J1 = 1:4:length(t3);
I2 = 1:4:length(x4);
J2 = 1:16:length(t4);
diff1 = u2-u3(I1,J1);
diff2 = u3(I1,J1)-u4(I2,J2);
Rat2 = diff1./diff2;
r2 = max(abs(diff1))./max(abs(diff2));

I1 = 1:2:length(x4);
J1 = 1:4:length(t4);
I2 = 1:4:length(x5);
J2 = 1:16:length(t5);
diff1 = u3-u4(I1,J1);
diff2 = u4(I1,J1)-u5(I2,J2);
Rat3 = diff1./diff2;
r3 = max(abs(diff1))./max(abs(diff2));

I1 = 1:2:length(x5);
J1 = 1:4:length(t5);
I2 = 1:4:length(x6);
J2 = 1:16:length(t6);
diff1 = u4-u5(I1,J1);
diff2 = u5(I1,J1)-u6(I2,J2);
Rat4 = diff1./diff2;
r4 = max(abs(diff1))./max(abs(diff2));

I1 = 1:2:length(x6);
J1 = 1:4:length(t6);
I2 = 1:4:length(x7);
J2 = 1:16:length(t7);
diff1 = u5-u6(I1,J1);
diff2 = u6(I1,J1)-u7(I2,J2);
Rat5 = diff1./diff2;
r5 = max(abs(diff1))./max(abs(diff2));

I1 = 1:2:length(x7);
J1 = 1:4:length(t7);
I2 = 1:4:length(x8);
J2 = 1:16:length(t8);
diff1 = u6-u7(I1,J1);
diff2 = u7(I1,J1)-u8(I2,J2);
Rat6 = diff1./diff2;
r6 = max(abs(diff1))./max(abs(diff2));

figure(6)
clf
plot(t1,r1,':')
hold on
plot(t2,r2,'--')
plot(t3,r3,'-.')
plot(t4,r4,'-')
plot(t5,r5,'r:')
plot(t6,r6,'r--')
xlabel('time')
title('plot the ratios of the norms of the solution differences versus time')

figure(7)
clf
plot(x1,Rat1(:,end),':')
hold on
plot(x2,Rat2(:,end),'--')
plot(x3,Rat3(:,end),'-.')
plot(x4,Rat4(:,end),'-')
plot(x5,Rat5(:,end),'r:')
plot(x6,Rat6(:,end),'r--')
axis([0,2*pi,-2,6])
xlabel('space')
title('plot the Richardson ratios (at the final time) versus space')




