Numsym L/ODE.html, Diffyhtälöitä, erityisesti numeriikkaa

Kirjallisuutta

[CvL] C. van Loan: Chapter 9, ODE jaossa 25.2.,
      Kirjan ODE-koodit
[Moler] Moler: Numerical Computing with Matlab
          Matlab-koodit, myös pdf-tiedostot luvuittain
[Ch-Kin] Cheney-Kincaid: Numerical Mathematics
[B-F] Burden-Faires: Numerical Analysis

Linkit aikaisempiin

• odemaple.html
• odemaple.mws
• ODEmbl1.html
• ode1.mws
• maple/ns04.mpl Maple-koodeja, teksti (read)

1. Perusteet

Diffyhtälön (alkuarvotehtävän) määritelmä:

 (AA)   y'=f(t,y), y(t_0)=y_0 

      y'(t)=f(t,y(t)) kaikilla t jollain välillä, ja joka
      toteuttaa alkuehdon: y(t_0)=y(0)

Käsittelimme luennolla heiluriyhtälön muuntamista 1. kertaluvun systeemiksi ja sen kirjoittamista Matlab-funktioksi ODE-ratkaisijaa varten tässä näin. Tiedosto hei.m

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

  heisys=inline('[y(2);-sin(y(1))]','t','y')

ya=[1 0];
[T,Y]=ode45(@hei,[0 10],ya);
% tai:
[T,Y]=ode45(heisys,[0 10],ya);

2. Yksiaskelmenetelmät

2.1. Eulerin menetelmä

function [T,Y]=eulerskalaari(fnimi,Tspan,ya,n)
% Tämä vain kehittely- ja opettelutarkoituksessa. 
% Funktio eulerV hoitaa niin skalaari- kuin vektoriversion.
%  (24.2.04)
a=Tspan(1);b=Tspan(2);
h=(b-a)/n;
Y=zeros(1,n+1);T=a:h:b;
Y(1)=ya;
for j=1:n
   Y(j+1)=Y(j)+h*feval(fnimi,T(j),Y(j));
end;

Eulerin menetelmä Matlabin ODE-tyylillä

function [T,Y]=eulerV(fnimi,Tspan,ya,n)
% fnimi : funktiomääritys, inline tai m-tiedosto, määritellään kuten
%         ode45:lle ym.  (help ode45)
%         Tspan on muotoa [a b]  ya:alkuarvovektori (y(a)=ya)
%          n: välin (a,b) jakovälien lukumäärä
% Tulos: T: sarake, aikapisteet.
%        Y: Matriisi, jonka j:s sarake ilmaisee ratkaisufunktion
%           y_j arvot T-aikapisteissä.
%           i:s rivi: T(i)  Y(i,1), Y(i,2), ..., Y(i,N)
% Siis jos on 2x2-systeemi, niin Y:ssä on 2 saraketta. 
%
% Esim: Fsys=inline('[y(1)-y(1)*y(2)/2;-2*y(2)+y(1)*y(2)]','t','y'); 
%       y0=[20;10];
%       [T,Y]=eulerV(Fsys,[0,1],y0,20); plot(T,Y)
     %  Huom! Jos Fsys määritellään m-funktioksi, on kutsussa oltava 
     %  'Fsys' tai @Fsys.
%
a=Tspan(1);b=Tspan(2);
h=(b-a)/n;
N=length(ya);
Y=zeros(n+1,N);  % i:s rivi:  Y(i,1), Y(i,2), ..., Y(i,N)
T=a:h:b;
Y(1,:)=ya';
for i=1:n
   Y(i+1,:)=Y(i,:)+h*(feval(fnimi,T(i),Y(i,:)))'; 
   % Ffunktio 'fnimi' palauttaa pystyvektorin, siksi pitää transponoida.
end;

2.2 Lokaali ja globaali virhe

Olkoon t->y(t) AA-tehtävän ratkaisufunktio, ts. y'(t)=f(t,y(t)).
Taylorin kaavan mukaan:

      y(t+h)=y(t)+h y'(t) + h2/2 y''(xi) = y(t) + h f(t,y(t))+O(h2).

ln = u(tn)-yn = h2/2 u''(xi).

[CvL]-kirjassa merkitään LTE_n.

Globaali katkaisuvirhe, g_n=y(t_n)-y_n , missä t->y(t) on AA-tehtävän y'=f(t,y), y(t₀)=y₀ tarkka ratkaisu.

CvL-tiedostot: ODE-luvun 9 m-fileet.
Ovat myös hakemistossa: /p/edu/mat-1.192/

Tässä Euleriin liittyvät.

-rw-r--r--    1 apiola   a-matlai      690 Feb  6 12:59 EulerRoundoff.m
-rw-r--r--    1 apiola   a-matlai      710 Feb  6 12:59 FixedEuler.m
-rw-r--r--    1 apiola   a-matlai      510 Feb  6 12:59 ShowEuler.m
-rw-r--r--    1 apiola   a-matlai      626 Feb  6 12:59 ShowEulerRoundoff.m
-rw-r--r--    1 apiola   a-matlai      326 Feb  6 12:59 ShowFixedEuler.m

..