% matlab/harj4tietok.m
%  Nimi ja opintonumero: 
%  
% Lineaarinen 2 x 2- systeemi, diagonalisoituva matriisi
% close all
%% Määrittele matriisi A tehtävän harj 4 AV mukaisesti.
%% Muodosta matriisi X, jonka sarakkeina ovat ominaisvektorit ja 
% matriisi D, jonka lävistäjällä vastaavat ominaisarvot.
%% Muodosta vektori d, jossa on D:n lävistäjä.
%% Kirjoita yleinen ratkaisukaava tähän tyyliin:
% C(1)=1; C(2) = 1; % sama kuin C=[1 1];  
% t=0.2;    
% Y=C(1)*exp(xxx)*xxx + C(2)*exp(xxx)*xxx
%% Ratkaise vakiot C alkuarvojen avulla tyyliin:
Y0=[1;0]   % Muuteltava alkuarvo
C= xxx     % Ratkaise yhtälöryhmä  % C=M\b, missä M=xxx ja b=xxx.
%% Nyt haluaisit piirtää ratkaisun aikakuvan ja ennenkaikkea faasikuvan.
% Mutta voi, eihän se tuosta Y:n lausekkeesta onnistu, siihen ei saa
% mitenkään järkevästi t:n paikalle t-vektoria.
%   hm - hm - hm  
%
%  Mikä se juttu olikaan, joka on mainittu luennolla 27 kertaa?
%  Lineaarikombinaatiot matriisitulona, riviajattelu vs. sarakeajattelu:
% Y = X*[C(1)*exp(xxx);C(2)*exp(xxx)]  % Kokeile, saatko saman kuin yllä,
%    kun t:llä on skalaariarvo.
% Jälkimmäisen muodon suuri etu on tässä se, että kerrottavia vektoreita
% voi olla paljon, ne kootaan 2-riviseksi matriisiksi, ja homma hoituu yhdellä
% matriisitulolla.
%% No niin:
Y0=[1;1]; C=X\xxx;
t=linspace(-3,0.5,50);
Y = X*[C(1)*exp(xxx);C(2)*exp(xxx)];
Y(:,1:10)
 plot(t,Y(1,:),t,Y(2,:)); grid on; shg  % Aikakuva
 figure   % seuraava ikkuna
 plot(Y(1,:),Y(2,:))  % Faasikuva
% Miten piirretään ominaisvektorit? Nehän ovat sarakkeina matriisissa X.
hold on
plot([0 X(1,1)],[0,X(2,1)],[0,X(1,2)],[0,X(2,2)])

%% Olispa hauskaa, jos voitaisiin piirtää suuntakenttä ja siinä valita
%  alkupiste hiirellä, kuten Java-pplane:ssa.
%  No näin vähäisellä koodilla se käy tässä:
[y1,y2]=meshgrid(-5:0.5:5,-5:0.5:5); % Koordinaattihilapisteistö
dy1=A(1,1)*y1+A(1,2)*y2; dy2=A(2,1)*y1+A(2,2)*y2;  % Derivaattavektori
% saadaan lineaariselle näin: dy=A*y, jos y on jokin 2 x 1-vektori. 
% Tässä on kirjoitettava matriisitulo auki, koska y1 ja y2 ovat
% hilapisteistön määrääviä matriiseja, joten dy=A*[y1;y2] on "mahdoton yhtälö". 
quiver(y1,y2,dy1,dy2) % Suuntanuolet hilapisteisiin derivaattojen suuntaan.
axis([-5 5 -5 5]); grid on;hold on
Y0=ginput(1)'; % Matlab antaa hiusristikon, jonka avulla voit hiirellä valita
                % alkupisteen.
C=xxx\Y0        % Ratkaise vakiot.
t=linspace(-3,1,50);   % Kokeile eri arvoja
Y = X*[C(1)*exp(xxx);C(2)*exp(xxx)];
plot(Y(1,:),Y(2,:),Y0(1),Y0(2),'o')    % Faasikuva ja alkuarvopiste
%axis([-5 5 -5 5]); grid on; shg

%%
% Tästä eteenpäin olevia kommentteja EI SAA POISTAA!
%
% Ajatellaan, että olet saanut homman toimimaan ja olet tallettanut sen
% tiedostoon linsys1.m
% Voit ajaa Matlab-istunnossa komentamalla  linsys1, 
% mutta ota yllä ja alla olevat kiellot kiltisti onkeesi.
% Alla lisäehdotuksia, jos haluat piirtää useita trajektoreita.
%   
%% HUOM! HUOM! HUOM! 
%  Oletetaan, että yllä skripti on tiedostossa  linsys1.m
%       ÄLÄ IHMEESSÄ poista seuraavia kommentteja ja suorita
%       komentoa linsys1  ---  silloin skripti kutsuu ikuisesti itseään ---
%       Silloin paina CTR-C  ja toivo parasta!
%% Matlab-istunnossa voit komentaa tähän tapaan, kun haluat piirtää
%     3 tai 8 trajektoria:
% 
%  close all  % Suljetaan grafiikkaikkunat
%  linsys1;linsys1;linsys1
%  Tai esim:
% for i=1:8
%   linsys1
% end

