Differentiaaliyhtälö(systeemi)t

• Luentoja diffyhtälöiden numeriikasta (KP3/2008)
• Moler odes
• help -> Mathematics -> ODEs -> Initial value problems
• Ohje Matlabin ODE-rutiinien käyttötekniikkaan (KP3/2008 HA)

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ä

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;

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

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.

Eulerskriptin ajo html-kauniiksi

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

Runge-Kutta muunnelmineen

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)

   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   

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.

    tpisteet=linspace(0,4);
    [Y,T]=ode45(f,tpisteet,y0);

    all(T==tpisteet')
    sum(T~=tpisteet')   % vaihtoehtoisesti

Suuntakenttä

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

y''' = F(t,y,y',y''),

       y1' = y2
       y2' = y3
       y3' = F(t,y1,y2,y3)
y=y1

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

  
                y'' = -k sin(y)