%% linDYStehtaviin.m
% Muunna ohjeet "moderneiksi" dokun osalta (ei enää diary-tiedostoa)
% Opettaja voi muokata tästä monenlaista hauskaa
%  
% Kurssi kp3-ii 2008  http://math.aalto.fi/opetus/kp3-ii/08/matlab/harj4tietok.m
%
% Lineaarinen 2 x 2- systeemi, diagonalisoituva matriisi
% 
%  Nimi ja opintonumero: 
%  
% Tee samaan pohjaan AV tehtävät 4,2,6 
% Kirjoita kommentin taakse muutettavat osat ja merkitse, mihin tehtävään
% mikin syöte kuuluu.
%  Palauta:
%    - tämä ajotiedosto editoituna
%    - editoitu versio diary-tiedostosta mielellään sopivin kuvin 
%      ja/tai niiden kuvauksin höystettynä kunkin tehtävän osalta.
% 
% close all
%% Määrittele matriisi A aluksi AV tehtävän 4 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;    
% Käytä muuttujia, älä Matlabin laskemia lukuja (selkeyden vuoksi ja 
% jotta voit vaivatta soveltaa samaa ajoa muihin tehtäviin).
% Y=C(1)*exp(xxx)*xxx + C(2)*exp(xxx)*xxx
%% Ratkaise vakiot C alkuarvojen avulla tyyliin:
Y0=[1;0]   % Muuteltava alkuarvo
% Y0=[25;0]  % suolatehtävässä
C= xxx\xxx    % Ratkaise yhtälöryhmä  % C=M\b, missä M=xxx ja b=xxx.
%  Selitä myös, miksi AA-tehtävän vakiot saadaan näin.
%% Nyt haluaisit piirtää ratkaisun aikakuvan ja faasikuvan.
% Mutta voi, eihän se tuosta Y:n lausekkeesta onnistu, siihen ei saa
% mitenkään järkevästi t:n paikalle t-vektoria. (Toki voitaisiin 
% muodostaa erikseen y1 ja y2 ihan kohtuullisesti, vrt L/suolafaasi.html, 
% mutta mietitään parempaa ... )
%
%   hm - hm - hm  
%
%  Mikä se juttu olikaan, joka on mainittu luennolla 27 kertaa?
%  Lineaarikombinaatiot matriisitulona, riviajattelu
%  vs. sarakeajattelu.
% (Itse asiassa, mikäli oikein ymmärsit, sitähän käytit juuri äsken
%  AA:ssa.)
% Y = X*[C(1)*exp(xxx);C(2)*exp(xxx)]  % Kokeile, saatko saman kuin yllä,
% kun t:llä on skalaariarvo.
% (Huomaa, että lienet muuttanut C:tä sen jälkeen, kun laskit ekan
% muodon, joka siis 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,30)  % AV 2 suolatehtävässä puoli tuntia, voit kokeilla muutakin.
Y = X*[C(1)*exp(xxx);C(2)*exp(xxx)];
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])  % Sama koordinaatisto kuin pplanessa.
% axis([0 25 0 25]) % Suolatehtävässä, muuta myös pplanessa.
% Miten piirretään ominaissuorat? Nehän ovat sarakkeina matriisissa X.
hold on; grid 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, saadaan ominaiskoordinaattiakselit 
      % Voit joutua muuttamaan kertoimia eri tehtävissä.

%%
% Tulin kehittäneeksi pienimuotoisen suuntakenttähommankin mukaan
% Matlabin tehokkailla työkaluilla. Mutta jätetään se nyt tähän, jos
% jotakuta sattuisi kiehtomaan. 
% Meille tässä riittäköön, kun tehdään suuntakentät pplanella.
%       
%%%%%%% THE END %%%%%%%      
%% 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\xxx         % Ratkaise vakiot.
t=linspace(-3,1,50);   % Kokeile eri arvoja
Y = X*[C(1)*exp(d(1)*t);C(2)*exp(d(2)*t)]; % Kai muistit ratkaista C:n.
plot(Y(1,:),Y(2,:),Y0(1),Y0(2),'o')    % Faasikuva ja alkuarvopiste

%%
%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  
