dynumer.html

Diff. yht. numeriikkaa

K3/p3 22.11.05

Euler, Heun, eli parannettu Euler (RK2)

Esim 1 y' = t+y

> restart:

> with(plots):

Warning, the name changecoords has been redefined

> f:=(t,y)->t+y; # Diff. yhtälön määrittelevä funktio f(t,y)

> h:=0.4:T[0]:=0:y[0]:=0: N:=3: # Alkuasetukset.

> T[1]:=T[0]+h;
y[1]:=y[0]+h*f(T[0],y[0]); # 1. Euler-askel

> T[2]:=T[1]+h;
y[2]:=y[1]+h*f(T[1],y[1]); # 2. Euler-askel

Nyt riitti manuaalinen toiminta. Aloitetaan alusta ja rakennetaan for-silmukka.

> h:=0.4:T[0]:=0:y[0]:=0: N:=3: # Alkuasetukset uudestaan.

>    for k from 0 to N do
  T[k+1]:=T[k]+h;
  y[k+1]:=y[k]+h*f(T[k],y[k])
end do:

> polku1:=[seq([T[k],y[k]],k=0..N)]; # Tässä ovat pisteet [t[k],y[k]]:

> plot(polku1,axes =box,legend="h=0.4");

> Ekuva1:=%:

> h:=h/2;T[0]:=0;y[0]:=0: N:=2*N; # Puolitetaan askel h ja kaksinkertaistetaan N.

>    for k from 0 to N do
  T[k+1]:=T[k]+h:
  y[k+1]:=y[k]+h*f(T[k],y[k])
end do:

> polku2:=[seq([T[k],y[k]],k=0..N)]: plot(polku2,color=green,axes=box,legend="h=0.2"); Ekuva2:=%:

> h:=h/2;T[0]:=0;y[0]:=0: N:=2*N; T1:=N*h; # Ja vielä kerran.

>    for k from 0 to N do
  T[k+1]:=T[k]+h:
  y[k+1]:=y[k]+h*f(T[k],y[k])
end do:

> polku3:=[seq([T[k],y[k]],k=0..N)]; plot(polku3,color=blue,axes=box,legend="h=0.1"); Ekuva3:=%:

> display(Ekuva1,Ekuva2,Ekuva3,title="Eulerin polut");

Katsotaan kunkin polun viimeinen y-arvo, eli ratkaisuapproksimaatio pisteessä t=1.2:

> yLoppu[1]:=polku1[-1][2];yLoppu[2]:=polku2[-1][2];yLoppu[3]:=polku3[-1][2];

>

> dsolve({diff(Y(t),t)=t+Y(t),Y(0)=0},Y(t));

> ytarkka:=rhs(%);

> display(Ekuva1,Ekuva2,Ekuva3,plot(ytarkka,t=0..1.2,color=black,legend="Tarkka ratkaisu"));

> ytarkka12:=eval(subs(t=1.2,ytarkka));

> virheet:=seq(ytarkka12-yLoppu[k],k=1..3);

Aika lähellä on O(h)-käytös, eli virhe puolittuu, kun askel puolittuu.

Esim. 2 y'=t+y Parannetulla E:lla

> #restart:

> with(plots):

> f:=(t,y)->t+y;

> h:=0.4;T[0]:=0:y[0]:=0: N:=3:

>    for n from 0 to N-1 do
  s1:=f(T[n],y[n]);
  s2:=f(T[n]+h,y[n]+h*s1);
  y[n+1]:=y[n]+h*(s1+s2)/2;
  T[n+1]:=T[n]+h;
end do;

> polku1:=[seq([T[k],y[k]],k=0..N)];

> plot(polku1,color=blue,legend="Heun, h=0.4");

> Hkuva1:=%:

> display(Ekuva1,Hkuva1,title="Punainen Euler, sininen Heun");

Eulerin ja Heunin menetelmät systeemille

Esim 1

> restart: with(plots): with(LinearAlgebra):

Warning, the name changecoords has been redefined

> diff(y(t),t,t)+y(t)=0;

> A:=<<0,-1>|<1,0>>; Id:=<<1,0>|<0,1>>;

> Y[0]:=<-1,0>;

> h:=0.1: N:=ceil(2*Pi/h);

> for k from 0 to N do
Y[k+1]:=(Id+h*A).Y[k];
end do:

> polku1:=[seq(convert(Y[k],list),k=0..N)]: polku1[1..4];

>

> display(plot(polku1,color=blue),complexplot(exp(I*t),t=0..2*Pi),scaling=constrained);

>

> h:=h/2: N:=ceil(2*Pi/h);

> for k from 0 to N do
Y[k+1]:=(Id+h*A).Y[k];
end do:

> polku2:=[seq(convert(Y[k],list),k=0..N)]:

>

> display(plot(polku2,color=blue),complexplot(exp(I*t),t=0..2*Pi),scaling=constrained);

On luonnollista, että Euler on huono ympyrälle ja yleisemmin keskus-tyypille (kuten kettu-jänis). Joka askeleella virhe menee samaan suuntaan,

kun itsepintaisesti mennään tangentin suuntaan. Ajaudutaan väkisin ulospäin menevälle spiraalille.

Sama Heunin menetelmällä

> h:=0.4: N:=ceil(2*Pi/h);

> for n from 0 to N do
Ytahti:=(Id+h*A).Y[n];
K:=1/2*A.(Y[n]+Ytahti);
Y[n+1]:=Y[n]+h*K;
end do:

> polku1:=[seq(convert(Y[k],list),k=0..N)]:

>

> display(plot(polku1,color=blue),complexplot(exp(I*t),t=0..2*Pi),scaling=constrained);

> h:=0.2: N:=ceil(2*Pi/h);

> for n from 0 to N do
Ytahti:=(Id+h*A).Y[n];
K:=1/2*A.(Y[n]+Ytahti);
Y[n+1]:=Y[n]+h*K;
end do:

> polku1:=[seq(convert(Y[k],list),k=0..N)]:

>

> display(plot(polku1,color=blue),complexplot(exp(I*t),t=0..2*Pi),scaling=constrained);

>

Askeleella h=0.2 saadaan jo kuvan tarkkuudella miltei tarkka. Huomattava parannus. Mitä sitten tuokaan Runge-Kutta? (Harj. teht.)