Heiluriyhtälö ja suuntakentät

Sovitetaan omat skriptit ja kirjan esitykset yhteen.

clear all
a=-5; b=5; c=-3; d=3;
h=0.5;
y1=a:h:b;
y2=c:h:d;

% Määritellään diffyhtsys-funktio f inline-tyylillä. Kestämme sen, että
% joudumme kirjoittamaan 2 eri koodia f:lle. Tämä on pikainen kertakäyttö-
% koodi. Varsinainen m-tiedosto on se, jota odexxx kutsuu.
%
% Koska suuntakenttiä ei juuri kannata piirtää muille kuin autonomisille
% systeemeille, suunnittelemme jatkot siten, että f on kirjoitettu ilman
% t-parametria.
%
% f=inline('[y2; -sin(y1)]','y1','y2')   % Puretaan sittenkin komponenteiksi,
%                                        % kts. alla.

%                                 % Erikseen tämä tiedostossa hei.m
%                                   %%%%%%%%%%%%%%%%%%%%%%%
f1=inline('y2','y1','y2')           %function z=hei(t,y)
f2=inline('-sin(y1)','y1','y2')     %  z=[y(2);-sin(y(1))];
                                    %%%%%%%%%%%%%%%%%%%%%%%%
[Y1,Y2]=meshgrid(y1,y2);

dy1=f1(Y1,Y2); dy2=f2(Y1,Y2);
clf
%quiver(Y1,Y2,dy1,dy2)  
normit=sqrt(dy1.^2+dy2.^2);
quiver(Y1,Y2,0.7*h*dy1./normit,0.7*h*dy2./normit,0,'c') % Väri 'cyan'
% Argumentti 0 lopussa poistaa automaattisen skaalauksen
axis equal
axis([a,b,c,d])
hold on
plot(Y1(:),Y2(:),'xr')
% Piirretään kokeeksi hilapisteet näkyviin (punaisella ristillä).

%%%%%%%%% hei.m %%%%%%%%%%%%%%%%%
%function z=hei(t,y)
%     z=[y(2);-sin(y(1))];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ya=ginput(1)


 [T,Y]=ode45('hei',[0 10],ya);
 plot(Y(:,1),Y(:,2),'g',Y(1,1),Y(1,2),'ok') 

 figure(2)  % Toiseen ikkunaan aikakuva.

 plot(T,Y)
 xlabel t
 ylabel 'y1,y2'
 title('Aikakuva')

Annetaan mennä luupissa

Aloita tästä yhden kerran

clear all

% Piirtoikkuna:
a=-5; b=5; c=-3; d=3;
h=0.5;
y1=a:h:b;
y2=c:h:d;

% Systeemi (autonominen)

f1=inline('y2','y1','y2')
f2=inline('-sin(y1)','y1','y2')


% Suuntakenttä:
[Y1,Y2]=meshgrid(y1,y2);
dy1=f1(Y1,Y2); dy2=f2(Y1,Y2);
normit=sqrt(dy1.^2+dy2.^2);
normit(normit==0)=1;   % Poistetaan mahdollinen 0:lla jako.
clf
quiver(Y1,Y2,0.7*h*dy1./normit,0.7*h*dy2./normit,0,'c')
% Argumentti 0 lopussa poistaa automaattisen skaalauksen
axis equal
axis([a,b,c,d])
hold on
plot(Y1(:),Y2(:),'xr')   

% Ratkaisukäyriä: Valitaan hiirellä, enter kuvaikkunassa lopettaa.

varit='mrgbkmr'
vari=1;

(Lisää) ratkaisukäyriä samaan suuntakenttään, jatka tästä:

ya=[0;0];
while length(ya)>0 
  [T,Y]=ode45('hei',[0,10],ya);
  plot(Y(:,1),Y(:,2),varit(vari),ya(1),ya(2),'ok') 
  vari=vari+1; vari=1+mod(vari,7);
  ya=ginput(1);
end;

shg  % Show graphics