% This file solves the Forward Heat Conduction Problem
%
% u_t = u_xx,  0 < x < pi (+ boundary & initial conditions)
% 
% for (quite different) initial temperature distributions.
% 
%
% Jennifer Mueller, Samuli Siltanen and Esa Niemi.


% Choose relative amplitude of noise (0.01 equals 1 percent noise)
noiselevel = 0.1;

% Numerical parameters
L    = pi;
D    = 1;

Tmax = 1;    % use xdim = 50, tdim = 500
% Tmax = .1;   % use xdim = 20, tdim = 50
% Tmax = .2;   % use xdim = 40, tdim = 100
% Tmax = .4;  % use xdim = 24, tdim =108
%Tmax = .1;   % use xdim = 48, tdim = 216

xdim = 50;
tdim = 500;


x    = linspace(0,L,xdim);
t    = linspace(0,Tmax,tdim);
dt   = t(2)-t(1);
dx   = x(2)-x(1);
% check the stability criterion dt/(dx^2)<.5
if (dt/(dx^2)>.5)
    display('Stability criterion violated')
    display('dt/(dx^2) = ')
    dt/(dx^2)
end


% Synthesize data by time stepping up to Tmax. Rows represent a fixed
% time, columns represent spatial location
ZM1 = zeros(xdim,tdim);   % allocate space for the solution as a matrix
ZM2 = zeros(xdim,tdim);
ZM1(:,1) = 10*sin(2*x);   % the initial temperature
ZM2(:,1) = 12*nonDiffSin2(x);   % another initial temperature
% ZM2(:,1) = ZM1(:,1)+noiselevel*max(abs(ZM1(:,1)))*randn(size(ZM1(:,1)));   % noisy initial temperature

ZM1(1,:) = 0;             % the BC u(0,t)
ZM1(xdim,:) = 0;          % the BC u(L,t)
ZM2(1,:) = 0;             % the BC u(0,t)
ZM2(xdim,:) = 0;          % the BC u(L,t)

mmax = max(max(ZM1(:,1)),max(ZM2(:,1)));

figure(1)
clf
plot(x,ZM1(:,1),'b', x,ZM2(:,1),'r')
legend('temp. distr. 1', 'temp. distr. 2')
axis([0, L, -mmax, mmax])
title('Two initial temperature distributions')
pause

s = D*(dt/dx^2);
n=length(x) - 2;
A = (1-2*s)*eye(n,n);
for ii =1:n-1
    A(ii,ii+1)=s;
end
for ii = 1:n-1
    A(ii+1,ii)=s;
end
v=2:(xdim-1);
for tin=2:tdim
	tip = tin-1;
	ZM1(v,tin) =  A*ZM1(v,tip);  % If BCs are not zero, need to  add a column vector to the result
    ZM2(v,tin) =  A*ZM2(v,tip);
    if (mod(tin,round(tdim/10)) == 0) || (tin==tdim)
        figure(1)
        clf
        plot(x,ZM1(:,tin),'b', x,ZM2(:,tin),'r')
        legend('temp. distr. 1', 'temp. distr. 2')
        title(['Temperature distributions at time t=', num2str(t(tin))])
        axis([0, L, -mmax, mmax])
        pause
    end
end


