Differentiaaliyhtälö(systeemi)t

ODE = Ordinary Differential Equation

MATLABissa on tavallisten kattava valikoima 1. kertaluvun differentiaaliyhtälöryhmien alkuarvotehtävien ratkaisufunktioita. Tässä muutama: ode45, ode23 ja ode15s. Korkeampaa kertalukua olevat tavalliset dy:t on muutettavan 1. kl:n muotoon ennen ratkaisua

Kaikkien kutsusyntaksi on sama.

Ensin vähän numeerisen ratkaisemisen teoriaa ja Matlab-esimerkkejä.

Eulerin menetelmä

Tästä tarina aina alkaa, siitä voit lukaista vaikka yllä olevasta luentoprujulinkistä. Matlab-toteutus voidaan yksinkertaisimmillaan tehdä tällä funktiolla:
function [T,Y]=eulerS(f,Tspan,ya,n)
% Tämä vain kehittely- ja opettelu- ja havainnollistustarkoituksessa. 
% Funktio eulerV hoitaa niin skalaari- kuin vektoriversion.
%  HA  21.8.2010
% Esim: y'=t+y, y(0)=1
%       f=@(t,y)t+y
%       [T,Y]=eulerS(f,[0 4],1,6), plot(T,Y,T,Y,'.r');shg  
a=Tspan(1);b=Tspan(2);
h=(b-a)/n;
Y=zeros(n+1,1);T=(a:h:b)'; %Pystyvektorit yhdenmukaisesti ode45:n
Y(1)=ya;                   %  kanssa 
for j=1:n
   Y(j+1)=Y(j)+h*f(T(j),Y(j));
end;
Funktion argumenttilista ja tulosten muoto on järjestetty juuri samaksi kuin Matlabin omien ratkaisijoiden. (Paitsi tässä tuo havainnollistukseen käytettävä viimeinen argumentti n.)

Menetelmän toimintaan voit omakohtaisesti tutustua tämän skriptin avulla. Huomaa, että sitä ei ajeta kokonaisena skriptinä, vaan leikataan/liimataan osissa ja nn-iteroidaan. (nn=nuolinäppäin)

% lyhyt/euleresim1.m
% Eulerin menetelmä: alkuarvothtävälle y'=t+y, y(0)=1
% Nuolinäppäinversio, havainnollistus
% Lineaarinen epähomogeeninen, tarkka ratkaisu: 
%              y(t)=2exp(t) - t -1
clear; close all
% y'=t+y, y(0)=1, ratk. välillä [0,1.2]
f=@(t,y)t+y
YL=[]   % Talletetaan lasketut arvot välin loppupisteessä
n=3;
[T,Y]=eulerS(f,[0 4],1,n); [T Y],plot(T,Y,T,Y,'.r');shg;YL=[YL Y(end)]
grid on; hold on
%% Suorita tätä riviä nuolinäppäiniteroiden (tai editorissa: CTR-ENTER):
n=2*n,[T,Y]=eulerS(f,[0 4],1,n); [T Y];plot(T,Y,T,Y,'.r');shg;YL=[YL Y(end)]
%% Iteroinnin jälkeen:
ytarkka=@(t) 2*exp(t)-t-1;
t=linspace(0,4);
plot(t,ytarkka(t),'k','LineWidth',3);shg
%fplot(ytarkka,[0 4])  %plot-optiot eivät kaikki toimi fplot:lla                    
format compact
['Loppupiste Tarkka ratk.  Euler-approksimaatio']
[T(end) ytarkka(T(end)) YL(end)]
virheetloppupisteesa=ytarkka(T(end))-YL
Tehtävä
  1. Tutustu ja kokeile skriptin ajamista. Muuttele parametreja, voit muuttaa myös diffyhtälöä.
  2. Nuolinäppäinversio on helppo muuttaa while 1-rakenteiseksi ikuiseksi silmukaksi, joka katkeaa käyttäjän päätöksellä. Kts. esim/euleresim1b.m. Tutustu ja kokeile tätäkin.
Edellinen on joustavampi, jälkimmäinen on kokonainen skripti, joka voidaan ajaa sellaisenaan ilman leikkau/liimaus-vaihetta.
Itse kallistun tässä yhteydessä (lievästi) 1. vaihtoehdon puolelle.

Eulerskriptin ajo html-kauniiksi

euleresim1.m html-julkaisuksi ajettuna ("publish").

Runge-Kutta muunnelmineen

Eulerin menetelmä on hyödyllinen ennen kaikkea ideatasolla, se kun on riisuttu ylimääräisistä konstailuista ja tarjoaa yksinkertaisen ajatusmallin edistyneempien pohjaksi. Matlab-käyttäjän kannattaa yleensä ensialkuun käydä ratkaisufunktion ode45 kimppuun, ellei ongelman erityisluonteesta (huomattava tarkkuusvaatimus, "kankeus", ym.) ole tarkempaa tietoa.

Yllä oleva tehtävä voitaisiin ratkaista näin:

   close all; clear
   f=@(t,y)t+y
   tspan=[0 4];
   y0=1;
   [T,Y]=ode45(f,tspan,y0);
   plot(T,Y,'*--')
   xlabel t, ylabel y(t)
Jatketaan katsomalla tarkemmin tulosvektoreita T ja Y sekä virhettä.
   size(T)
   diff(T)  % Vektori (T(k+1) - T(k)), nähdään, että askeleet alussa ja lopussa
            % n. puolet keskivaiheen askelista.
   ytarkka=@(t) 2*exp(t)-t-1;
   subplot(2,1,1)
     plot(T,Y)   
     title('Runge-Kutta 45-ratkaisu: y''=t+y, y(0)=1') 
                   % (Merkkijonossa y' hipsukan kahdennus.)
   subplot(2,1,2)
      plot(T,ytarkka(T)-Y)
      title('Virhe: ytarkka(t)-y_{RK}(t)')
   shg   
Komento diff(T) ilmaisee, että ratkaisija ode45 valitsee diskretointipisteet omin päin. Tässä alku- ja loppupään askelpituudet ovat n. puolet keskiosan askelpituudesta. Kyseessä on ns. adaptiivinen askelpituuden valinta, joka tehdään lokaalien virhearvioiden perusteella. (Kts. [Moler], [CvL]) Tämä näkyy niin, että askeleet ovat yleensä pienempiä alueella, jolla ratkaisufunktio kasvaa voimakkaammin.

Entä, jos halutaan laskea arvoja muissa kuin ode45:n valitsemissa pisteissä.

  1. Annetaan argumentiksi välin (kuten [0 4]) sijasta pisteistö, jossa tulokset halutaan.
  2. Suoritetaan splinisovitus.
Katsotaan 1. vaihtoehtoa. Tällöin voitaisiin kutsua näin:
    tpisteet=linspace(0,4);
    [Y,T]=ode45(f,tpisteet,y0);
Tarkista, että ode45 palauttaa todella samat, syötteenä annetut tpisteet.
    all(T==tpisteet')
    sum(T~=tpisteet')   % vaihtoehtoisesti
Tämä ei merkitse, että annettuja tpisteet-arvoja käytettäisiin laskennan jakopisteinä.

Suuntakenttä

Kts. [Hig-Hig] s. 176-177

Korkeamman kertaluvun yhtälöt, diffyhtälösysteemit

Matlabin ODE-ratkaisijat käsittelevät 1. kertaluvun systeemejä. Tämä ei ole mikään rajoitus, sillä korkeamman kertaluvun yhtälöt voidaan muuntaa helpolla muunnoksella 1. kertaluvun systeemiksi. Jos yhtälö on vaikka 3. kertaluvun:
y''' = F(t,y,y',y''),
merkitään: y=y1, y'=y2, y''=y3, jolloin yhtälö muuntuu systeemiksi:
       y1' = y2
       y2' = y3
       y3' = F(t,y1,y2,y3)
y=y1

Luentokalvolla no 14 käsitellään heiluriyhtälöä

  
                y'' = -k sin(y)
Tässä y(t) on heilahduskulma, y'(t) kulmanopeus ja y''(t) kulmakiihtyvyys.
[Edellinen] [Seuraava] [Alkusivu]