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

от Уикипедия, свободната енциклопедия
Направо към: навигация, търсене

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

Метод на Рунге-Кута от 4-ти ред[редактиране | edit source]

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

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

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

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

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

\begin{align}
y_{n+1} &= y_n + \tfrac{1}{6} \left(k_1 + 2k_2 + 2k_3 + k_4 \right) \\
t_{n+1} &= t_n + h \\
\end{align}

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


\begin{align}
k_1 &= hf(t_n, y_n),
\\
k_2 &= hf(t_n + \tfrac{1}{2}h , y_n +  \tfrac{1}{2} k_1),
\\
k_3 &= hf(t_n + \tfrac{1}{2}h , y_n +   \tfrac{1}{2} k_2),
\\
k_4 &= hf(t_n + h , y_n + k_3).
\end{align}
[1]
(Забележка: горните уравнения имат различни, но равнозначни дефиниции в различните текстове).[2]

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

  • k_1 е нарастване, основаващо се на наклона от началото на интервала, използвайки  y_n , Метод на (Ойлер) ;
  • k_2 е нарастване, основаващо се на наклона от междинна точка от интервала, използвайки  y_n + \tfrac12 k_1  ;
  • k_3 е отново нарастване, основаващо се на налкона от междинна точка, но използвайки  y_n + \tfrac12 k_2  ;
  • k_4 е нарастване, основаващо се на наклона в края на интервала, използвайки  y_n + k_3 .

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

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

Експлицитни методи на Рунге-Кута[редактиране | edit source]

Семейството на експлицитните методи на Рунге-Кута са обобщение на РК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_2h, y_n+a_{21}k_1), \,
 k_3 = hf(t_n+c_3h, y_n+a_{31}k_1+a_{32}k_2), \,
 \vdots
 k_s = hf(t_n+c_sh, y_n+a_{s1}k_1+a_{s2}k_2+\cdots+a_{s,s-1}k_{s-1}). [4]
(Забележка: горните уравнения имат различни, но равнозначни дефиниции в различните текстове).[2]

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

0
 c_2  a_{21}
 c_3  a_{31}  a_{32}
 \vdots  \vdots  \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{for}\ i=2, \ldots, s.

Има и други съпътстващи изисквания, ако искаме метода да има определен ред p, което ознавача, че локалното грешката от закръгляване е O(hp+1). Това може да бъде извлечено от дефиницията за грешка от закръгляване. Например, 2-етапен метод е от 2-ри ред ако b1 + b2 = 1, b2c2 = 1/2, и a21 = c2.[6]

Примери[редактиране | edit source]

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

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-ри ред с два етапа[редактиране | edit source]

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

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

Таблицата е:

0
1/2 1/2
0 1

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

 y_{n+1} = y_n + h\bigl( (1-\tfrac1{2\alpha}) f(t_n, y_n) + \tfrac1{2\alpha} f(t_n + \alpha h, y_n + \alpha h f(t_n, y_n)) \bigr). [8]

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

0
\alpha \alpha
1-\tfrac1{2\alpha} \tfrac1{2\alpha}

В това семейство при частният случай \alpha=\tfrac12 получаваме метода на междинната точка, а при \alpha=1 метод на Хойн.[3]

Използване[редактиране | edit source]

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

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

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

 \begin{align}
k_1 &= f(t_n,y_n), \\
k_2 &= f(t_n + \tfrac{2}{3}h, y_n + \tfrac{2}{3}h k_1), \\
y_{n+1} &= y_n + h\left(\tfrac{1}{4}k_1+\tfrac{3}{4}k_2\right).
\end{align}

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

 y' = \tan(y)+1,\quad y(1)=1,\ t\in [1, 1.1]

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

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

t_0=1 \colon
y_0=1
t_1=1.025 \colon
y_0 = 1 k_1=2.557407725 k_2 = f(t_0 + \tfrac23h ,y_0 + \tfrac23hk_1) = 2.7139
y_1=y_0+h(\tfrac14k_1 + \tfrac34k_2)=\underline{1.066869388}
t_2=1.05 \colon
y_1 = 1.066869388 k_1=2.813524695 k_2 = f(t_1 + \tfrac23h ,y_1 + \tfrac23hk_1)
y_2=y_1+h(\tfrac14k_1 + \tfrac34k_2)=\underline{1.141332181}
t_3=1.075 \colon
y_2 = 1.141332181 k_1=3.183536647 k_2 = f(t_2 + \tfrac23h ,y_2 + \tfrac23hk_1)
y_3=y_2+h(\tfrac14k_1 + \tfrac34k_2)=\underline{1.227417567}
t_4=1.1 \colon
y_3 = 1.227417567 k_1=3.796866512 k_2 = f(t_3 + \tfrac23h ,y_3 + \tfrac23hk_1)
y_4=y_3+h(\tfrac14k_1 + \tfrac34k_2)=\underline{1.335079087}.

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

Адаптивни методи на Рунге-Кута[редактиране | edit source]

Адаптивните методи са съставени, така че да оценяват локалната грешката от закръгляване на всяка РК-стъпка. Това става използвайки два метода в таблици, единия от ред 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}
 \vdots  \vdots  \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

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

Имплицитни методи на Рунге-Кута[редактиране | edit source]

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

 y_{n+1} = y_n + \sum_{i=1}^s b_i k_i,

където

 k_i = hf\bigl( t_n + c_i h, y_n + \sum_{j=1}^s a_{ij} k_j \bigr), \quad i = 1, \ldots, s. [10]

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


\begin{array}{c|cccc}
c_1    & a_{11} & a_{12}& \dots & a_{1s}\\
c_2    & a_{21} & a_{22}& \dots & a_{2s}\\
\vdots & \vdots & \vdots& \ddots& \vdots\\
c_s    & a_{s1} & a_{s2}& \dots & a_{ss} \\
\hline
       & b_1    & b_2   & \dots & b_s\\
\end{array} = 

\begin{array}{c|c}
\mathbf{c}& A\\
\hline
          & \mathbf{b^T} \\
\end{array}

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

Примери[редактиране | edit source]

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

y_{n + 1} = y_n + h f(t_n + h, y_{n + 1}). \,

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


\begin{array}{c|c}
1 & 1 \\
\hline
  & 1 \\
\end{array}

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

 k_1 = f(t_n + h, y_n + h k_1) \quad\text{and}\quad y_{n+1} = y_n + h k_1,

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

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


\begin{array}{c|cc}
0 & 0& 0\\
1 & \frac{1}{2}& \frac{1}{2}\\
\hline
  &  \frac{1}{2}&\frac{1}{2}\\
\end{array}

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


\begin{array}{c|cc}
\frac12 - \frac16 \sqrt3 & \frac14                  & \frac14 - \frac16 \sqrt3 \\
\frac12 + \frac16 \sqrt3 & \frac14 + \frac16 \sqrt3 & \frac14 \\ 
\hline
                         & \frac12                  & \frac12 
\end{array}
[11]

Устойчивост[редактиране | edit source]

Предимството на имплицитните методи на Рунге-Кута пред експлицитните е по-голямата им устойчивост, особено когато се прилагат върху твърди уравнения. Разглеждаме линейното тестово уравнение 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+zeb^T)}{\det(I-zA)}, [13]

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

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

Извеждане на метода на Рунге-Кута от 4-ти ред[редактиране | edit source]

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

y_{t + h} = y_t + h \cdot \sum_{i=1}^s a_i k_i +\mathcal{O}(h^{s+1})

където:

k_i = f\left(y_t + h \cdot \sum_{j = 1}^s \beta_{ij} k_j, t_n + \alpha_i h \right)

са нараствания получени чрез пресмятане на производните на 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{align}
y^1_{t+h} &= y_t + hf\left(y_t, t\right) \\
y^2_{t+h} &= y_t + hf\left(y^1_{t+h/2}, t+\frac{h}{2}\right) \\
y^3_{t+h} &= y_t + hf\left(y^2_{t+h/2}, t+\frac{h}{2}\right)
\end{align}

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

\begin{align}
k_1 &= f(y_t, t) \\
k_2 &= f\left(y^1_{t+h/2}, t + \frac{h}{2}\right) \\
k_3 &= f\left(y^2_{t+h/2}, t + \frac{h}{2}\right) \\
k_4 &= f\left(y^3_{t+h}, t + h\right)
\end{align}

и имайки в предвид горните изрази можем да покажем, че следните равенства имат точност \mathcal{O}(h^2):

\begin{align}
k_2 &= f\left(y^1_{t+h/2}, t + \frac{h}{2}\right) = f\left(y_t + \frac{h}{2} k_1, t + \frac{h}{2}\right) \\
&= f\left(y_t, t\right) + \frac{h}{2} \frac{d}{dt}f\left(y_t,t\right) \\
k_3 &= f\left(y^2_{t+h/2}, t + \frac{h}{2}\right) = f\left(y_t + \frac{h}{2} f\left(y_t + \frac{h}{2} k_1, t + \frac{h}{2}\right), t + \frac{h}{2}\right) \\
&= f\left(y_t, t\right) + \frac{h}{2} \frac{d}{dt} \left[ f\left(y_t,t\right) + \frac{h}{2} \frac{d}{dt}f\left(y_t,t\right) \right] \\
k_4 &= f\left(y^3_{t+h}, t + h\right) = f\left(y_t + h f\left(y_t + \frac{h}{2} k_2, t + \frac{h}{2}\right), t + h\right) \\
&= f\left(y_t + h f\left(y_t + \frac{h}{2} f\left(y_t + \frac{h}{2} f\left(y_t, t\right), t + \frac{h}{2}\right), t + \frac{h}{2}\right), t + h\right)  \\
&= f\left(y_t, t\right) + h \frac{d}{dt} \left[ f\left(y_t,t\right) + \frac{h}{2} \frac{d}{dt}\left[ f\left(y_t,t\right) + \frac{h}{2} \frac{d}{dt}f\left(y_t,t\right) \right]\right]
\end{align}

където:

\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{align}
y_{t+h} &= y_t + h  \left\lbrace a \cdot f(y_t, t) + b \cdot \left[ f\left(y_t, t\right) + \frac{h}{2} \frac{d}{dt}f\left(y_t,t\right) \right] \right.+ \\
&+ c \cdot \left[ f\left(y_t, t\right) + \frac{h}{2} \frac{d}{dt} \left[ f\left(y_t,t\right) + \frac{h}{2} \frac{d}{dt}f\left(y_t,t\right) \right] \right] + \\
&+ d \cdot \left[f\left(y_t, t\right) + h \frac{d}{dt} \left[ f\left(y_t,t\right) + \frac{h}{2} \frac{d}{dt}\left[ f\left(y_t,t\right) 
+ \left. \frac{h}{2} \frac{d}{dt}f\left(y_t,t\right) \right]\right]\right]\right\rbrace + \mathcal{O}(h^5) \\
&= y_t + a \cdot h f_t + b \cdot h f_t + b \cdot \frac{h^2}{2} \frac{df_t}{dt}   + c \cdot h f_t+ c \cdot \frac{h^2}{2} \frac{df_t}{dt} + \\
&+ c \cdot \frac{h^3}{4} \frac{d^2f_t}{dt^2} + d \cdot h f_t + d \cdot h^2 \frac{df_t}{dt} + d \cdot \frac{h^3}{2} \frac{d^2f_t}{dt^2} + d \cdot \frac{h^4}{4} \frac{d^3f_t}{dt^3} + \mathcal{O}(h^5)
\end{align}

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

\begin{align}
 y_{t+h} &= y_t + h \dot y_t + \frac{h^2}{2} \ddot y_t + \frac{h^3}{6} y^{(3)}_t + \frac{h^4}{24} y^{(4)}_t + \mathcal{O}(h^5) = \\
&= y_t + h f(y_t, t) + \frac{h^2}{2} \frac{d}{dt}f(y_t, t) + \frac{h^3}{6} \frac{d^2}{dt^2}f(y_t, t) + \frac{h^4}{24} \frac{d^3}{dt^3}f(y_t, t) 
\end{align}

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

\left\{
 \begin{align}
 &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{align}\right.

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


Източници[редактиране | edit source]

  1. Press et al. 2007, с. 908
    Süli & Mayers 2003, с. 328
  2. а б Шаблон:Harvtxt, Шаблон:Harvtxt, Шаблон:Harvtxt and Шаблон:Harvtxt leave out the factor h in the definition of the stages. Шаблон:Harvtxt, Шаблон:Harvtxt and Шаблон:Harvtxt use the y-values as stages.
  3. а б Süli & Mayers 2003, с. 328
  4. Press et al. 2007, с. 907
  5. Iserles 1996, с. 38
  6. Iserles 1996, с. 39
  7. а б Süli & Mayers 2003, с. 352
  8. Süli & Mayers 2003, с. 327
  9. Süli & Mayers 2003, с. 349–351
  10. Iserles 1996, с. 41
    Süli & Mayers 2003, с. 351–352
  11. а б Süli & Mayers 2003, с. 353
  12. Iserles 1996, с. 47
  13. Hairer & Wanner 1996, с. 40–41
  14. Hairer & Wanner 1996, с. 40
  15. а б Iserles 1996, с. 60

Литература[редактиране | edit source]

Външни препратки[редактиране | edit source]

Шаблон:Numerical integrators

Криейтив Комънс - Признание - Споделяне на споделеното Лиценз за свободна документация на ГНУ Тази страница частично или изцяло представлява превод на страницата „Runge–Kutta methods“ в Уикипедия на английски. Оригиналният текст, както и този превод, са защитени от Лиценза „Криейтив Комънс - Признание - Споделяне на споделеното“, а за съдържание, създадено преди юни 2009 година — от Лиценза за свободна документация на ГНУ. Прегледайте историята на редакциите на оригиналната страница, както и на преводната страница, за да видите списъка на съавторите.