H5T2vaiheiluri

Tehdään mieluummin suoraan mukaillen ohjetta:

Contents

http://math.aalto.fi/opetus/MatOhjelmistot/2014kevat/L/ODEohje.html

ODEfunktio "heiluri" voidaan määritellä myös funktiokahvana:

k=1;
c=0.1;  % Vaimennusvakio, Vaihdellaan.
c=0.2;  %
vaiheiluri=@(t,y) [y(2);-k*sin(y(1))-c*y(2)]
Tvali=[0 20]; y0=[1;1];        % Alkukulma=1, alkukulmanopeus=1.
[T,Y]=ode45(vaiheiluri,Tvali,y0); % ode45 havaitsee funktiokahvan, sille
                               % ei tarvitse sitä enää @-merkillä
                               % ilmoittaa. (Muistisääntö 1 @-merkki
                               % joko määrittelyssä tai kutsussa, ei molemmissa.)
size(T)  % Aikapistevektorin pisteiden lkm.
size(Y)  % Kaksisarakkeinen matriisi, ratkaisut y1 ja y2 aikapisteissä t.
plot(T,Y(:,1),T,Y(:,2))  % Aikakuva: heilahduskulma ja kulmanopeus ...
                         %              ajan funktiona
figure
plot(Y(:,1),Y(:,2))     % Faasikuva: Kulmanopeus kulman "funktiona" ...
                        % eli käyrä (Theta(t),omega(t)), 0< t < 20).

vaiheiluri = 

    @(t,y)[y(2);-k*sin(y(1))-c*y(2)]


ans =

   101     1


ans =

   101     2

Kriittset pisteet samat kuin vaimentamattomalla

sin(y1)=0, y2=0, siis y1-akselilla n*pi.

3 eri trajektorityypppiä:

close all
ya0=[1;1]; % Pieni alkukulma/nopeus -> heilahtelee edestakaisin
yb0plus=[-pi;0.01]; % Kriittinen alkutilanne (Heiluri pysytyasennossa+)
yb0miinus=[-pi;-0.01]; % Kriittinen alkutilanne (Heiluri pysytyasennossa-)
yc0=[-5;2];   % Suuret alkuarvot, heiluri pyörii ympäri
yd0=[5;-2];   %    ... toisinpäin
Tvali=[0 10];
[ta,Ya]=ode45(vaiheiluri,Tvali,ya0);
[tbplus,Ybplus]=ode45(vaiheiluri,Tvali,yb0plus);
[tbmiinus,Ybmiinus]=ode45(vaiheiluri,Tvali,yb0miinus);
[tc,Yc]=ode45(vaiheiluri,Tvali,yc0);
[td,Yd]=ode45(vaiheiluri,Tvali,yd0);
hold on
plot(Ya(:,1),Ya(:,2),Ybplus(:,1),Ybplus(:,2),Ybmiinus(:,1),Ybmiinus(:,2),...
    Yc(:,1),Yc(:,2),Yd(:,1),Yd(:,2))
axis equal; axis([-6 6 -4 4])
xlabel y_1(t); ylabel y_2(t)
grid on

Hiiri avuksi, samalla suuntakentta:

[y1,y2]=meshgrid(-6:0.5:6,-6:0.5:6); % Koordinaattihilapisteistö
dy1=y2; dy2=-sin(y1);  % Heiluridiffyhtälöryhmän derivaatat hilapisteissä
quiver(y1,y2,dy1,dy2)  % Suuntanuolet hilapisteisiin derivaattojen suuntaan.
% Toden totta, näillä 3:lla koodirivillä saat suuntakentän!

Tahan jatkoksi se valmis ginput-skripti pienin virityksin ...

Odottakoon toistaiseksi. (Toki voit tehdä pplanella)

Published with MATLAB® R2013a