Metódy Runge-Kutta

Metódy Runge-Kutta sú jednokrokové metódy vyšších rádov, ktoré je možné odvodiť podobným spôsobom ako Eulerovu metódu pre riešenie začiatočných úloh y'=f(t,y),y((Subscript[t, 0]))=(Subscript[y, 0]). Počas výpočtu nie je potrebné počítať  derivácie vyšších rádov ako napríklad pri Taylorových metódach, čo výrazným spôsobom umožňuje urýchliť výpočet rovnice. 
Z teórie vieme, že všeobecný výraz pre explicitné metódy Runge-Kutta je

(Subscript[y, n+1])=(Subscript[y, n])+h (∑_(i=1)^s) (Subscript[b, i])(Subscript[k, i])

(Subscript[k, i])= f ((Subscript[t, n])+(Subscript[c, i])h,(Subscript[y, n])+h(∑_(j=1)^s) (Subscript[a, ij])(Subscript[k, j]))

kde (Subscript[a, ij])=0 toj≥i and (∑_(j=1)^s) (Subscript[a, ij])=(Subscript[c, i]).
	 Klasická metóda Runge-Kutta druhého rádu s dvoma stupňami má nasledujúcu Butcherovu tabuľku:

({{0, 0, .}, {(Subscript[c, 2]), (Subscript[c, 2]), 0}, {., (Subscript[b, 1]), (Subscript[b, 2])}})

kde koeficienty diagramu overuje systém rovníc:

(Subscript[b, 1])+(Subscript[b, 2])=1
(Subscript[b, 2])(Subscript[c, 2])= (1/2)

Existuje preto nekonečne veľa metód Runge-Kutta druhého rádu. Najčastejšie sa používajú:

a) Eulerova modifikovaná metóda, ktorá zodpovedá voľbe koeficientov  (Subscript[b, 1])=0,(Subscript[b, 2])=1,(Subscript[c, 2])=(1/2). Jej iteračná schéma má tvar (Subscript[t, n+1])=(Subscript[t, n])+h,  (Subscript[y, n+1])=(Subscript[y, n])+ h (Subscript[k, 2]), kde   (Subscript[k, 1])=f((Subscript[t, n]),(Subscript[y, n])) a (Subscript[k, 2])=f((Subscript[t, n])+(h/2),(Subscript[y, n])+((h (Subscript[k, 1]))/2)).

V systéme Mathematica naprogramujeme dve procedúry, eulermod a eulermodgraf. Tieto procedúry nám umožnia vykonať výpočty tabuliek hodnôt a nakresliť grafy približných riešení získaných modifikovanou Eulerovou metódou:

eulermod [f_,  h_,  ini_,a_,b_]:=
Module [ {yrk, t, y, rktable1,c},
  c=(b-a)/h;
  yrk[0]= ini;
  t[n_]:= a+n h; 
  yrk[n_]:= 
    Module[{k1, k2}, 
      k1=f[t[n-1], yrk[n-1]];
      k2=f[t[n-1]+(h/2), yrk[n-1]+(h/2)  k1];
      yrk[n]=yrk[n-1] +h k2];
  rktable1 = Table[yrk[i], {i,0,c}];
Table[{t[i],rktable1[[i+1]]},{i,0,c}]//TableForm]

eulermodgraf[f_,h_,ini_,a_,b_]:=Module[{yrk,t,y,rktable1,c},c=((b-a)/h);yrk[0]=ini;t[n_]:=a+n h;yrk[n_]:=Module[{k1,k2},k1=f[t[n-1],yrk[n-1]];k2=f[t[n-1]+(h/2),yrk[n-1]+((h k1)/2)];yrk[n]=yrk[n-1]+h k2];rktable1=Table[yrk[i],{i,0,c}];ListPlot[Table[{t[i],rktable1[[i+1]]},{i,0,c}],Joined->True,PlotStyle->{RGBColor[1,0,0]},PlotRange->All];Print["y[",t[c],"]=",rktable1[[c+1]]]]

kde f je funkcia zodpovedajúca diferenciálnej rovnici, h je dĺžka kroku, ini je hodnota začiatočnej podmienky a a a b sú hranice intervalu.

b) Heunova metóda, ktorá zodpovedá voľbe koeficientov (Subscript[b, 1])=(1/2),(Subscript[b, 2])=(1/2),(Subscript[c, 2])=1. Jej iteračná schéma má tvar (Subscript[t, n+1])=(Subscript[t, n])+h,  (Subscript[y, n+1])=(Subscript[y, n])+ ((h ((Subscript[k, 1])+(Subscript[k, 2])))/2) s (Subscript[k, 1])=f((Subscript[t, n]),(Subscript[y, n])) a (Subscript[k, 2])=f((Subscript[t, n])+h,(Subscript[y, n])+h (Subscript[k, 1])).

V systéme Mathematica naprogramujeme dve procedúry, mejoreuler a mejoreulergraf. Tieto vypočítajú tabuľku hodnôt a poskytnú grafické riešenie použitej Heunovej metódy:

mejoreuler [f_,  h_,  ini_,a_,b_]:=
Module [ {yrk, t, y, rktable1,c},
  c=(b-a)/h;
  yrk[0]= ini;
  t[n_]:= a+n h; 
  yrk[n_]:= 
    Module[{k1, k2}, 
      k1=f[t[n-1], yrk[n-1]];
      k2=f[t[n-1]+h, yrk[n-1]+h k1];
      yrk[n]=yrk[n-1] +(h/2) (k1+k2)];
  rktable1 = Table[yrk[i], {i,0,c}];
Table[{t[i],rktable1[[i+1]]},{i,0,c}]//TableForm]

mejoreulergraf[f_,h_,ini_,a_,b_]:=Module[{yrk,t,y,rktable1,c},c=((b-a)/h);yrk[0]=ini;t[n_]:=a+n h;yrk[n_]:=Module[{k1,k2},k1=f[t[n-1],yrk[n-1]];k2=f[t[n-1]+h,yrk[n-1]+h k1];yrk[n]=yrk[n-1]+(1/2) h (k1+k2)];rktable1=Table[yrk[i],{i,0,c}];ListPlot[Table[{t[i],rktable1[[i+1]]},{i,0,c}],Joined->True,PlotStyle->{RGBColor[1,0,0]},PlotRange->All];Print["y[",t[c],"]=",rktable1[[c+1]]]]

kde f je funkcia zodpovedajúca diferenciálnej rovnici, h je dĺžka kroku, ini je hodnota začiatočnej podmienky a a a b sú hranice intervalu.

Obidve metódy je možné použiť na nájdenie riešenia začiatočnej úlohy y'=-t y+((4t)/y), y(0)=1 na intervalu [0,1], s dĺžkou kroku 0.1.

f[t_,y_]=-t y+((4 t)/y);

eulermod[f,0.1,1,0,1]

eulermodgraf[f,0.1,1,0,1]

mejoreuler[f,0.1,1,0,1]

mejoreulergraf[f,0.1,1,0,1]

V predchádzajúcej kapitole 9 - opisujúcej jednokrokovú Eulerovu metódu bolo ukázané, že presné riešenie začiatočnej úlohy pre  t=1 je 1.70187. Preto v tomto prípade riešenia, ktoré získame pomocou modifikácie tejto metódy, sú lepšie ako aproximácia získaná pôvodným postupom. Vypovedá o tom aj rád metódy.

Klasická metóda Runge-Kutta tretieho rádu s tromi stupňami má nasledujúcu Butcherovu tabuľku

({{0, 0, ., .}, {(Subscript[c, 2]), (Subscript[c, 2]), 0, .}, {(Subscript[c, 3]), (Subscript[c, 3])-(Subscript[a, 32]), (Subscript[a, 32]), 0}, {., (Subscript[b, 1]), (Subscript[b, 2]), (Subscript[b, 3])}})

kde koeficienty iteračnej schémy musia spĺňať sústavu rovníc:

(Subscript[b, 1])+(Subscript[b, 2])+(Subscript[b, 3])=1
(Subscript[b, 2])(Subscript[c, 2])+(Subscript[b, 3])(Subscript[c, 3])= (1/2)
(Subscript[b, 2])((Subscript[c, 2])^2)+(Subscript[b, 3])((Subscript[c, 3])^2)= (1/3)
(Subscript[b, 3])(Subscript[c, 2])(Subscript[a, 32])= (1/6)

Rovnice sú závislé, preto existuje nekonečný počet metód Runge-Kutta tretieho rádu. Najpoužívanejšia je metóda:

({{0, 0, ., .}, {(1/2), (1/2), 0, .}, {1, -1, 2, 0}, {., (1/6), (2/3), (1/6)}})

Jej iteračná schéma má tvar (Subscript[t, n+1])=(Subscript[t, n])+h,  (Subscript[y, n+1])=(Subscript[y, n])+ ((h ((Subscript[k, 1])+4(Subscript[k, 2])+(Subscript[k, 3])))/6), kde (Subscript[k, 1])=f((Subscript[t, n]),(Subscript[y, n])), (Subscript[k, 2])=f((Subscript[t, n])+(h/2),(Subscript[y, n])+((h (Subscript[k, 1]))/2)) a (Subscript[k, 3])=f((Subscript[t, n])+h,(Subscript[y, n])+2h (Subscript[k, 2])-h (Subscript[k, 1])).

V systéme Mathematica naprogramujeme dve procedúry Runge3 a Runge3graf. Tieto vypočítajú tabuľku hodnôt a poskytnú grafickú reprezentáciu riešenia pomocou metódy Runge-Kutta tretieho rádu:

Runge3 [f_,  h_, ini_, a_, b_]  := 
Module [ {yrk3, t, rktable3,c },
  c=(b-a)/h;
  yrk3[0]= ini;
  t[n_]:=a+n h; 
  yrk3[n_]:= 
    Module[{k1, k2, k3}, 
      k1 = f[t[n-1], yrk3[n-1]];
      k2 = f[t[n-1] + h/2, yrk3[n-1] + (h/2) k1];
      k3 = f[t[n-1] + h, yrk3[n-1] -h k1 + 2 h k2];
    yrk3[n]= yrk3[n-1] + h ((1/6) k1 +(2/3) k2 +(1/6) k3)];
  rktable3 = Table[yrk3[i], {i,0,c}];
Table[{t[i],rktable3[[i+1]]},{i,0,c}]//TableForm]

Runge3graf[f_,h_,ini_,a_,b_]:=Module[{yrk3,t,rktable3,c},c=((b-a)/h);yrk3[0]=ini;t[n_]:=a+n h;yrk3[n_]:=Module[{k1,k2,k3},k1=f[t[n-1],yrk3[n-1]];k2=f[t[n-1]+(h/2),yrk3[n-1]+((h k1)/2)];k3=f[t[n-1]+h,yrk3[n-1]-h k1+2 h k2];yrk3[n]=yrk3[n-1]+h ((k1/6)+((2 k2)/3)+(k3/6))];rktable3=Table[yrk3[i],{i,0,c}];ListPlot[Table[{t[i],rktable3[[i+1]]},{i,0,c}],Joined->True,PlotStyle->{RGBColor[1,0,0]},PlotRange->All];Print["y[",t[c],"]=",rktable3[[c+1]]]]

kde f je funkcia zodpovedajúca diferenciálnej rovnici, h je dĺžka kroku, ini je hodnota začiatočnej podmienky a a a b sú hranice intervalu.

Obidve metódy sa použijú na vyriešenie začiatočnej úlohy y'=-t y+((4t)/y), y(0)=1  na intervale  [0,1] s dĺžkou kroku 0,1.

f[t_,y_]=-t y+((4 t)/y);

Runge3[f,0.1,1,0,1]

Runge3graf[f,0.1,1,0,1]

Najčastejšie používaná metóda je Runge-Kutta metóda štvrtého rádu, ktorá je reprezentovaná Butcherovou tabuľkou:

({{0,  0,  ,  ,  }, {(1/2), (1/2),  0,  ,  }, {(1/2), 0, (1/2),  0,  }, {1, 0, 0, 1,  0}, { , (1/6), (1/3), (1/3), (1/6)}})

a jej iteračná schéma má tvar

(Subscript[y, n+1])=(Subscript[y, n])+(h/(6 ))((Subscript[K, 1])+2(Subscript[K, 2])+2(Subscript[K, 3])+(Subscript[K, 4])) 
        (Subscript[K, 1])=f((Subscript[t, n]),(Subscript[y, n]))
        (Subscript[K, 2])=f((Subscript[t, n])+(h/(2 )),(Subscript[y, n])+(h/(2 ))(Subscript[K, 1]))
          (Subscript[K, 3])=f((Subscript[t, n])+(h/(2 )),(Subscript[y, n])+(h/(2 ))(Subscript[K, 2])) 
          (Subscript[K, 4])=f((Subscript[t, n])+h,(Subscript[y, n])+h (Subscript[K, 3]))
          (Subscript[t, n+1])=(Subscript[t, n])+h .

Znovu v systéme Mathematica naprogramujeme procedúry Runge4 a Runge4graf . Tieto vypočítajú tabuľku hodnôt a poskytnú grafickú reprezentáciu riešenia pomocou metódy Runge-Kutta štvrtého rádu:

Runge4 [f_,  h_, ini_, a_, b_]  := 
Module [ {yrk4, t, rktable4,c },
  c=(b-a)/h;
  yrk4[0]= ini;
  t[n_]:=a+n h; 
  yrk4[n_]:= 
    Module[{k1, k2, k3, k4 }, 
      k1 = f[t[n-1], yrk4[n-1]];
      k2 = f[t[n-1] + h/2, yrk4[n-1] + (h/2) k1];
      k3 = f[t[n-1] + h/2, yrk4[n-1] + (h/2) k2];
      k4 = f[t[n-1]+h, yrk4[n-1] + h k3];
      yrk4[n]= 
        yrk4[n-1] + 
          (h/6) (k1+2 k2+2 k3+k4)];
  rktable4 = Table[yrk4[i], {i,0,c}];
Table[{t[i],rktable4[[i+1]]},{i,0,c}]//TableForm]

Runge4graf[f_,h_,ini_,a_,b_]:=Module[{yrk4,t,rktable4,c},c=((b-a)/h);yrk4[0]=ini;t[n_]:=a+n h;yrk4[n_]:=Module[{k1,k2,k3,k4},k1=f[t[n-1],yrk4[n-1]];k2=f[t[n-1]+(h/2),yrk4[n-1]+((h k1)/2)];k3=f[t[n-1]+(h/2),yrk4[n-1]+((h k2)/2)];k4=f[t[n-1]+h,yrk4[n-1]+h k3];yrk4[n]=yrk4[n-1]+(1/6) h (k1+2 k2+2 k3+k4)];rktable4=Table[yrk4[i],{i,0,c}];ListPlot[Table[{t[i],rktable4[[i+1]]},{i,0,c}],Joined->True,PlotStyle->{RGBColor[1,0,0]},PlotRange->All];Print["y[",t[c],"]=",rktable4[[c+1]]]]

kde f je funkcia zodpovedajúca diferenciálnej rovnici, h je dĺžka kroku, ini je hodnota začiatočnej podmienky a a a b sú hranice intervalu.

Obidve metódy použijeme na nájdenie riešenia začiatočnej úlohy y'=-t y+((4t)/y), y(0)=1 na intervale  [0,1], s dĺžkou kroku 0,1.

f[t_,y_]=-t y+((4 t)/y);

Runge4[f,0.1,1,0,1]

Runge4graf[f,0.1,1,0,1]

Príklad 1.

Použime metódu Runge-Kutta  štvrtého rádu na nájdenie riešenia začiatočnej úlohy y'=(((y^2)-3(t^2)-2t y)/((t^2)+2t y)), y(1)=2, na intervale [1,2] s krokom h=0.1.

Riešenie

Definujme funkciu zodpovedajúcu pravej strane rovnice v zadanej začiatočnej úlohe a použime pripravené procedúry na nájdenie riešenia.

f[t_,y_]:=(((y^2)-3(t^2)-2t y)/((t^2)+2t y))

Runge4[f,0.1,2,1,2]

Runge4graf[f,0.1,2,1,2]

Príklad 2.

Riešte začiatočnú úlohu y'=(Sqrt[y])-((20(e^(-100(((t-2))^2))))/(Sqrt[π])), y(1)=1 . Nájdite približnú hodnotu riešenia v bode t=3, použitím procedúry Runge4 a s krokom dĺžky h=0.01. Urobte analýzu grafu a porovnajte ju so začiatočnou úlohou.

Riešenie

Pri pokuse o nájdenie presného riešenia začiatočnej úlohy vidíme, že systém Mathematica túto začiatočnú úlohu nedokáže vyriešiť:

DSolve[{y'[t]==(Sqrt[y[t]])-((20(E^(-100(((t-2))^2))))/(Sqrt[π])),y[1]==1},y[t],t]

Definujme funkciu zodpovedajúcu pravej strane začiatočnej úlohy a aplikujeme metódu Runge-Kutta 4. rádu

f[t_,y_]:=(Sqrt[y])-((20(E^(-100(((t-2))^2))))/(Sqrt[π]))

Runge4graf[f,0.01,1,1,3]

Teraz vyriešime pomocou rovnakej procedúry aj druhú začiatočnú úlohu.

f[t_,y_]:=(Sqrt[y])

Runge4graf[f,0.01,1,1,3]

Grafy sa približne zhodujú až do hodnoty t = 2. Potom riešenie druhej začiatočnej úlohy začne rýchlo divergovať. Náhly pokles v grafe riešenia prvej začiatočnej úlohy je spôsobený impulzom z pravej strany rovnice, ktorý zmení správanie riešení tejto rovnice. Impulz môžeme nakresliť nasledovne:

Plot[((20 (E^(-100 (((t-2))^2))))/(Sqrt[π])),{t,1,4},PlotRange->All]