%----------------------------------------
% Numeerinen ja symbolinen laskenta 2005
% Harjoitus 6, tehtävä 3
%----------------------------------------

% Nämä ovat samat kuin Maplessa.
%-------------------------------
clear; close all
c=1;
L=10;
Tol=1e-6;
T=[0.2 0.1 0.05 0.025 0.01];
%TolN=[15 20 27 37 57];  
TolN=[25 34 46 63 96];

nx=100;
nt=20;
Tmax=1;
dt=Tmax/(nt-1);

x=linspace(0,L,nx);
t=linspace(0,Tmax,nt);
ZS=struct('Z',{},'ERR',{});

figure(1), clf
for i=1:5
    ZS(i).Z=lamposin(c,x,t,'sin100.txt',TolN(i));
%    ZS(i).Z=lamposin(c,x,t,'hab.txt',TolN(i));
    subplot(2,3,i)
    surf(x,t,ZS(i).Z);
end

u0=inline('-sin(pi*x/5).*exp(x/2)','x');
ga=inline('zeros(size(x))','x');
gb=inline('zeros(size(x))','x');
[cx,ct,cZ]=CrNich(u0,ga,gb,[0,L],nx-2,nt-1,dt);

figure(2), clf
for i=1:5
    ZS(i).ERR=abs(ZS(i).Z-cZ);
    subplot(2,3,i)
    surf(x,t,ZS(i).ERR);
    axis([0 L 0 Tmax 0 1])
end
figure(1)
subplot(2,3,6)
surf(x,t,cZ)

maxVirheet=zeros(5,1);
for i=1:5
    maxVirheet(i)=max(max(ZS(i).ERR));
end

[TolN(:),maxVirheet(:)]

% Yllä oleva on vähän karkeaa. Parempi olisi katsoa ao. aikasiivua tähän tapaan:

find(abs(t-0.2)==min(abs(t-0.2)))  % Haetaan t-vektorista lähinnä arvoa 0.2 
                                   % olevan indeksi ==>  5
U=ZS(1).Z;
figure
plot(x,U(5,:))   Analyyttinen ratkaisu "5:nnellä ajanhetkellä"
% Vastaavasti numeerisen ratkaisun ja analyyttisen vertailu.
virhe=ZS(1).ERR;
plot(x,virhe(5,:))
max(virhe(5,:))
% Muissa aikapisteissä vastaavasti.

