Menu







Równania różniczkowe - Metody Rungego-Kutty

PL
Data dodania: 2011-09-18, Autor: Karol, Wyświetleń: 1363

Witam wszystkich bardzo serdecznie.

Dzisiaj chciałbym omówić metody Rungego-Kutty. Metody te tyczą się rozwiązywania równan różniczkowych zwyczajych. A mowiąc precyzyjniej - problemu Cauchy'ego, czyli zagadnienia z określonymi warunkami początkowymi.

Problem Cauchy'ego w postaci różniczkowej wygląda mniej więcej tak:

$$ \left\{ \begin{array}{l} x'(t)=f(t,x(t)) \\ x(t_0)=x_0 \end{array}\right. $$

Dla tych co widzą cos takiego pierwszy raz, to powiem, że jest to układ równań różniczkowych ze zdefiniowanymi warunkami początkowymi.

Kolejnym pytaniem jakie może się pojawić to gdzie jest ten układ równań rózniczkowy, skoro ja widzę tylko jedno:

$$ x'(t)=f(t,x(t)) $$

Otóż, dlatego nazywamy to układem, gdyż funkcja f może być funkcją macierzową, np.

$$ \left\{ \begin{array}{l} x'(t)=Ax(t) \\ x(0)=x_0 \end{array} \right. $$

gdzie A jest macierzą $(d \times d)$ Jest to klasyczne równanie różniczkowe liniowe, jednorodne o stałych współczynnikach :)

Wówczas $x(t)=(x_1(t), \cdots , x_d(t) )$

Ponadto, gdy równanie $x'(x) = f(t,x) = Ax$ spełnia warunek Lipschitza, tzn.

$$ |Ax-Ay|=|A(x-y)| \leq \|A\| \cdot |x-y| $$

gdzie $\|A\|$ jest normą macierzy A, a zarazem w tym wypadku stałą Lipschitza; Wówczas problem Cauchy'ego na mocy Twierdzenia Picaarda ma jednoznaczne rozwiązanie

$$ x\colon\lbrack t_0,T\rbrack \to \mathbb{R}^n $$

Rozwijąc funkcję w ciąg przybliżonych rozwiązan (który to też nadaje sie na metode numeryczna), dochodzimy do rozwiązania

$$ x(t)=exp(tA) \cdot x_0 $$

OK, tyle teorii. Teraz zajmijmy się zanalezieniem tego rozwiązania numerycznie. Wiemy, ze problem Cauchy'ego można przedstawić w postaci całkowej:

$$ x(t)=x_0 + \int_{t_0}^T f(s,x(s))ds \qquad t\in \lbrack t_0,T \rbrack $$

Wiec, aby znaleźć x(t) należałoby rozwiązać podaną całkę.

O najprostszym sposobie rozwiązania takiej całki już pisałem. Dzisiaj jednak chciałbym przedstawić najmocniejszą metodę stałokrokową jaką znam.

Wiemy, że aby policzyć całke oznaczaną numerycznie należy jak najlepiej przybliżyć przyrost funkcji po kroku h. Innymi słowy, jeżeli funkcja y(t) jest funkcją pierwotną f(t) to szukamy:

$$ y(t+h)-y(t)=\delta y $$

Ja podam jedynie wzory Ryngego-Kytty na obliczenie owego przyrostu. Szczegóły rachunku można znaleźć w książce Mineura, "Techniques du calcul numerique", str. 367

Zacznijmy od wzorów stopnia trzeciego, które przybliżająa nam przyrost z dokładnaścią $h^4$

$$ \begin{array}{rcl} k_1 & = & h f(x_0, y_0) \\ k_2 & = & h f(x_0+ah, y_0+\alpha k_1)\\ k_3 & = & h f(x_0+bh, y_0+\beta k_1 + \gamma k_2)\\ \\ y_1-y_0=\Delta y & = & R_1 k_1 + R_2 k_2 + R_3 k_3 \end{array} $$

Teraz należy dobrać wartości $a, b; \alpha, \beta, \gamma; R_1, R_2, R_3$ tak, by różnica pomiędzy dokładną wartością y i jego wartością przybliżoną była nieskończenie małą rzędu $h^4$. a oraz b możemy wybrać dowolnie otrzymując pozostałe stałe ze wzorów:

$$ \begin{array}{rcl} \alpha & = & a \\ \beta & = & b \frac{3a^2-3a+b}{a(3a-2)}\\ \gamma & = & \frac{b(a-b)}{a(3a-2)}\\ R_1 & = & \frac{6ab-3(a+b)+2}{6ab}\\ R_2 & = & \frac{2-3b}{6a(a-b)}\\ R_3 & = & \frac{3a-2}{6b(a-b)} \end{array} $$

Otrzymujemy zatem całą rodzinę wyrazów zależnych od dwóch parametrów a oraz b. Przyjmijmy przykładowo $a=\frac{1}{2}$, b=1, co prowadzi do formuł Rungego-Kutty stopnia 3

$$ \begin{array}{rcl} k_1 & = & h f(x_0, y_0) \\ k_2 & = & h f(x_0+h/2, y_0+ k_1/2)\\ k_3 & = & h f(x_0+h, y_0-k_1+ k_2/2)\\ \\ \Delta y & = & \frac{1}{6}(k_1 + 4k_2 + k_3) \end{array}{rcl} $$

Aby nie było niejasności objaśnię, o co chodzi w tych wzorkach. Otóż najpierw obliczamy wartości $k_i$ które są potrzebne do obliczenia kolejnych k i ostatecznie przyrostu funkcji. x jest zmienną całkowania (czyli np. czasem), a y funkcją. Dla funkcji zależnych tylko od czasu (całkowanych po czasie) istotny bedzie tylko pierwszy parametr. Dla funkcji zależnych tylko od położenia (całkowanych po czasie) istotny będzie tylko drugi parametr (przyklłdem takiej funkcji może być na przykład funkcja opisująca ruch cząstki w polu potencjału, bo pole potecjału zależy od położenia, a ruch od czasu - o tym bedzie więcej w przykladzie, za jakieś 2 strony). Dla funkcji zależnych i od położenia i od czasu (np. $x'(t)=x*\sin t$) istotne sa obydwa parametry.

A teraz bardzo prosty przykład wykorzystujący tą metodę (schemat)

function f(t:real):real; { funkcja zalezna tylko od czasu }

(...)

begin
 { zmienna czasowa (po tej calkujemu) ustawiamy na czas poczatkowy }
 t:=t_0; 

 wynik:=f_0; {ustawiamy warunki poczatkowe}

 while t<T_e do
  begin

   {do wyniku dodajemy przyrost obliczony ze wzorów} 
   wynik:=wynik+delta(t);  

   {w tym miejscu mozemy wypisywac wartosci funkcji otrzymujac
    dyskretny obrazy ewolucji ukladu}

  end;
 
 {w tym miejscu w zmiennej wynik mamy wartosc calki oznaczonej z
  funkcji f na przedziale [t_0,T]}

end.

Mam nadzieję, że to troszkę naświetliło do czego te wzorki służą. Ok, ale my idziemy dalej. Czas poprawić metodę, gdyż błąd $h^4$ zdecydowanie nam nie odpowiada.

Podam teraz formuły czwartego rzędu, które są bardzo wygodne w programowaniu ze względu na swoją prostotę. To właśnie te poniższe wzory noszą nazwę klasycznej metody Rungego-Kutty.

$$ \begin{array}{rcl} k_1 & = & h f(x_0, y_0) \\ k_2 & = & h f(x_0+h/2, y_0+ k_1/2)\\ k_3 & = & h f(x_0+h/2, y_0+ k_2/2)\\ k_4 & = & h f(x_0+h, y_0+k_3)\\ \\ \Delta y & = & \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{array} $$

Dobra, teraz będzie ciekawie. Milne w swojej książce "Numerical solutions of differential equations", na stronie 73 wyprowadził formuły stopnia szóstego. Błąd w jednym kroku jest rzędu $h^7$, czyli już nie najgorzej :). Oto te wzorki:

$$ \begin{array}{rcl} k_1 & = & h f(x_0, y_0) \\ k_2 & = & h f(x_0+\frac{h}{3}, y_0+ \frac{k_1}{3})\\ k_3 & = & h f(x_0+\frac{2h}{5}, y_0+ \frac{6k_2+4k_1}{25})\\ k_4 & = & h f(x_0+h, y_0 + \frac{15k_3-12k_2+k_1}{4})\\ k_5 & = & h f(x_0+\frac{2h}{3}, y_0+\frac{8k_4-50k_3+90k_2+6k_1}{81})\\ k_6 & = & h f(x_0+\frac{4h}{5}, y_0+\frac{8k_4+10k_3+36k_2+6k_1}{75})\\ \\ \Delta y & = & \frac{1}{192}(23k_1 + 125k_3 -81k_5+125k_6) \end{array} $$

Prawda, że ładne :) Ale przygotowałem jeszcze jeden zestaw wzorków :) Tym razem autorem ich wyprowadzenia jest nie jaki Huta. Jest to metoda stopnia (uwaga, uwaga) ósmego. Czyli błąd jaki popełniamy przy jednym kroku to $h^9$. Aby łatwiej można było sobie to wyobrazić, całkując funkcje krokiem 0.1 błąd jaki popełniamy to 0.000000001, a to jest już bardzo ładnie :). Jest to najsilniejsz metoda stałokrokowa jaką znam.

Dobra, pora na wzorki :]

$$ \begin{array}{rcl} k_1 & = & h f(x_0, y_0) \\ k_2 & = & h f(x_0+\frac{h}{9}, y_0+ \frac{k_1}{9})\\ k_3 & = & h f(x_0+\frac{h}{6}, y_0+ \frac{3k_2+k_1}{24})\\ k_4 & = & h f(x_0+\frac{h}{3}, y_0+ \frac{4k_3-3k_2+k_1}{6})\\ k_5 & = & h f(x_0+\frac{h}{2}, y_0+ \frac{6k_4-24k_3+27k_2-5k_1}{8})\\ \end{array} $$
$$ \begin{array}{rcl} k_6 & = & h f(x_0+\frac{2h}{3}, y_0+ \frac{k_5-102k_4+867k_3-981k_2+221k_1}{9})\\ k_7 & = & h f(x_0+\frac{5h}{6}, y_0+ \frac{3k_6+80k_5-66k_4-472k_3+678k_2-183k_1}{48})\\ k_8 & = & h f(x_0+h, y_0+ \frac{72k_7-9k_6-454k_5+834k_4+1002k_3-2072k_2+716k_1}{82})\\ \\ \Delta y & = & \frac{41(k_1+k_8)+216(k_3+k_7)+27(k_4+k_6)+272k_5}{840} \end{array} $$

OK. Przejdźmy teraz do przykładu jaki przygotowałem. Zagadnienie jakie idealnie nadaje się do rozwiązania tą metodą jest problem ruchu cząstki w polu potencjału.

Jednak, aby udało się rozwiązać to zadanie muszę powiedzieć jeszcze jak policzyć numerycznie pochodną.

Najłatwiejszym spodsobem na policzenie pochodnej w punkcie jest skorzystanie z definicji pochodnej, która wygląda tak:

$$ \lim_{\Delta x \to \infty}\frac{f(x_0+\Delta x)-f(x_0)}{\Delta x} $$

To jest definicja pochodnej funkcji f w punkcie $x_0$

Dla funkcji wielu zmiennych wygląda to tak

$$ \lim_{\Delta x \to \infty}\left(\frac{f(x_1,\cdots,x_i+\Delta x_i, \cdots, x_n)}{\delta x_i}-\frac{f(x_i, \cdots, x_n)}{\Delta x_i}\right) $$

i nazywa się pochodną cząstkową funkcji n zmiennych lub granicą ilorazu różnicowego funkcji f w punkcie $(x_1, \cdots, x_n)$.

Oznacza się ją przez $\frac{\partial f}{\partial x_i}$

Programuje się to dokładnie tak jak wygląda, czyli:

function pochodnaFwX1(x1,x2,x3:real):real;
 var delta:real;
begin
 delta:=x1/1e10; {wybieramy sobie jakas mala delte}

 {proste, łatwe i przyjemne :) }
 pochodnaFwX1:=F(x1+delta,x2,x3)/delta; 
end;

No to jak wiemy już jak policzyć pochodną, to do dzieła!

ZADANIE:

Przecałkować numerycznie na okres 0.7 dnia z krokiem 0.00001 dnia układ równań różniczkowych

$$ \ddot{x}=\frac{\partial U}{\partial x} \qquad , \qquad \ddot{y}=\frac{\partial U}{\partial y} \qquad , \qquad \ddot{z}=\frac{\partial U}{\partial z} $$

opisujących ruch cząstki w potencjału:

$$ U=\frac{\mu}{r} \left[ 1- I_4 \frac{R^4}{7r^4} \left( 3- 30 \frac{z^2}{r^4} + 35 \frac{z^4}{r^4} \right) \right] $$

dla wartości początkowych:

$$ \begin{array}{lll} x(0)=-2897.37km & y(0)=5911.49km & z(0)=1063.50km\\ \dot{x}(0)=41625.0\frac{km}{doba} & \dot{y}(0)=-266841.2\frac{km}{doba} & \dot{z}(0)=151165.7\frac{km}{doba} \end{array} $$

Parametry użyte w równaniach mają wartości:

$$ R=6378.14km \quad , \quad I_4=0.000001593\quad,\quad \mu=398600.47\frac{km^3}{s^2} $$
$$ r=\sqrt{x^2+y^2+z^2} $$

Tak BTW to było moje zadanie na zaliczenie metod numerycznych :)

Oto program:

{$A+,B-,D+,E+,F-,G+,I+,L+,N+,O-,P-,Q-,R-,S+,T-,V+,X+,Y+}
{$M 16384,0,655360}
uses crt;

type float=extended;    {typ zmiennych do obliczen}

var rysuj:boolean;      {rysowac czy nie}
const ile=6;            {ilosc rownan do wycalkowania}
type wektor=array[1..ile]of float; {typ do zmiennych funkcji}

const
      {STALE UZYWANE WE WZORZE}
      C_r  :float = 6378.14; {km}
      C_r4 :float = 6378.14*6378.14*6378.14*6378.14; {km^4}
      C_ni :float = 398600.47{km^3/s^2}*3600*24*3600*24; {km^3/doba^2}
      C_I4 :float = -0.000001593; {bezwymiarowe}

      {stala do policzenia epsilona (uzywane przy pochodnych)}
      Depsd=1e10;  

      swiatelko=15; {bajerek}

{funkcja opisujaca potencjal, zalezna od punktu przestrzeni}
function U(x,y,z:float):float; 
var r,r2,r4,z2,z4:float;
begin
 r:=sqrt(x*x+y*y+z*z);
 r2:=r*r;
 r4:=r2*r2;
 z2:=z*z;
 z4:=z2*z2;
 U:=(C_ni/r) * (1-  C_I4*(C_r4/(7*r4))*(3-30*(z2/r2)+35*(z4/r4))  );
end;

function Udx(x,y,z:float):float; {pochodna funkcji U po x}
var epsd:float;
begin
 epsd:=x/Depsd; {wyznaczamy sobie epsilona}
 Udx:=(U(x+epsd,y,z)-U(x,y,z))/epsd; {i liczymy pochodna z definicji}
end;

function Udy(x,y,z:float):float; {pochodna funkcji U po y}
var epsd:float;
begin
 epsd:=y/Depsd;
 Udy:=(U(x,y+epsd,z)-U(x,y,z))/epsd;
end;

function Udz(x,y,z:float):float;  {pochodna funkcji U po z}
var epsd:float;
begin
 epsd:=z/Depsd;
 Udz:=(U(x,y,z+epsd)-U(x,y,z))/epsd;
end;

procedure f(t:float;wej:wektor;var wyj:wektor);
{tutaj definiujemy funkcje, jakie beda calkowane}
{czyli po prosty prawe strony rownan rozniczkowych}
begin
 wyj[1]:=wej[4];
 wyj[2]:=wej[5];
 wyj[3]:=wej[6];

 wyj[4]:=Udx(wej[1],wej[2],wej[3]);
 wyj[5]:=Udy(wej[1],wej[2],wej[3]);
 wyj[6]:=Udz(wej[1],wej[2],wej[3]);
end;

{inicjuje tryb 13h potrzebny do narysowania wykresu}
procedure init13h;assembler; 
asm
  mov ax,0013h
  int 10h
end;

procedure close13h;assembler; {przywraca tryb tekstowy}
asm
  mov ax,0003h
  int 10h
end;

{rysuje pixel na pulpicie}
procedure putpixel(x,y:word;color:byte);assembler; 
{jak ktos nie wie o co ty chodzi to niech przeczyta artykul 
Tryb 13h w dziale pascal}
asm
 mov ax,0a000h
 mov es,ax
 mov di,x
 mov dx,y
 shl dx,6
 add di,dx
 shl dx,2
 add di,dx
 mov cl,color
 mov es:[di],cl
end;

procedure delta(t,h:float;wej:wektor; var wyj:wektor);
{liczy przyrost funkcji zaleznej od t i zmiennych 
w wektorze wej na przedzile czasu h}
var d,w:wektor;
    k1,k2,k3,k4,k5,k6,k7,k8:wektor;
    i:byte;

begin

{Metoda Rungego-Kutty stopnia 8}

{szczegoly -> patrz tresc artykulu}

 f(t,wej,d);
 for i:=1 to ile do k1_:=h*d_;

 for i:=1 to ile do w_:=wej_+k1_/9;
 f(t+h/9,w,d);
 for i:=1 to ile do k2_:=h*d_;

 for i:=1 to ile do w_:=wej_+(k1_+3*k2_)/24;
 f(t+h/6,w,d);
 for i:=1 to ile do k3_:=h*d_;

 for i:=1 to ile do w_:=wej_+(4*k3_-3*k2_+k1_)/6;
 f(t+h/3,w,d);
 for i:=1 to ile do k4_:=h*d_;

 for i:=1 to ile do w_:=
    wej_+(6*k4_-24*k3_+27*k2_-5*k1_)/8;
 f(t+h/2,w,d);
 for i:=1 to ile do k5_:=h*d_;

 for i:=1 to ile do 
    w_:=wej_+(k5_-102*k4_+867*k3_-981*k2_+221*k1_)/9;
 f(t+2*h/3,w,d);
 for i:=1 to ile do k6_:=h*d_;

 for i:=1 to ile do 
    w_:=wej_+(3*k6_+80*k5_-66*k4_-472*
          k3_+678*k2_-183*k1_)/48;
 f(t+5*h/6,w,d);
 for i:=1 to ile do k7_:=h*d_;

 for i:=1 to ile do 
    w_:=wej_+(72*k7_-9*k6_-454*k5_+834*k4_+1002*
          k3_-2072*k2_+716*k1_)/82;
 f(t+h,w,d);
 for i:=1 to ile do k8_:=h*d_;

 for i:=1 to ile do 
    wyj_:=(41*(k1_+k8_)+216*(k3_+k7_)+27*
            (k4_+k6_)+272*k5_)/840;
end;

var t,h:float;    {czas i krok}
    wej,wyj:wektor;  {wejscie i wyjscie}
    i:byte;         {pomocnicza }
    zoom:float; {skala wykresu (1:zoom) }
    x,y:word; {pomocniczew do rysowania}

begin
 t:=0; {czas poczatkowy 0}

 {krok musi byc taki maly, bo funkcja bardzo szybko sie zmienia}
 h:=0.00001;  
{x}       wej[1]:=-2897.37; {km}
{y}       wej[2]:=5911.49;  {km}
{z}       wej[3]:=1063.50;  {km}
{x. = Vx} wej[4]:=41625.0;  {km/doba}
{y. = Vy} wej[5]:=-266841.2;{km/doba}
{z. = Vz} wej[6]:=151165.7; {km/doba}

 zoom:=100; {skala wykresu 1:100, tj. 1punkt=100km }

 rysuj:=true; {rysujemy}

 if rysuj then  {jezeli rysujemy to : }
   begin
     init13h;   {inicjujemy tryb graficzny}
     for x:=1 to swiatelko do    {i rysujemy w gornym lewym rogu}
       for y:=1 to swiatelko do  {czerwony kwadracik, ktory mowi, }
         putpixel(x,y,4);        {ze liczenie jeszcze sie nie skonczylo}
   end;

 putpixel(160,100,4); {rysujemy czerwana kropke w srodku naszego ukladu}

 {liczymy ruch czastki na przedziale 0.7dnia, tj. niecale 17 godzin}
 while t<0.7 do 
 begin
  delta(t,h,wej,wyj); {liczymy nasze przyrosty zmiennych}

  {dodajemy przyrosty do zmiennych}
  for i:=1 to ile do wej_:=wej_+wyj_; 

  t:=t+h; {zwiekszamy czas o krok}

  if rysuj then  {rysyjemy punkcik}
     putpixel(160+round((wej[1]-wej[3])/zoom),100-round((wej[2]+wej[3])/zoom),15)
  else  {albo wypisujemy co wyliczylismy}
  writeln(t:0:2,' x=',wej[1]:0:2,' y=',wej[2]:0:2,' z=',wej[3]:0:2,
        ' Vx=',wej[4]:0:2,' Vy=',wej[5]:0:2,' Vz=',wej[6]:0:2);

  {jak ktos wcisnie [ESC] konczymy dzialanie programu}
  if keypressed then if readkey=#27 then 
  begin 
    if rysuj then close13h; 
    halt(1); 
    exit; 
  end;

 end;

 {jak rysujemy to czerwone swiatelko zmienia sie w zielone}
 if rysuj then  
  for x:=1 to swiatelko do
   for y:=1 to swiatelko do
    putpixel(x,y,2);

  repeat until keypressed;readkey; {czekamy az ktos wcisnie klawisz}

if rysuj then close13h; {i jak rysowalismy to wracamy do trybu tekstowego}
end.

To by było na tyle. Pozdrawiam

PS. LaTeX Rulez!!!

Oceń ten artykuł:

 

Aby dodawać komentarze musisz być zalogowany!


Polecane
Czytelnia
Kontakt

Jeśli chcesz się z nami skontaktować napisz na adres: info(at)binboy.org lub odwiedź nasz profil na Facebooku!

O Nas

Serwis binboy.org to kopalnia wiedzy dla wszystkich z branży IT, w szczególności dla programistów i webmasterów. To duży zbiór kursów programowania, tutoriali, darmowych ebooków, setki kodów źródłowych itp.

Bądź w kontakcie

Panel użytkownika

Zaloguj się do panelu użytkownika.
Nie masz konta? Zarejestruj się!
Zapomniałeś hasła?