Metoda Runge-Kutta

Aktuální verze stránky ještě nebyla zkontrolována zkušenými přispěvateli a může se výrazně lišit od verze recenzované 12. května 2021; kontroly vyžadují 2 úpravy .

Runge-Kuttovy metody ( v literatuře je název Runge-Kuttovy metody ) je rozsáhlá třída numerických metod pro řešení Cauchyho úlohy pro obyčejné diferenciální rovnice a jejich soustavy. První metody této třídy navrhli kolem roku 1900 němečtí matematici K. Runge a M. V. Kutta .

Třída metod Runge-Kutta zahrnuje explicitní Eulerovu metodu a modifikovanou Eulerovu metodu s přepočtem, což jsou metody prvního a druhého řádu přesnosti. Existují standardní explicitní metody třetího řádu přesnosti, které nejsou široce používány. Nejčastěji používanou a implementovanou v různých matematických balíčcích ( Maple , MathCAD , Maxima ) je klasická metoda Runge-Kutta , která má čtvrtý řád přesnosti. Při provádění výpočtů se zvýšenou přesností se stále častěji používají metody pátého a šestého řádu přesnosti [1] [2] . Konstrukce obvodů vyššího řádu je spojena s velkými výpočetními obtížemi [3] .

Metody sedmého řádu musí mít alespoň devět stupňů a metody osmého řádu musí mít alespoň 11 stupňů. U metod devátého a vyšších řádů (které však nemají velký praktický význam) není známo, kolik stupňů je potřeba k dosažení odpovídajícího řádu přesnosti [3] .

Klasická metoda Runge-Kutta čtvrtého řádu

Metoda čtvrtého řádu Runge-Kutta pro výpočty s konstantním integračním krokem je tak rozšířená, že se často nazývá jednoduše Runge-Kutta metoda.

Zvažte Cauchyho problém pro systém obyčejných diferenciálních rovnic prvního řádu. (Dále a ). ${\displaystyle \mathbf {y} ,\mathbf {f} ,\mathbf {k} _{i}\in \mathbb {R} ^{n))$ $x,h\in {\mathbb {R}}^{1}$

{\textbf {y}}'={\textbf {f}}(x,{\textbf {y}}),\quad {\textbf {y}}(x_{0})={\textbf {y} }_{0}.

Potom se přibližná hodnota v následujících bodech vypočítá pomocí iteračního vzorce:

\textbf{y}_{n+1} = \textbf{y}_n + {h \over 6}(\textbf{k}_1 + 2\textbf{k}_2 + 2\textbf{k}_3 + \ textbf{k}_4)

Výpočet nové hodnoty probíhá ve čtyřech fázích:

{\textbf {k}}_{1}={\textbf {f}}\left(x_{n},{\textbf {y}}_{n}\right),

{\textbf {k}}_{2}={\textbf {f}}\left(x_{n}+{h \over 2},{\textbf {y}}_{n}+{h \over 2}{\textbf {k}}_{1}\right),

{\textbf {k}}_{3}={\textbf {f}}\left(x_{n}+{h \over 2},{\textbf {y}}_{n}+{h \over 2}{\textbf {k}}_{2}\right),

{\textbf {k}}_{4}={\textbf {f}}\left(x_{n}+h,{\textbf {y}}_{n}+h\ {\textbf {k}} _{3}\vpravo).

kde je hodnota kroku mřížky v . $h$ $X$

Tato metoda má čtvrtý řád přesnosti. To znamená, že chyba v jednom kroku je řádná a celková chyba v posledním integračním intervalu je řádná . $O(h^{5})$ $O(h^4)$

Explicitní metody Runge-Kutta

Rodina explicitních metod Runge-Kutta je zobecněním jak explicitní Eulerovy metody, tak klasické metody Runge-Kutta čtvrtého řádu. Je to dáno vzorci

\textbf{y}_{n+1} = \textbf{y}_n + h\sum_{i=1}^s b_i \textbf{k}_i,

kde je hodnota kroku mřížky a výpočet nové hodnoty probíhá v následujících krocích: $h$ $X$ $s$

{\begin{array}{ll}{\textbf {k}}_{1}=&{\textbf {f}}(x_{n},{\textbf {y}}_{n}) ,\\{\textbf {k}}_{2}=&{\textbf {f}}(x_{n}+c_{2}h,{\textbf {y}}_{n}+a_{21 }h{\textbf {k}}_{1}),\\\cdots &\\{\textbf {k}}_{s}=&{\textbf {f}}(x_{n}+c_{ s}h,{\textbf {y}}_{n}+a_{s1}h{\textbf {k}}_{1}+a_{s2}h{\textbf {k}}_{2}+ \cdots +a_{s,s-1}h{\textbf {k}}_{s-1})\end{array}}

Konkrétní metoda je určena počtem a koeficienty a . Tyto koeficienty jsou často uspořádány do tabulky (nazývané Butcher table ): $s$ $b_{i},a_{{ij}}$ $c_{i}$

{\begin{array}{c|ccccc}0&&&&&\\c_{2}&a_{{21}}&&&&\\c_{3}&a_{{31}}&a_{{32}}&&&\\\vdots &\ vdots &\vdots &\ddots &&\\c_{s}&a_{{s1}}&a_{{s2}}&\tečky &a_{{ss-1}}&\\\hline &b_{1}&b_{2} &\tečky &b_{{s-1}}&b_{s}\end{array}}

Pro koeficienty metody Runge-Kutta musí být splněny podmínky pro . Pokud chcete, aby metoda měla řád , měli byste uvést i podmínku $\sum _{{j=1}}^{{i-1}}a_{{ij}}=c_{i}$ $i=2,\ldots ,s$ $p$

{\bar {{\textbf {y}}}}(h+x_{0})-{{\textbf {y}}}(h+x_{0})=O(h^{{p+1} }),

kde je aproximace získaná metodou Runge-Kutta. Po vícenásobné derivaci je tato podmínka transformována do soustavy polynomických rovnic s ohledem na koeficienty metody. ${\bar {{\textbf {y))))(h+x_{0})$

Implicitní metody Runge-Kutta

Všechny dosud uvedené metody Runge-Kutta jsou explicitní metody . Bohužel explicitní metody Runge-Kutta zpravidla nejsou vhodné pro řešení tuhých rovnic kvůli malé oblasti jejich absolutní stability [4] . Nestabilita explicitních metod Runge-Kutta také vytváří velmi vážné problémy při numerickém řešení parciálních diferenciálních rovnic .

Nestabilita explicitních metod Runge-Kutta motivovala vývoj implicitních metod. Implicitní metoda Runge-Kutta má tvar [5] [6] :

y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i},

kde

k_{i}=f{\bigl (}x_{n}+c_{i}h,y_{n}+h\sum _{j=1}^{s}a_{ij}k_{j }{\bigr )},\quad i=1,\ldots ,s.

Explicitní metoda se vyznačuje tím, že matice koeficientů pro ni má nižší trojúhelníkový tvar (včetně nulové hlavní diagonály) - na rozdíl od implicitní metody, kde má matice libovolný tvar. To je také vidět v Butcher's table . $a_{{ij}}$

{\begin{array}{c|cccc}c_{1}&a_{{11}}&a_{{12}}&\tečky &a_{{1s}}\\c_{2}&a_{{21}}&a_{ {22}}&\dots &a_{{2s}}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{s}&a_{{s1}}&a_{{s2}}&\tečky &a_{{ss}}\\\hline &b_{1}&b_{2}&\tečky &b_{s}\\\end{array}}={\begin{array}{c|c}{\mathbf {c }}&A\\\hline &{\mathbf {b^{T}}}\\\end{array}}

Důsledkem tohoto rozdílu je potřeba řešit v každém kroku soustavu rovnic pro , kde je počet stupňů. To zvyšuje výpočetní náklady, nicméně při dostatečně malé hodnotě lze aplikovat princip kontrakčních mapování a řešit tento systém jednoduchou iterací [7] . V případě jedné iterace to zvyšuje výpočetní náklady pouze dvakrát. $k_{i},i=1,2,...,s$ $s$ $h$

Na druhé straně Jean Kunzman (1961) a John Butcher (1964) ukázali, že pro libovolný počet fází existuje implicitní metoda Runge-Kutta s řádem přesnosti . To například znamená, že pro výše popsanou explicitní čtyřstupňovou metodu čtvrtého řádu existuje implicitní analog s dvojnásobnou přesností. $s$ $p=2s$

Implicitní metoda Runge-Kutta druhého řádu

Nejjednodušší implicitní metodou Runge-Kutta je upravená Eulerova metoda „s přepočtem“. Je dán vzorcem:

$\mathbf {y} _{n+1}=\mathbf {y} _{n}+h{\frac {\mathbf {f} (x_{n},\mathbf {y} _{n} )+\mathbf {f} (x_{n+1},\mathbf {y} _{n+1})}{2}}$ .

K jeho implementaci v každém kroku jsou vyžadovány alespoň dvě iterace (a dvě vyhodnocení funkcí).

Předpověď:

{\tilde {\mathbf {y} }}_{n+1}=\mathbf {y} _{n}+h\mathbf {f} (x_{n},\mathbf {y} _{ n})

Oprava:

\mathbf {y} _{n+1}=\mathbf {y} _{n}+h{\frac {\mathbf {f} (x_{n},\mathbf {y} _{n} )+\mathbf {f} (x_{n+1},{\tilde {\mathbf {y} }}_{n+1})}{2}}

Druhý vzorec je jednoduchou iterací řešení soustavy rovnic vzhledem k , zapsanou ve formě zobrazení kontrakce. Pro zlepšení přesnosti lze iterační korekci provést několikrát nahrazením . Modifikovaná Eulerova metoda „s přepočtem“ má druhý řád přesnosti. $\mathbf {y} _{n+1}$ ${\tilde {\mathbf {y} }}_{n+1}=\mathbf {y} _{n+1}$

Udržitelnost

Výhodou implicitních metod Runge-Kutta ve srovnání s explicitními je jejich větší stabilita, což je důležité zejména při řešení tuhých rovnic . Uvažujme jako příklad lineární rovnici y' = λ y . Obvyklá metoda Runge-Kutta použitá na tuto rovnici se redukuje na iteraci , přičemž r se rovná $y_{{n+1}}=r(h\lambda)\,y_{n}$

r(z)=1+zb^{T}(I-zA)^{{-1}}e={\frac {\det(I-zA+zeb^{T})}{\det(I- zA)}},

kde označuje sloupcový vektor jednotek [8] . Funkce se nazývá funkce stability [9] . Je to vidět ze vzorce, který je poměrem dvou polynomů stupně , pokud má metoda stupně. Explicitní metody mají přísně nižší trojúhelníkovou matici , z čehož vyplývá, že a že funkce stability je polynom [10] . $E$ $r$ $r$ $s$ $s$ $A,$ $\det(I-zA)=1,$

Numerické řešení tohoto příkladu konverguje k nule za podmínky c . Množina takových se nazývá oblast absolutní stability . Konkrétně se metoda nazývá A-stabilní , pokud jsou všechna c v oblasti absolutní stability. Stabilizační funkce explicitní metody Runge-Kutta je polynom, takže explicitní metody Runge-Kutta nemohou být v zásadě A-stabilní [10] . $|r(z)|<1$ $z=h\lambda$ $z$ $z$ ${\textrm {Re}}(z)<0$

Pokud má metoda pořadí , pak funkce stability splňuje podmínku pro . Zajímavý je tedy poměr polynomů daného stupně, který nejlépe aproximuje exponenciální funkci. Tyto vztahy jsou známé jako Padé aproximanty . Padého aproximant s čitatelem stupňů a jmenovatelem stupňů je A-stabilní tehdy a jen tehdy, když [11] $p$ $r(z)={\textrm {e}}^{z}+O(z^{{p+1}})$ $z\na 0$ $m$ $n$ $m\leq n\leq m+2.$

$s$ Metoda -stage Gauss-Legendre má řád , takže její stabilizační funkcí je Padého aproximace . To znamená, že metoda je A-stabilní [12] . To ukazuje, že A-stabilní metody Runge-Kutta mohou být libovolně vysokého řádu. Naproti tomu řád A-stability Adamsových metod nemůže překročit dva. $2s$ $m=n=s$

Výslovnost

Podle gramatických norem ruského jazyka se příjmení Kutta odmítá, takže říkají: "Metoda Runge-Kutta." Pravidla ruské gramatiky předepisují skloňovat všechna příjmení (včetně mužských) končící na -а, -я, před kterými je souhláska. Jedinou výjimkou jsou příjmení francouzského původu s přízvukem na poslední slabice jako Dumas, Zola [13] . Někdy však existuje nesklonná verze „Metody Runge-Kutta“ (například v knize [14] ).

Příklad řešení v algoritmických programovacích jazycích

y''+4y=\cos 3x,\quad y(0)=0,8,\y'(0)=2,\ x\in [0,1],\ h=0,1

provedením výměny a přenesením na pravou stranu získáme systém: $y'=z$ $4r$

{\begin{cases}y'=z=g(x,y,z)\\z'=\cos(3x)-4y=f(x,y,z)\end{cases))

Java kód pro řešení systému diferenciálních rovnic metodou Runge-Kutta public class MainClass { public static void main ( String [] args ) { int k = 2 ; dvojité Xo , Yo , Y1 , Zo , Z1 ; dvojité k1 , k2 , k4 , k3 , h ; dvojité ql , q2 , q4 , q3 ; /* *Počáteční podmínky */ Xo = 0 ; Yo = 0,8 ; Zo = 2 ; h = 0,1 ; // krok Systém . ven . println ( "\tX\t\tY\t\tZ" ); pro (; r ( Xo , 2 ) < 1,0 ; Xo += h ){ k1 = h * f ( Xo , Yo , Zo ); ql = h * g ( Xo , Yo , Zo ); k2 = h * f ( Xo + h / 2,0 , Yo + ql / 2,0 , Zo + kl / 2,0 ); q2 = h * g ( Xo + h / 2,0 , Yo + ql / 2,0 , Zo + kl / 2,0 ); k3 = h * f ( Xo + h / 2,0 , Yo + q2 / 2,0 , Zo + k2 / 2,0 ); q3 = h * g ( Xo + h / 2,0 , Yo + q2 / 2,0 , Zo + k2 / 2,0 ); k4 = h * f ( Xo + h , Yo + q3 , Zo + k3 ); q4 = h * g ( Xo + h , Yo + q3 , Zo + k3 ); Z1 = Zo + ( k1 + 2,0 * k2 + 2,0 * k3 + k4 ) / 6,0 ; Y1 = Yo + ( ql + 2,0 * q2 + 2,0 * q3 + q4 ) / 6,0 ; Systém . ven . println ( "\t" + r ( Xo + h , k ) + "\t\t" + r ( Y1 , k ) + "\t\t" + r ( Z1 , k )); Yo = Y1 ; Zo = Z1 ; } } /** * funkce pro zaokrouhlení a vyřazení "ocasu" */ public static double r ( double value , int k ){ return ( double ) Math . kolo (( Math . pow ( 10 , k ) * hodnota )) / Math . pow ( 10 , k ); } /** * funkce, které jsou získány ze systému */ public static double f ( double x , double y , double z ){ return ( Math . cos ( 3 * x ) - 4 * y ); } public static double g ( double x , double y , double z ){ return ( z ); } } Kód v C# pomocí System ; pomocí System.Collections.Generic ; jmenný prostor PRJ_RungeKutta { /// <summary> /// Implementace metody Runge-Kutta pro obyčejnou diferenciální rovnici /// </summary> public abstract class RungeKutta { /// <summary> /// Aktuální čas /// </summary > public double t ; /// <summary> /// Požadované řešení Y[0] je samotné řešení, Y[i] je i-tá derivace řešení /// </summary> public double [] Y ; /// <summary> /// Interní proměnné /// </summary> double [] YY , Y1 , Y2 , Y3 , Y4 ; chráněný dvojitý [] FY ; /// <summary> /// Konstruktor /// </summary> /// <param name="N">systémová dimenze</param> public RungeKutta ( uint N ) { Init ( N ); } /// <summary> /// Konstruktor /// </summary> public RungeKutta (){} /// <summary> /// Alokace paměti pro pracovní pole /// </summary> /// <param name="N">Rozměry polí</param> protected void Init ( uint N ) { Y = new double [ N ]; YY = nový double [ N ]; Y1 = nový dvojitý [ N ]; Y2 = nový double [ N ]; Y3 = nový dvojitý [ N ]; Y4 = nový dvojitý [ N ]; FY = new double [ N ]; } /// <summary> /// Nastavit počáteční podmínky /// </summary> /// <param name="t0">Čas zahájení</param> /// <param name="Y0">Počáteční podmínka </param> public void SetInit ( double t0 , double [] Y0 ) { t = t0 ; if ( Y == null ) Init (( uint ) Y0 . Délka ); for ( int i = 0 ; i < Y . Délka ; i ++) Y [ i ] = Y0 [ i ]; } /// <summary> /// Výpočet správných částí systému /// </summary> /// <param name="t">aktuální čas</param> /// <param name=" Y">vektorová řešení</param> /// <returns>pravá strana</returns> abstract public double [] F ( double t , double [] Y ); /// <summary> /// Další krok metody Runge-Kutta /// </summary> /// <param name="dt">aktuální časový krok (může být proměnný)</param> public void DalšíKrok ( double dt ) { int i ; if ( dt < 0 ) return ; // výpočet Y1 Y1 = F ( t , Y ); for ( i = 0 ; i < Y. Délka ; i ++) YY [ i ] = Y [ i ] + Y1 [ i ] * ( dt / 2,0 ) ; // výpočet Y2 Y2 = F ( t + dt / 2,0 , YY ); for ( i = 0 ; i < Y. Délka ; i ++) YY [ i ] = Y [ i ] + Y2 [ i ] * ( dt / 2,0 ) ; // výpočet Y3 Y3 = F ( t + dt / 2,0 , YY ); for ( i = 0 ; i < Y. Délka ; i ++ ) YY [ i ] = Y [ i ] + Y3 [ i ] * dt ; // výpočet Y4 Y4 = F ( t + dt , YY ); // výpočet řešení v novém kroku pro ( i = 0 ; i < Y . Délka ; i ++) Y [ i ] = Y [ i ] + dt / 6,0 * ( Y1 [ i ] + 2,0 * Y2 [ i ] + 2,0 * Y3 [ i ] + Y4 [ i ]); // výpočet aktuálního času t = t + dt ; } } class TMyRK : RungeKutta { public TMyRK ( uint N ) : základ ( N ) { } /// <summary> /// příklad matematického kyvadla /// y''(t)+y(t)=0 /// </summary> /// <param name="t">Čas</param > /// <param name="Y">Řešení</param> /// <returns>Pravá strana</returns> veřejné přepsání double [] F ( double t , double [] Y ) { FY [ 0 ] = Y [ 1 ]; FY [ 1 ] = - Y [ 0 ]; návrat FY ; } /// <summary> /// Příklad použití /// </summary> static public void Test () { // Časový krok double dt = 0,001 ; // TMyRK metoda objekt task = new TMyRK ( 2 ); // Definujte počáteční podmínky y(0)=0, y'(0)=1 úlohy double [] Y0 = { 0 , 1 }; // Nastavte počáteční podmínky pro úlohu task . SetInit ( 0 , Y0 ); // řešení do 15 sekund while ( task . t <= 15 ) { Console . WriteLine ( "Čas = {0:F5}; Func = {1:F8}; d Func / dx = {2:F8}" , úkol . t , úkol . Y [ 0 ], úkol . Y [ 1 ]); // výstup t, y, y' // výpočet dalšího kroku, integrační krok task . Další Krok ( dt ); } Konzole . ReadLine (); } } class Program { static void Main ( string [] args ) { TMyRK . test (); } } }

Program C# používá abstraktní třídu RungeKutta , která by měla přepsat abstraktní metodu F, která určuje pravé strany rovnic.

Příklad řešení v prostředí MATLAB

Řešení systémů diferenciálních rovnic metodou Runge-Kutta je jednou z nejběžnějších metod numerického řešení ve strojírenství. V prostředí MATLAB je implementována jedna z jeho variant - metoda Dorman-Prince . Chcete-li vyřešit soustavu rovnic, musíte nejprve napsat funkci, která vypočítá derivace, tedy funkce y = g ( x ,  y ,  z ) a z = cos(3 x ) − 4 y = f ( x ,  y ,  z ), jak je popsáno výše. V jednom z adresářů, které jsou přístupné ze systému MATLAB , musíte vytvořit textový soubor s názvem (například) runge.m s následujícím obsahem (pro MATLAB verze 5.3):

MATLAB , runge.m funkce Dy = rozsah ( x, y ) Dy = y (:); Dy ( 1 ) = y ( 2 ); Dy ( 2 ) = cos ( 3 * x ) -4 * y ( 1 ) ;

Název souboru a název funkce se musí shodovat, ale může to být cokoliv, co nebylo dříve použito.

Poté je potřeba vytvořit hlavní soubor s názvem např. main.m , který provede základní výpočty. Tento hlavní soubor bude obsahovat následující text:

MATLAB , hlavní.m jasný ; clc ; % Vymažte paměť a obrazovku h = 0,1 ; % Integrační krok x_fin = 8 ; % Finální integrační čas yo = 0,8 ; % Počáteční hodnota funkce Dy0 = 2 ; % Počáteční hodnota derivace funkce [ x , y ] = ode45 ( 'rozsah' , [ 0 : h : x_fin ], [ y0 Dy0 ]); % metoda Runge-Kutta plot ( x , y , 'LineWidth' , 2 ); mřížka ; % Plot a mřížka legenda ( 'y(x)' , 'y''(x)' , 0 ); % Legenda na grafu

Protože je MATLAB maticově orientovaný, lze řešení Runge-Kutta velmi snadno provést pro celý rozsah x , jako například ve výše uvedeném příkladu programu. Zde je řešením graf funkce v časech od 0 do x_fin .

Proměnné x a y vrácené funkcí ODE45 jsou hodnotové vektory. Je zřejmé, že řešením výše uvedeného příkladu je druhý prvek x , protože první hodnota je 0, integrační krok je h = 0,1 a sledovaná hodnota x = 0,1. Následující záznam v příkazovém okně MATLABu poskytne požadované řešení:

řešení MATLAB y1 = y ( najít ( x == 0,1 ))

Odpověď: y1 = 0,98768

Poznámky

↑ Bakhvalov N. S. , Zhidkov N. P. , Kobelkov G. M. . Numerické metody. - M . : Laboratoř základních znalostí, 2001. - 630 s. — ISBN 5-93208-043-4 . - S. 363-375.
↑ Ilyina V. A., Silaev P. K. . Numerické metody pro teoretické fyziky. Část 2. - Moskva-Iževsk: Institut počítačového výzkumu, 2004. - 118 s. — ISBN 5-93972-320-9 . - S. 16-30.
↑ 12 Butcher J. C. . Numerické metody pro obyčejné diferenciální rovnice. — New York: John Wiley & Sons , 2008. — ISBN 978-0-470-72335-7 .
↑ Süli & Mayers, 2003 , str. 349-351.
↑ Iserles, 1996 , str. 41.
↑ Süli & Mayers, 2003 , str. 351-352.
↑ Implicitní metody Runge – Kutta archivován 6. března 2019 na Wayback Machine .
↑ Hairer & Wanner, 1996 , s. 40-41.
↑ Hairer & Wanner, 1996 , s. 40.
↑ 1 2 Iserles, 1996 , str. 60.
↑ Iserles, 1996 , str. 62-63.
↑ Iserles, 1996 , str. 63.
↑ Jak odmítnout příjmení (obtížné případy) - "Gramota.ru" - referenční a informační internetový portál "Ruský jazyk" . Získáno 5. července 2016. Archivováno z originálu 23. května 2018. (neurčitý)
↑ Demidovich B. P., Maron I. A., Shuvalova E. Z. . Numerické metody analýzy. 3. vyd. — M .: Nauka , 1967.

Literatura

Hairer E., Wanner G. . Řešení obyčejných diferenciálních rovnic II: Tuhé a diferenciálně-algebraické úlohy. 2. vyd. - Berlin, New York: Springer-Verlag , 1996. - ISBN 978-3-540-60452-5 .
Iserles A. První kurz numerické analýzy diferenciálních rovnic. - Cambridge: Cambridge University Press , 1996. - ISBN 978-0-521-55655-2 .
Suli E., Mayers D. Úvod do numerické analýzy. - Cambridge: Cambridge University Press , 2003. - ISBN 0-521-00794-1 .

Metoda konečných rozdílů
Obecné články	Metoda konečných rozdílů rozdílové schéma diferenční rovnice
Typy rozdílových schémat	Eulerova metoda Metoda Runge-Kutta Adamsova metoda Integrace Werlet Metoda konečného rozdílu v časové doméně