Методи на Рунге-Кута

В математическия анализ, методите на Рунге-Кута са важна част от имплицитните и експлицитните итеративни методи за получаване на приближено решение на обикновени диференциални уравнения. Тези похвати са били въведени около 1900 г. от немските математици Карл Рунге и Мартин Вилхелм Кута.

Метод на Рунге-Кута от 4-ти ред

Методът на Рунге-Кута от 4-ти ред е най-често използваният, също известен като „РК“, „Класически метод на Рунге-Кута“ или просто „Метод на Рунге-Кута“.

Дадена е следната начална задача (още се нарича „задача на Коши“):

$\dot{y}=f(t,y),y(t_0)=y_0.$

Казано с думи, това означава, че стойностите, които получава y са функции на y и на t(time). В началото времето е $t_0$, a y е $y_0$. В уравнението y може да е скаларна или векторна променлива.

Методът на РК от 4-ти ред за тази задача е получен от следните формули:

$\begin{array}{cc}{y}_{n+1}& =y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\\ t_{n+1}& =t_n+h\\ \end{array}$

където $y_{n+1}$ приближението с РК4 на $y(t_{n+1})$, и

$\begin{array}{cc}{k}_1& =hf(t_n,y_n),\\ k_2& =hf(t_n+\frac{1}{2}h,y_n+\frac{1}{2}{k}_1),\\ k_3& =hf(t_n+\frac{1}{2}h,y_n+\frac{1}{2}{k}_2),\\ k_4& =hf(t_n+h,y_n+k_3).\end{array}$

(Забележка: горните уравнения имат различни, но равнозначни дефиниции в различните текстове).

Така новата стойност ($y_{n+1}$) се пресмята с помощта на текущата стойност($y_n$) плюс стойност на четирите нараствания с определени тегла, където всяко нарастване е произведение от големината на интервала h и оценения наклон, зададен от функцията f от дясната страна на диференциалното уравнение.

• $k_1$ е нарастване, основаващо се на наклона от началото на интервала, използвайки $y_n$, Метод на (Ойлер);
• $k_2$ е нарастване, основаващо се на наклона от междинна точка от интервала, използвайки $y_n+\frac{1}{2}{k}_1$;
• $k_3$ е отново нарастване, основаващо се на налкона от междинна точка, но използвайки $y_n+\frac{1}{2}{k}_2$;
• $k_4$ е нарастване, основаващо се на наклона в края на интервала, използвайки $y_n+k_3$.

При усредняване на четирите нараствания се дава по-голямо тегло на нарастванията в междинните точки. Теглата са избрани така, че ако $f$ не зависи от $y$ и диференциалното уравнение има пръв интеграл, тогава РК4 e „правило на Симпсън“.

РК4 е метод от 4-ти ред, което означава, че грешката на всяка стъпка е от порядъка на $h^5$, докато крайната натрупана грешка има ред $h^4$.

Явни методи на Рунге-Кута

Семейството на експлицитните методи на Рунге-Кута са обобщение на РК4 метода, споменат преди малко. Те са представени чрез:

$y_{n+1}=y_n+\sum _{i=1}^s{b}_i{k}_i,$

където

$k_1=hf(t_n,y_n),$

$k_2=hf(t_n+c_{2}h,y_n+a_{21}{k}_1),$

$k_3=hf(t_n+c_{3}h,y_n+a_{31}{k}_1+a_{32}{k}_2),$

$⋮$

$k_s=hf(t_n+c_{s}h,y_n+a_{s1}{k}_1+a_{s2}{k}_2+\cdots +a_{s,s-1}{k}_{s-1}).$

(Забележка: горните уравнения имат различни, но равнозначни дефиниции в различните текстове).

За да се специфицира определен метод, е нужно да се зададе цяло число s (брой етапи) и коефициентите a_ij (за 1 ≤ j < i ≤ s), b_i (за i = 1, 2, ..., s) и c_i (за i = 2, 3, ..., s). Матрицата [a_ij] се нарича матрица на Рунге-Кута, където b_i и c_i са познати като тегла и възли. Таблицата на Бутчер е следната:

	0
	$c_2$	$a_{21}$
	$c_3$	$a_{31}$	$a_{32}$
	$⋮$	$⋮$		$\ddots$
	$c_s$	$a_{s1}$	$a_{s2}$	$\cdots$	$a_{s,s-1}$
		$b_1$	$b_2$	$\cdots$	$b_{s-1}$	$b_s$

Методът на Рунге-Кута е коректен, ако

$\sum _{j=1}^{i-1}{a}_{ij}=c_i\mathrm{f}\mathrm{o}\mathrm{r}i=2,\dots ,s.$

Има и други съпътстващи изисквания, ако искаме метода да има определен ред p, което ознавача, че локалното грешката от закръгляване е O(h^p+1). Това може да бъде извлечено от дефиницията за грешка от закръгляване. Например, 2-етапен метод е от 2-ри ред ако b₁ + b₂ = 1, b₂c₂ = 1/2, и a₂₁ = c₂.

Примери

РК4 методът има следната структура. Неговата таблицата е както следва:

	0
	1/2	1/2
	1/2	0	1/2
	1	0	0	1
		1/6	1/3	1/3	1/6

Най-простият метод на Рунге-Кута е методът на Ойлер, представен с формулата $y_{n+1}=y_n+hf(t_n,y_n)$. Това е единственият коректно поставен експлицитен метод на Рунге-Кута от един етап. Таблицата, която отговаря на тази формула е:

	0
		1

Методи от 2-ри ред

Методът на междинната точка представлява пример за метод от 2-ри ред с два етапа:

$y_{n+1}=y_n+hf(t_n+\frac{1}{2}h,y_n+\frac{1}{2}hf(t_n,y_n)).$

Таблицата е:

	0
	1/2	1/2
		0	1

ТМетодът на междинната точка не е единствен метод на РК от 2-ри ред с два етапа. Всъщност, има семейство от такива методи, параметризирани с α и представени с формула

$y_{n+1}=y_n+h((1-\frac{1}{2\alpha })f(t_n,y_n)+\frac{1}{2\alpha }f(t_n+\alpha h,y_n+\alpha hf(t_n,y_n))).$

Таблицата на Бутчер е:

	0
	$\alpha$	$\alpha$
		$1-\frac{1}{2\alpha }$	$\frac{1}{2\alpha }$

В това семейство при частния случай $\alpha =\frac{1}{2}$ получаваме метода на междинната точка, а при $\alpha =1$ метод на Хойн.

Използване

Като пример, да разгледаме метода на РК от 2-ри ред с 2 етапа с α=2/3. Той се представя с таблицата:

	0
	2/3	2/3
		1/4	3/4

със съответните уравнения:

$\begin{array}{cc}{k}_1& =f(t_n,y_n),\\ k_2& =f(t_n+\frac{2}{3}h,y_n+\frac{2}{3}h{k}_1),\\ y_{n+1}& =y_n+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2).\end{array}$

Методът се използва за да реши началната задача:

$y^{\prime }=\tan (y)+1,y(1)=1,t\in [1,1.1]$

със стъпка h=0.025, следователно методът трябва да направи 4 стъпки.

Изпълнението на метода е както следва:

$t_0=1:$
	$y_0=1$
$t_1=1.025:$
	$y_0=1$	$k_1=2.557407725$	$k_2=f(t_0+\frac{2}{3}h,y_0+\frac{2}{3}h{k}_1)=2.7139$
	$y_1=y_0+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2)=\underset{\_}{1.066869388}$
$t_2=1.05:$
	$y_1=1.066869388$	$k_1=2.813524695$	$k_2=f(t_1+\frac{2}{3}h,y_1+\frac{2}{3}h{k}_1)$
	$y_2=y_1+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2)=\underset{\_}{1.141332181}$
$t_3=1.075:$
	$y_2=1.141332181$	$k_1=3.183536647$	$k_2=f(t_2+\frac{2}{3}h,y_2+\frac{2}{3}h{k}_1)$
	$y_3=y_2+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2)=\underset{\_}{1.227417567}$
$t_4=1.1:$
	$y_3=1.227417567$	$k_1=3.796866512$	$k_2=f(t_3+\frac{2}{3}h,y_3+\frac{2}{3}h{k}_1)$
	$y_4=y_3+h(\frac{1}{4}{k}_1+\frac{3}{4}{k}_2)=\underset{\_}{1.335079087}.$

Числените решения отговарят на подчертаните стойности.

Адаптивни методи на Рунге-Кута

Адаптивните методи са съставени, така че да оценяват локалната грешката от закръгляване на всяка РК-стъпка. Това става използвайки два метода в таблици, единия от ред $p$, а другия от ред $p-1$.

Стъпката от по-нисък ред се получава от:

$y_{n+1}^{*}=y_n+\sum _{i=1}^s{b}_i^{*}{k}_i,$

където $k_i$ са същите като за метод от по-висок ред. Тогава грешката е:

$e_{n+1}=y_{n+1}-y_{n+1}^{*}=h\sum _{i=1}^s(b_i-b_i^{*})k_i,$

което е$O(h^p)$.

Таблицата на Бутчер за тови вид метод е разширена за да се въведат стойностите на $b_i^{*}$:

	0
	$c_2$	$a_{21}$
	$c_3$	$a_{31}$	$a_{32}$
	$⋮$	$⋮$		$\ddots$
	$c_s$	$a_{s1}$	$a_{s2}$	$\cdots$	$a_{s,s-1}$
		$b_1$	$b_2$	$\cdots$	$b_{s-1}$	$b_s$
		$b_1^{*}$	$b_2^{*}$	$\cdots$	$b_{s-1}^{*}$	$b_s^{*}$

Методът на Рунге-Кута-Фелберг използва два метода от 5-и и 4-ти ред. Неговата разширена таблица е:

	0
	1/4	1/4
	3/8	3/32	9/32
	12/13	1932/2197	−7200/2197	7296/2197
	1	439/216	−8	3680/513	-845/4104
	1/2	−8/27	2	−3544/2565	1859/4104	−11/40
		16/135	0	6656/12825	28561/56430	−9/50	2/55
		25/216	0	1408/2565	2197/4104	−1/5	0

Най-простият адаптивен метод на РК включва комбинирането на метод на Хойн, който е от 2-ри ред, с метод на Ойлер, който е от 1-ви ред. Неговата разширена таблица е:

	0
	1	1
		1/2	1/2
		1	0

Оценката на грешката се използва за контролиране на големината на стъпката.

Неявни методи на Рунге-Кута

Всички споменати методи на Рунге-Кута са експлицитни. Експлицитните методи на Рунге-Кута са неподходящи за решаване на така наречените „stiff equations“ (твърди уравнения), защото обхвата на абсолютната им устойчивост е малък, по-точно ограничен. Този проблем е особено важен при решаването на частни диференциални уравнения.

Неустойчивостта на експлицитните методи подбужда създаването на имплицитните методи. Имплицитен метод на РК има формата:

$y_{n+1}=y_n+\sum _{i=1}^s{b}_i{k}_i,$

където

$k_i=hf(t_n+c_{i}h,y_n+\sum _{j=1}^s{a}_{ij}{k}_j),i=1,\dots ,s.$

Разликата с експлицитния метод е, че при експлицитния метод сумата на j стига само до i – 1. Това се вижда и в таблицата. За имплицитен метод, матрицата на коефициенти $a_{ij}$ не е непременно долно-триъгълна:

$\begin{array}{ccccc}{c}_1& a_{11}& a_{12}& \dots & a_{1s}\\ c_2& a_{21}& a_{22}& \dots & a_{2s}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ c_s& a_{s1}& a_{s2}& \dots & a_{ss}\\ & b_1& b_2& \dots & b_s\end{array}=\begin{array}{cc}𝐜& A\\ & {𝐛}^{𝐓}\end{array}$

Последиците от тази разлика са, че на всяка стъпка се решава система от алгебрични уравнения. Това значително повишава изчислителния разход. Ако метод със s етапа се използва за решаване на диференциално уравнение с m компонента, системата от алгебрични уравнения ще има ms компонента. За сравнение, при имплицитен s-стъпков линеен метод се решава система от алгебрични уравнения със само s компонента.

Примери

Най-простият имплицитен метод на РК е представен с метод на Ойлер назад:

$y_{n+1}=y_n+hf(t_n+h,y_{n+1}).$

Таблицата е просто:

$\begin{array}{cc}1& 1\\ & 1\end{array}$

Тази таблица отговаря на формулите:

$k_1=f(t_n+h,y_n+h{k}_1)\text{and}{y}_{n+1}=y_n+h{k}_1,$

които могат да бъдат пренаредени за да се получи формула за метод на Ойлер назад, описан по-горе.

Друг пример за имплицитен метод на Рунге-Кута е правилото на трапеца. Съответната таблица е:

$\begin{array}{ccc}0& 0& 0\\ 1& \frac{1}{2}& \frac{1}{2}\\ & \frac{1}{2}& \frac{1}{2}\end{array}$

Методите на Гаус-Льожандър формират семейство методи, базирани на квадратурата на Гаус. Метод на Гаус-Льожандър със s етапа е от ред 2s. Метод с два етапа (и 4-ти ред) има следната таблица:

$\begin{array}{ccc}\frac{1}{2}-\frac{1}{6}\sqrt{3}& \frac{1}{4}& \frac{1}{4}-\frac{1}{6}\sqrt{3}\\ \frac{1}{2}+\frac{1}{6}\sqrt{3}& \frac{1}{4}+\frac{1}{6}\sqrt{3}& \frac{1}{4}\\ & \frac{1}{2}& \frac{1}{2}\end{array}$

Устойчивост

Предимството на имплицитните методи на Рунге-Кута пред експлицитните е по-голямата им устойчивост, особено когато се прилагат върху твърди уравнения. Разглеждаме линейното тестово уравнение y' = λy. Методът на РК, приложен за това уравнение, свежда решаването му до изпълняване на итерациите $y_{n+1}=r(h\lambda )y_n$, където r е дадено чрез

$r(z)=1+z{b}^T(I-zA{)}^{-1}e=\frac{\det (I-zA+ze{b}^T)}{\det (I-zA)},$

където е е единичния вектор. Функцията r се нарича функция на устойчивостта. Експлицитните методи имат строга ниско-триъгълна матрица А, която предполага, че det(I − zA) = 1 и функцията на устойчивостта е полином.

Численото решение на линейното тестово уравнение се разпада до нула ако |r(z)| < 1 при z=hλ. Множеството от такива стойности на z се нарича област на абсолютна устойчивост. В частност, методът трябва да бъде A-стабилен, ако всяко z с Re(z)<0 се намира в областта на абсолютната устойчивост. Функцията на устойчивостта на експлицитния метод на РК е полином, следователно експлицитните методи на РК никога не могат да бъдат A-стабилни.

Извеждане на метода на Рунге-Кута от 4-ти ред

Методът на Рунге-Кута от ред $s$ може да бъде написан така:

$y_{t+h}=y_t+h\cdot \sum _{i=1}^s{a}_i{k}_i+𝒪(h^{s+1})$

където:

$k_i=f(y_t+h\cdot \sum _{j=1}^s{\beta }_{ij}{k}_j,t_n+{\alpha }_{i}h)$

са нараствания получени чрез пресмятане на производните на $y_t$ в $i$-тия ред.

Извеждането на метода на Рунге-Кута от 4-ти ред се постига чрез използване на общата формула с $s=4$, както е описано по-горе, в началната точка, междинна точка и последната точка на всеки интервал $(t,t+h)$, затова избираме:

${\alpha }_i$	${\beta }_i$
${\alpha }_1=0$	${\beta }_{21}=\frac{1}{2}$
${\alpha }_2=\frac{1}{2}$	${\beta }_{32}=\frac{1}{2}$
${\alpha }_3=\frac{1}{2}$	${\beta }_{43}=1$
${\alpha }_4=1$

и ${\beta }_{ij}=0$. Започваме с определянето на следните величини:

$\begin{array}{cc}{y}_{t+h}^1& =y_t+hf(y_t,t)\\ y_{t+h}^2& =y_t+hf(y_{t+h/2}^1,t+\frac{h}{2})\\ y_{t+h}^3& =y_t+hf(y_{t+h/2}^2,t+\frac{h}{2})\end{array}$

където $y_{t+h/2}^1=\frac{y_t+y_{t+h}^1}{2}$ and $y_{t+h/2}^2=\frac{y_t+y_{t+h}^2}{2}$ Ако задедем:

$\begin{array}{cc}{k}_1& =f(y_t,t)\\ k_2& =f(y_{t+h/2}^1,t+\frac{h}{2})\\ k_3& =f(y_{t+h/2}^2,t+\frac{h}{2})\\ k_4& =f(y_{t+h}^3,t+h)\end{array}$

и имайки предвид горните изрази можем да покажем, че следните равенства имат точност $𝒪(h^2)$:

$\begin{array}{cc}{k}_2& =f(y_{t+h/2}^1,t+\frac{h}{2})=f(y_t+\frac{h}{2}{k}_1,t+\frac{h}{2})\\ & =f(y_t,t)+\frac{h}{2}\frac{d}{dt}f(y_t,t)\\ k_3& =f(y_{t+h/2}^2,t+\frac{h}{2})=f(y_t+\frac{h}{2}f(y_t+\frac{h}{2}{k}_1,t+\frac{h}{2}),t+\frac{h}{2})\\ & =f(y_t,t)+\frac{h}{2}\frac{d}{dt}[f(y_t,t)+\frac{h}{2}\frac{d}{dt}f(y_t,t)]\\ k_4& =f(y_{t+h}^3,t+h)=f(y_t+hf(y_t+\frac{h}{2}{k}_2,t+\frac{h}{2}),t+h)\\ & =f(y_t+hf(y_t+\frac{h}{2}f(y_t+\frac{h}{2}f(y_t,t),t+\frac{h}{2}),t+\frac{h}{2}),t+h)\\ & =f(y_t,t)+h\frac{d}{dt}[f(y_t,t)+\frac{h}{2}\frac{d}{dt}[f(y_t,t)+\frac{h}{2}\frac{d}{dt}f(y_t,t)]]\end{array}$

където:

$\frac{d}{dt}f(y_t,t)=\frac{\partial }{\partial y}f(y_t,t){\dot{y}}_t+\frac{\partial }{\partial t}f(y_t,t)=f_y(y_t,t)\dot{y}+f_t(y_t,t):={\ddot{y}}_t$

е пълната производна на $f$ по отношение на времето t.

Ако изразим общата формула, използвайки това, което изведохме по-горе, получаваме:

$\begin{array}{cc}{y}_{t+h}& =y_t+h\{a\cdot f(y_t,t)+b\cdot [f(y_t,t)+\frac{h}{2}\frac{d}{dt}f(y_t,t)]+\\ & +c\cdot [f(y_t,t)+\frac{h}{2}\frac{d}{dt}[f(y_t,t)+\frac{h}{2}\frac{d}{dt}f(y_t,t)]]+\\ & +d\cdot [f(y_t,t)+h\frac{d}{dt}[f(y_t,t)+\frac{h}{2}\frac{d}{dt}[f(y_t,t)+\frac{h}{2}\frac{d}{dt}f(y_t,t)]]]\}+𝒪(h^5)\\ & =y_t+a\cdot h{f}_t+b\cdot h{f}_t+b\cdot \frac{h^2}{2}\frac{d{f}_t}{dt}+c\cdot h{f}_t+c\cdot \frac{h^2}{2}\frac{d{f}_t}{dt}+\\ & +c\cdot \frac{h^3}{4}\frac{d^2{f}_t}{d{t}^2}+d\cdot h{f}_t+d\cdot h^2\frac{d{f}_t}{dt}+d\cdot \frac{h^3}{2}\frac{d^2{f}_t}{d{t}^2}+d\cdot \frac{h^4}{4}\frac{d^3{f}_t}{d{t}^3}+𝒪(h^5)\end{array}$

И като сравним това с реда на Тейлър за $y_{t+h}$ около $y_t$:

$\begin{array}{cc}{y}_{t+h}& =y_t+h{\dot{y}}_t+\frac{h^2}{2}{\ddot{y}}_t+\frac{h^3}{6}{y}_t^{(3)}+\frac{h^4}{24}{y}_t^{(4)}+𝒪(h^5)=\\ & =y_t+hf(y_t,t)+\frac{h^2}{2}\frac{d}{dt}f(y_t,t)+\frac{h^3}{6}\frac{d^2}{d{t}^2}f(y_t,t)+\frac{h^4}{24}\frac{d^3}{d{t}^3}f(y_t,t)\end{array}$

получаваме система от изисквания за коефициентите:

$\{\begin{array}{cc}& a+b+c+d=1\\ & \frac{1}{2}b+\frac{1}{2}c+d=\frac{1}{2}\\ & \frac{1}{4}c+\frac{1}{2}d=\frac{1}{6}\\ & \frac{1}{4}d=\frac{1}{24}\end{array}$

която решена дава $a=\frac{1}{6},b=\frac{1}{3},c=\frac{1}{3},d=\frac{1}{6}$.

Литература

• Шаблон:Cite book.
• Шаблон:Cite book.
• Шаблон:Cite book.
• Шаблон:Cite book.
• Шаблон:Cite book.
• Шаблон:Cite book (see Chapter 6).
• Шаблон:Cite book.
• Шаблон:Cite book.
• Шаблон:Cite book.
• Шаблон:Cite book.
• Шаблон:Cite book. Also, Section 17.2. Adaptive Stepsize Control for Runge-Kutta Шаблон:Webarchive.
• Шаблон:Cite book.
• Шаблон:Cite book.

Външни препратки

• Runge–Kutta 4th Order Method
• Runge Kutta Method for O.D.E.'s
• DotNumerics: Ordinary Differential Equations for C# and VB.NET – Initial-value problem for nonstiff and stiff ordinary differential equations (explicit Runge–Kutta, implicit Runge–Kutta, Gear's BDF and Adams–Moulton).
• GafferOnGames Шаблон:Webarchive – A physics resource for computer programmers
• PottersWheel – Parameter calibration in ODE systems using implicit Runge–Kutta integration