

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% explicit timestepping

clear
t_0 = 0; t_f = 1; M = 200; N = 64;
L=2*pi;
h = L/N;
x = 0:h:L-h;

eps = 1;
u0 = 1./(eps^2 + sin(x/2).^2);

% phases = rand(1,10)*2*pi;
% amps = rand(1,10);
% nn = length(amps);
% u0 = zeros(size(x(1:N)));
% for j=1:nn
%   u0 = u0 + amps(j)*cos((j-1)*x(1:N)+phases(j));
% end

[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);

mm = max(max(abs(u)));

for j=0:15
plot(x,u(:,1+2*j))
axis([0,2*pi,0,2])
pause(1)
figure(1)
end
 
for j=0:15
plot(kk,amp(:,1+2*j),'o-')
axis([0,N/2,-17,2])
figure(1)
pause(1)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=2*M

[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
for j=0:15
plot(x,u(:,1+4*j))
axis([0,2*pi,0,2])
pause(1)
figure(1)
end
 
for j=0:15
plot(kk,amp(:,1+4*j),'o-')
axis([0,N/2,-17,2])
figure(1)
pause(1)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=2*M

[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
for j=0:15
plot(x,u(:,1+8*j))
axis([0,2*pi,0,2])
pause(1)
figure(1)

for j=0:15
plot(kk,amp(:,1+8*j),'o-')
axis([0,N/2,-17,2])
figure(1)
pause(1)
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% convergence study

% go back to unstable computation:
M = M/4;
[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
Err(1) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
Err(2) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
Err(3) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
Err(4) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
Err(5) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
Err(6) = max(abs(err(:,M+1)));
% M = 2*M;
% [u,err,x,t,kk,amp] = heat4(t_0,t_f,M,N,u0);
% Err(7) = max(abs(err(:,M+1)));

Err
Err(1:5)./Err(2:6)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% implicit time stepping
% 
% clear
t_0 = 0; t_f = 1; M = 100; N = 64;
% h = 2*pi/N;
% x = 0:h:2*pi;
% 
% phases = rand(1,10)*2*pi;
% amps = rand(1,10);
% nn = length(amps);
% u0 = zeros(size(x(1:N)));
% for j=1:nn
%   u0 = u0 + amps(j)*cos((j-1)*x(1:N)+phases(j));
% end

[u,err,x,t,kk,amp] = heat5(t_0,t_f,M,N,u0);

mm = max(max(abs(u)));

for j=0:15
plot(x,u(:,1+2*j))
axis([0,2*pi,0,1.5])
pause(1)
figure(1)
end
 
for j=0:15
plot(kk,amp(:,1+2*j),'o-')
pause(1)
figure(1)
axis([0,N/2,-17,2])
end

Err(1) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat5(t_0,t_f,M,N,u0);
Err(2) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat5(t_0,t_f,M,N,u0);
Err(3) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat5(t_0,t_f,M,N,u0);
Err(4) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat5(t_0,t_f,M,N,u0);
Err(5) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat5(t_0,t_f,M,N,u0);
Err(6) = max(abs(err(:,M+1)));
% M = 2*M;
% [u,err,x,t,kk,amp] = heat5(t_0,t_f,M,N,u0);
% Err(7) = max(abs(err(:,M+1)))

Err
Err(1:5)./Err(2:6)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Crank-Nicolson time stepping
% 
% clear
t_0 = 0; t_f = 1; M = 100; N = 64;
% h = 2*pi/N;
% x = 0:h:2*pi;
% 
% phases = rand(1,10)*2*pi;
% amps = rand(1,10);
% nn = length(amps);
% u0 = zeros(size(x(1:N)));
% for j=1:nn
%   u0 = u0 + amps(j)*cos((j-1)*x(1:N)+phases(j));
% end

[u,err,x,t,kk,amp] = heat6(t_0,t_f,M,N,u0);

mm = max(max(abs(u)));

for j=0:15
plot(x,u(:,1+2*j))
axis([0,2*pi,0,1.5])
pause(1)
figure(1)
end
 
for j=0:15
plot(kk,amp(:,1+2*j),'o-')
axis([0,N/2,-17,2])
pause(1)
figure(1)
end

Err(1) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat6(t_0,t_f,M,N,u0);
Err(2) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat6(t_0,t_f,M,N,u0);
Err(3) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat6(t_0,t_f,M,N,u0);
Err(4) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat6(t_0,t_f,M,N,u0);
Err(5) = max(abs(err(:,M+1)));
M = 2*M;
[u,err,x,t,kk,amp] = heat6(t_0,t_f,M,N,u0);
Err(6) = max(abs(err(:,M+1)));
% M = 2*M;
% [u,err,x,t,kk,amp] = heat5(t_0,t_f,M,N,u0);
% Err(7) = max(abs(err(:,M+1)))

Err
Err(1:5)./Err(2:6)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% exact solver
%
t_0 = 0; t_f = 1; M = 100; N = 32;
h = 2*pi/N;
x = 0:h:2*pi;

phases = rand(1,10)*2*pi;
amps = rand(1,10);
nn = length(amps);
u0 = zeros(size(x(1:N)));
for j=1:nn
  u0 = u0 + amps(j)*cos((j-1)*x(1:N)+phases(j));
end

[u,x,t,kk,amp] = heat7(t_0,t_f,M,N,u0);

MM = max(max(u));
mm = min(min(u));
for j=0:15
plot(x,u(:,1+2*j))
axis([0,2*pi,mm-.2,MM+.2])
pause(1)
figure(1)
end
 
for j=0:15
plot(kk,amp(:,1+2*j),'o-')
axis([0,N/2,-17,2])
pause(1)
figure(1)
end

