% ratk4tietok.m   
% Lineaarinen 2 x 2- systeemi, diagonalisoituva matriisi
% close all
format compact
%% AV 4 
A=[1/3 4/3;8/3 5/3]
%% AV 2
A=[-1/10 1/25;1/10 -1/10]
%% AV 6
A=[-1/2 -2; 1 3/2]
%% Muodosta matriisi X, jonka sarakkeina ovat ominaisvektorit ja 
% matriisi D, jonka lävistäjällä vastaavat ominaisarvot.
[X,D]=eig(A)
%% Muodosta vektori d, jossa on D:n lävistäjä.
d=diag(D)
%% 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(d(1)*t)*X(:,1) + C(2)*exp(d(2)*t)*X(:,2)
%% Ratkaise vakiot C alkuarvojen avulla tyyliin:
Y0=[1;0]   % AV 4 % Muuteltava alkuarvo
%Y0=[25;0]  % AV 2
Y0=[1;1];
C= X\Y0    % 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.
% (Itse asiassa käytit sitä juuri AA-kertoimissa, mikäli oikein ymmärsit.)
% Y = X*[C(1)*exp(d(1)*t);C(2)*exp(d(2)*t)]  % Kokeile, saatko saman kuin 
% yllä, kun t:llä on edellä annettu skalaariarvo.
% (Huomaa, että lienet muuttanut C:tä, joten edellinen lauseke pitää laskea
% uudestaan (cut/paste)).
% 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.
%% Nyt tositoimiin!

%t=linspace(-3,1);          % AV 4
%t=linspace(0,60);         % AV 2
t=linspace(-6*pi,pi);      % AV 6
Y = X*[C(1)*exp(d(1)*t);C(2)*exp(d(2)*t)];
Y(:,1:10)   % Katsotaan 10 ekaa pistettä.
 plot(t,Y(1,:),t,Y(2,:)); grid on; shg  % Aikakuva
 figure   % seuraava ikkuna
 plot(Y(1,:),Y(2,:))  % Faasikuva
 axis([-5 5 -5 5]) ; grid on % Sama koordinaatisto kuin pplanessa.
% Miten piirretään ominaissuorat? Nehän ovat sarakkeina matriisissa X.
hold on
plot(5*[0 X(1,1)],5*[0,X(2,1)],'k',5*[0,X(1,2)],5*[0,X(2,2)],'k')
plot(-5*[0 X(1,1)],-5*[0,X(2,1)],-5*[0,X(1,2)],-5*[0,X(2,2)]) % Myös vastavektorit

%% AV 6
a=real(X(:,1)), b=imag(X(:,1))
plot(7*[0 a(1)],7*[0,a(2)],'k',7*[0,b(1)],7*[0,b(2)],'k')
plot(-7*[0 a(1)],-7*[0,a(2)],'k',-7*[0,b(1)],-7*[0,b(2)],'k')

% Trajektorin koordinaatit vektorien a ja b määräämässä kannassa 
% ovat spiraalilla, jonka määrää kiertomatriisi kerrottuna
% exp-funktiolla. Kuvassa niiden merkitys ei valaistu yhtä selvästi
% kuin reaalisen tapauksen ominaisvektorien.

%%
% Tämä riittäköön, alla oleva jatko on "vapaaehtoista"
%
%% 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ö". 
% clf  % Puhdistetaan ikkuna 2.
quiver(y1,y2,dy1,dy2) % Suuntanuolet hilapisteisiin derivaattojen suuntaan.
axis([-5 5 -5 5]); grid on;hold on
plot(5*[0 X(1,1)],5*[0,X(2,1)],'k',5*[0,X(1,2)],5*[0,X(2,2)],'k')
plot(-5*[0 X(1,1)],-5*[0,X(2,1)],'k',-5*[0,X(1,2)],-5*[0,X(2,2)],'k')
disp('Kuljeta hiirtä grafiikkaikkunaan ja klikkaa')
Y0=ginput(1)'; % Matlab antaa hiusristikon, jonka avulla voit hiirellä valita
                % alkupisteen.
C=X\Y0;        % Ratkaise vakiot.
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

%%
%Talletetaan viimeinen kappale tiedostoon ratk4ajo.m
%Voit ajaa komennot komentamalla Matlab-ikkunassa:

close all    % grafiikkaruutujen poisto
%ratk4ajo
% Usempia trajektoreita:
%close all
%ratk4ajo;ratk4ajo;ratk4ajo
% tai:
close all
n=8;
for i=1:n
    disp(['Sinulla on ',num2str(n-i+1),' klikkaus(ta)'])
    ratk4ajo
 end  
