Tietokoneharjoitus 2 ohjeita

Differentiaaliyhtälöitä Matlab:lla

Matlab:n differentiaaliyhtälörartkaisijat ovat kappale hienointa Matlab:ia. Niihin on todella panostettu (ja panostetaan), ja jälki on sen mukaista. Ratkaisumenetelmät ovat numeerisia, kuten Matlab:ssa yleensäkin.
Matlab:n dokumentaatio: help ode45 tai doc ode45

ODE-ratkaisijan peruskäyttö heiluri-esimerkillä valaistuna

ode = "Ordinary Differential Equation"

Yleisratkaisija: ode45 , kaikkien ode-tyyppisten ratkaisijoiden kutsusyntaksi on sama, perusmuodossaan tässä tapauksessa
[t,Y]=ode45(@heiluri,Tvali,y0) ,
missä heiluri määrittelee diffyhtälösysteemin (esim. heiluri alla ja esim. Tvali=[0 10]; y0=[1;1];

Tee näin

• Mene Matlab-istunnon ylänauhasta valitsemalla omaan työskentelyhakemistoosi.
• Ota Matlab:n FILE-valikosta NEW->m-file ja kirjoita (tai cut/paste):

  function dy = heiluri(t,y)
  dy=[y(2);-sin(y(1))];
  

• Talleta funktio nimelle heiluri.m työskentelyhakemistoosi,
• Testaa funktiotasi tyyliin

  heiluri(1,[1,2]), heiluri(1,[pi/2,-5]), heiluri(1,[pi,-5])
  

Jos ihmettelet, mitä virkaa argumentilla t on, eihän sille tehdä mitään, niin vastaus on: "hyvä kysymys" ja sitten:
Sen on oltava yhdenmukaisuuden vuoksi, koska ode45 vaatii , että diffyhtälösysteemin argumenttilista on aina samanlainen, esiintyypä siinä t tai ei, ts. onko systeemi autonominen tai ei.
(Jos et heti sisäistä tätä, niin älä huoli, kunhan teet, mitä käsketään!)

Pikku kokeilu

Tvali=[0 10]; y0=[1;1]; 
[t,Y]=ode45(@heiluri,Tvali,y0);
size(t)  % Aikapistevektorin pisteiden lkm.
size(Y)  % Kaksisarakkeinen matriisi, ratkaisut y1 ja y2 aikapisteissä t.
plot(t,Y(:,1),t,Y(:,2))  % Aikakuva
figure
plot(Y(:,1),Y(:,2))     % Faasikuva

Yhdistetään mukaan suuntakentän piirto

Suuntakentän piirtoa varten tarvitaan diffyhtälösysteemin määrittelemät derivaatat dy1 ja dy2. Ne annetaan 2. rivillä. Siihen voit editoida haluamasi muun 2 x 2 (autonomisen) ryhmän määrittelyn.

close all
[y1,y2]=meshgrid(-5:0.5:5,-5:0.5:5); % Koordinaattihilapisteistö
dy1=y2; dy2=-sin(y1);  % Heiluridiffyhtälöryhmän derivaatat hilapisteissä
quiver(y1,y2,dy1,dy2)  % Suuntanuolet hilapisteisiin derivaattojen suuntaan.

Nyt voit piirtää lisää trajektoreita samaan kuvaan vaikka näin:

ya0=[1;1]; yb0=[-5;2]; yc0=[5;-2]; % Kolme eri alkuarvovektoria
Tvali=[0 10]; 
[ta,Ya]=ode45(@heiluri,Tvali,ya0);
[tb,Yb]=ode45(@heiluri,Tvali,yb0);
[tc,Yc]=ode45(@heiluri,Tvali,yc0);
hold on
plot(Ya(:,1),Ya(:,2),Yb(:,1),Yb(:,2),Yc(:,1),Yc(:,2))
axis equal; axis([-5 5 -5 5])
xlabel y_1(t); ylabel y_2(t)