% matlab/linsys1.m
% Lineaarinen systeemi, versio 1, Vrt. tehtävä av 4
% Diagonalisoituva matriisi
%close all
A=[1/3 4/3;8/3 5/3];  % Tähän voit vaihtaa matriisin.
[X,D]=eig(A);
d=diag(D);
%C1=0.1; C2 = .1;    % Näitä muuttelemalla voitaisiin piirrellä.
% Mutta näin on hallitumpaa.
Y0=[1;0]   % Muuttele tätä
C=X\Y0;    %%% Selitä, mitä tässä tapahtuu.

% Tekisi ehkä mieli kirjoittaa näin:
%t=1;
%Y=C(1)*exp(d(1)*t)*X(:,1)+C(2)*exp(d(2)*t)*X(:,2)
% No, mikä ettei, mutta entä jos t on vektori, kuten yleensä.
% Mitä onkaan opittu lineaarikombinaatoista ja matriisitulosta:
%Y = X*[C(1)*exp(d(1)*t);C(2)*exp(d(2)*t)]  % Kokeile, saatko saman!
%% Kuvaa lyhyesti sanallisesti, mikä matriisilaskentasääntö on taas
%  käytössä.
% Jälkimmäisen muodon suuri etu on tässä se, että kerrottavat vektorit
% voidaan koota 2-riviseksi matriisiksi ja homma hoituu yhdellä
% matriisitulolla.
%% Toista edellinen tähän tapaan:
t=linspace(-3,0.5,50);
Y = X*[C(1)*exp(d(1)*t);C(2)*exp(d(2)*t)];
Y(:,1:10)
%plot(Y(1,:),Y(2,:))
%figure   % seuraava ikkuna
%plot(t,Y(1,:),t,Y(2,:)); grid on; shg  % Aikakuva
%% 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 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)';
C=X\Y0
t=linspace(-3,1,50);   % Kokeile eri arvoja
Y = X*[C(1)*exp(d(1)*t);C(2)*exp(d(2)*t)];
plot(Y(1,:),Y(2,:),Y0(1),Y0(2),'o')    % Faasikuva ja alkuarvopiste
%axis([-5 5 -5 5]); grid on; shg
% Kokeile: MUTTA ÄLÄ IHMEESSÄ poista seuraavia kommentteja ja suorita
%   komentoa linsys  ---  silloin mennään ikuisen silmukkaan --- ja savu
%   nousee!
% linsys;linsys;linsys
%  Tai esim:
% for i=1:8
%   linsys
% end

