%% Harj. 2, teht. 15 Doku
% Kirjoitetaan differentiaaliyht.systeemin koodi tiedostoon diffsyst.m
%{
function dy  = diffsyst( t,x )
dy=[    -10*x(1)+28*x(2);       ...
        10*x(1)-x(2)-x(1)*x(3); ...
        -8/3*x(3)+x(1)*x(2)     ...
   ];
%}
% Huomaa, ett"a muuttuja t tarvitaan, vaikka sit"a ei tassa (autonomisen 
% systeemin tapauksessa k"a"aytet"a mihink"a"an. (Kutsusyntaksin oltava sama,
% jotta kutsuva ode-rutiini voi toimia kaikenlaisille yht.systeemeille.)
%% Ratkaistaan ode45-ratkaisijalla.
% Alkuarvot:
xyz0=[-10; -4; 12];
[T,XYZ]=ode45(@diffsyst,0:.01:40,xyz0);
% 
% Kun kutsussa annetaan T-valin sijasta T-vektori, niin tulokset
% palautetaan juuri annetun T-vektorin pisteissa. ode45 palauttaa
% T-vektorin sellaisenaan ja XYZ on 3-sarakkeinen matriisi, jonka kukin
% sarake ilmaisee ao. ratkaisufunkiton (approksimaation) ko. T-pisteessa.

%% Ratkaisufunktioiden piirto:
% Ratkaisufunktiot ajan funktiona. 

% Jaetaan grafiikkaruutu osiin subplot:lla.
subplot(3,1,1)
plot(T,XYZ(:,1))
title('x(t)')
grid on
subplot(3,1,2)
plot(T,XYZ(:,2))
title('y(t)')
grid on
subplot(3,1,3)
plot(T,XYZ(:,3))
title('z(t)')
grid on
shg

%% Faasipiirros eli k"ayr"a (x(t),y(t),z(t))
figure
plot3(XYZ(:,1),XYZ(:,2),XYZ(:,3))
title('Faasikuva')
grid on
%% Tutkitaan alkuarvojen pienen muutoksen vaikutusta ratkaisuihin.
%  Muutellaan alkuarvoja lis"a"am"all"a 0.1 alkup. vektoriin
%  kullakin kierroksella. 
%  Kolmen sarakkeen tulosvektorit on k"atev"a"a tallettaa
%  solutaulukkoon ("cell array"). Se voidaan ajatella 
%  indeksoiduksi muuttujaksi, jota indeksoidaan 
%  aaltosuluissa olevalla indeksilla.
N=2;

for k=0:N
  xyz0=(1+k*0.1)*ones(size(xyz0));
 [T,XYZ]=ode45(@diffsyst,0:.01:40,xyz0);
 XYZM{k+1}=XYZ;
end 


%% Verrataan pienen muutoksen vaikutusta ratkaisuun.
maxarvo=max(abs(XYZM{1}(:,1)))               % k=0, x(t)
erotuksenmax=max(abs(XYZM{1}(:,1)-XYZM{2}(:,1)))  % max |x(t)(k=0) - x(t)(k=1)|
subplot(2,1,1)
  plot(T,XYZM{1}(:,1))
title('x(t), kun k=0')
  subplot(2,1,2)
  plot(T,XYZM{2}(:,1))
title('x(t), kun k=1')  
shg  
%%  Ylla lasketun erotusfunktion kuva:
figure
plot(T,XYZM{2}(:,1)-XYZM{1}(:,1))
title('Erotusfunktio (x(t),k=1)-(x(t),k=0)')
shg
%%  Piirret"a"an kuvat viimeksi saaduilla arvoilla (k=2).
figure
subplot(3,1,1)
plot(T,XYZ(:,1))
title('k=2,x(t)')
subplot(3,1,2)
plot(T,XYZ(:,2))
title('k=2,y(t)')
subplot(3,1,3)
plot(T,XYZ(:,3))
title('k=2,z(t)')
shg
%% 
figure
plot3(XYZ(:,1),XYZ(:,2),XYZ(:,3))
title('k=2, faasikuva')
box;grid on
%% Mathworks "viimeisen paalle"-animaatio komennolla lorenz
% (Oma vaatimattomampi animaatio on suht. helppo tehd"a, mutta antaapa
% nyt olla.)
% Kokeile, lue doc lorenz ja Molerin NCM-kirjan sivuilta 
% 202-203 "Lorenz attractor"
%% Lue tarkemmin aiheesta: 
% 
%%
% <http://www.mathworks.se/moler/odes.pdf Moler NCM 7.8. Lorenz Attractor>

