Phương pháp Runge-Kutta

Trong giải tích số, các phương pháp Runge-Kutta là một họ của các phương pháp lặp ẩn (implicit) và hiện (explicit), trong đó bao gồm thường trình nổi tiếng được gọi là các phương pháp Euler, được sử dụng trong việc rời rạc hóa thời gian để tìm lời giải gần đúng cho các phương trình vi phân thường. Những phương pháp này được phát triển vào khoảng năm 1900 bởi hai nhà toán học người Đức C. Runge và M. W. Kutta.

Xem bài viết về các phương pháp số áp dụng cho việc tìm lời giải của các phương trình vi phân thường để hiểu hơn về nền tảng cơ sở và các phương pháp số khác. Xem thêm Danh sách các phương pháp Runge-Kutta.

Phương pháp Runge-Kutta
Thành viên được biết đến rộng rãi nhất của họ Runge-Kutta là "RK4", "phương pháp Runge-Kutta cổ điển" hoặc đơn giản là "phương pháp Runge-Kutta".

Cho một bài toán giá trị ban đầu được chỉ rõ như sau:
: \dot{y} = f(t, y), \quad y(t_0) = y_0.

Ở đây y là một hàm chưa biết (vô hướng hoặc vector) của thời gian t, mà chúng ta muốn tìm lời giải gần đúng; chúng ta biết rằng \dot{y}, tốc độ thay đổi của y, là một hàm của t và y. Tại thời điểm ban đầu t_0 giá trị y tương ứng là y_0. Hàm f và các dữ liệu t_0, y_0 được cho trước.

Bây giờ chọn một bước kích thước h> 0 và định nghĩa

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

với n = 0, 1, 2, 3,..., sử dụng
:
\begin{align}
k_1 &= f(t_n, y_n), \\
k_3 &= f(t_n + \frac{h}{2}, y_n + \frac{h}{2} k_2), \\
k_4 &= f(t_n + h, y_n + hk_3).
\end{align}

:(Lưu ý: các phương trình trên có những định nghĩa khác nhau nhưng tương đương trong các văn bản khác nhau).

Ở đây y_{n+1} là giá trị gần đúng (xấp xỉ) RK4 của y(t_{n+1}), và giá trị tiếp theo (y_{n+1}) được xác định bằng giá trị hiện tại (y_n) cộng với trung bình trọng lượng của bốn số gia, trong đó mỗi số gia là sản phẩm của kích cỡ của khoảng thời gian, h, và một độ dốc ước tính được chỉ rõ bởi hàm f ở phía bên tay phải của phương trình vi phân.
k_1 là số gia dựa trên độ dốc tại điểm bắt đầu của khoảng thời gian, sử dụng y (phương pháp Euler);
k_2 là số gia dựa trên độ dốc tại điểm giữa của khoảng thời gian, sử dụng y + \tfrac{h}{2}k_1 ;
k_3 cũng là số gia dựa trên độ dốc ở điểm giữa, nhưng sử dụng y + \tfrac{h}{2}k_2 ;
k_4 là số gia dựa trên độ dốc tại điểm cuối của khoảng thời gian, sử dụng y + hk_3 .

Trong việc trung bình của bốn số gia, trọng lượng lớn hơn được trao cho các số gia tại điểm giữa. Nếu f độc lập với y, để các phương trình vi phân tương đương với một tích phân đơn giản, thì RK4 là quy tắc Simpson.

Phương pháp RK4 là một phương pháp bậc 4, có nghĩa rằng sai số cắt cụt cục bộ bậc O(h^5), trong khi tổng lỗi tích lũy là bậc O(h^4).

Các phương pháp Runge-Kutta hiện
Họ các phương pháp Runge-Kutta hiện (explicit) là một sự tổng quát hóa của phương pháp RK4 đã đề cập bên trên. Nó được cho bởi công thức:
: y_{n+1} = y_n + h \sum_{i=1}^s b_i k_i,
trong đó
:
\begin{align}
k_1 & = f(t_n, y_n), \\
k_2 & = f(t_n+c_2h, y_n+h(a_{21}k_1)), \\
k_3 & = f(t_n+c_3h, y_n+h(a_{31}k_1+a_{32}k_2)), \\
& \ \ \vdots \\
k_s & = f(t_n+c_sh, y_n+h(a_{s1}k_1+a_{s2}k_2+\cdots+a_{s,s-1}k_{s-1})).
\end{align}

:(Lưu ý: các phương trình bên trên có các định nghĩa khác nhau nhưng tương đương trong các văn bản khác nhau). Những dữ liệu này thường được sắp xếp trong một thiết bị ghi nhớ, được biết đến như là một hoạt cảnh (tableau) Butcher (đặt theo tên của John C. Butcher):

:

Phương pháp Runge-Kutta là nhất quán nếu
:\sum_{j=1}^{i-1} a_{ij} = c_i \text{ for } i=2, \ldots, s.
Ngoài ra còn có các yêu cầu kèm theo nếu yêu cầu phương pháp phải có bậc p nào đó, có nghĩa rằng sai số cắt cụt cục bộ là O(hp+1). Những điều này có thể được rút ra từ định nghĩa của chính sai số cắt cụt. Cho ví dụ, một phương pháp hai giai đoạn có bậc 2 nếu b1 + b2 = 1, b2c2 = 1/2, và a21 = c2.

Nhìn chung, nếu một phương pháp Runge-Kutta hiện s-giai đoạn có s \ge p, và nếu p \ge 5, thì s > p. s nhỏ nhất cần thiết cho một phương pháp Runge-Kutta hiện s-giai đoạn để có bậc p là một vấn đề mở. Một số giá trị đã biết:
:
\begin{array}{c|cccccccc}
p & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\
\hline
\min s & 1 & 2 & 3 & 4 & 6 & 7 & 9 & 11
\end{array}

Ví dụ
Phương pháp RK4 trong khuôn khổ này. Hoạt cảnh của nó là

:

Một sự thay đổi đôi chút của phương pháp Runge-Kutta cũng là do Kutta tạo ra năm 1901 và được gọi là quy tắc-3/8. Ưu điểm chính phương pháp này là hầu như tất cả các hệ số lỗi là nhỏ hơn so với phương pháp thông thường, nhưng nó đòi hỏi nhiều FLOPS hơn đôi chút (các thao tác điểm-nổi) cho mỗi bước thời gian. Hoạt cảnh Butcher của nó là:

:

Tuy nhiên, phương pháp Runge-Kutta đơn giản nhất là phương pháp Euler (tiếp tới = forward), được cho bởi công thức y_{n+1} = y_n + hf(t_n, y_n) . Đây là phương pháp Runge-Kutte hiện nhất quán duy nhất với một giai đoạn. Hoạt cảnh tương ứng là

:

Các phương pháp bậc-hai với hai giai đoạn
Một ví dụ của một phương pháp bậc hai với hai giai đoạn là phương pháp điểm giữa:
: y_{n+1} = y_n + hf\left(t_n+\frac{1}{2}h, y_n+\frac{1}{2}hf(t_n,\ y_n)\right).
Hoạt cảnh tương ứng là

:

Phương pháp điểm giữa không phải là phương pháp Runge-Kutta bậc-hai duy nhất với hai giai đoạn; có một họ những phương pháp như vậy, được tham số hóa bởi α và cho bởi công thức

: 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).

Hoạt cảnh Butcher của nó là

:

Trong họ này, \alpha=\tfrac12 đối với phương pháp điểm giữa, và \alpha=1 đối với phương pháp Heun. nếu tất cả c_i,\,i=1,2,\ldots,s đều khác nhau.

Các phương pháp Runge-Kutta ẩn
Tất cả các phương pháp Runge-Kutta đã đề cập cho đến bây giờ đều là các phương pháp hiện (explicit). Các phương pháp Runge-Kutta hiện nói chung không phù hợp cho việc tìm lời giải của các phương trình cứng vì vùng ổn định tuyệt đối của chúng rất nhỏ; cụ thể là, bị bao bọc.
Vấn đề này đặc biệt quan trọng trong việc tìm lời giải cho các phương trình vi phân từng phần.

Sự bất ổn định của phương pháp Runge-Kutta hiện thúc đẩy sự phát triển của các phương pháp ẩn (implicit). Một phương pháp Runge-Kutta ẩn có dạng
: y_{n+1} = y_n + h \sum_{i=1}^s b_i k_i,
trong đó
: k_i = f\left(t_n + c_i h,\ y_{n} + h \sum_{j=1}^s a_{ij} k_j \right), \quad i = 1, \ldots, s.

Sự khác nhau so với phương pháp hiện là trong một phương pháp hiện, phép tính tổng đối với chỉ số j chỉ lên đến i − 1. Điều này cũng thấy rõ trong các hoạt cảnh Butcher: ma trận hệ số a_{ij} của một phương pháp hiện là ma trận tam giác dưới Trong một phương pháp ẩn, phép tổng đối với j tính đến s và ma trận hệ số không phải là ma trận tam giác, cho kết quả là bảng Butcher có dạng

Ví dụ
Ví dụ đơn giản của một phương pháp Runge-Kutta implicit là phương pháp Euler lùi (backward):

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

Hoạt cảnh Butcher chỉ đơn giản là:

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

Hoạt cảnh Butcher này tương ứng với các công thức

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

có thể được sắp xếp lại để có được công thức cho phương pháp Euler lùi được liệt kê ở trên.

Một ví dụ nữa của phương pháp Runge-Kutta ẩn là quy tắc hình thang. Hoạt cảnh Butcher của nó là:

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

Quy tắc hình thang là một phương pháp sắp đặt theo thứ tự (collocation) (như thảo luận trong bài viết về phương pháp này). Tất cả các phương pháp sắp đặt theo thứ tự đều là các phương pháp Runge-Kutta ẩn, nhưng không phải tất cả các phương pháp Runge-Kutta ẩn là các phương pháp sắp đặt theo thứ tự.

Các phương pháp Gauss-Legendre tạo thành một họ của các phương pháp sắp đặt theo thứ tự dựa trên phép cầu phương Gauss quadrature. Một phương pháp Gauss-Legendre với s giai đoạn s có bậc 2s (như vậy, các phương pháp với bậc cao tùy ý có thể được xây dựng). Phương pháp với hai giai đoạn (và do đó có bậc bốn) có hoạt cảnh Butcher:
:
\begin{array}{c|cc}
\frac12 - \frac16 \sqrt3 & \frac14 & \frac14 - \frac16 \sqrt3 \\
\frac12 + \frac16 \sqrt3 & \frac14 + \frac16 \sqrt3 & \frac14 \\
\hline
& \frac12 & \frac12 \\
& \frac12+\frac12 \sqrt3 & \frac12-\frac12 \sqrt3
\end{array}


trong đó e là viết tắt của vector với các thành phần đều bằng một. Hàm r được gọi là hàm tính ổn định. Từ công thức, ta có r là thương số của hai đa thức bậc s nếu phương pháp có s giai đoạn. Các phương pháp hiện có ma trận A là ma trận tam giác thấp với các thành phần bằng không dọc đường chéo, tức là det(I − zA) = 1 và rằng hàm tính ổn định là một đa thức.

Lời giải số cho phương trình thử nghiệm tuyến tính phân rã về không nếu | r(z) | y'=\lambda y. Dahlquist đã đề xuất việc nghiên cứu độ ổn định của các lược đồ số (numerical schemes) khi áp dụng cho các hệ phi tuyến thỏa mãn một điều kiện đơn điệu. Các khái niệm tương ứng được xác định như là độ ổn định-G đối với các phương pháp đa bước (và các phương pháp một-chân liên quan) và độ ổn định-B (Butcher, 1975) cho các phương pháp Runge-Kutta. Một phương pháp Runge-Kutta được áp dụng cho hệ phi tuyến tính y'=f(y), với \langle f(y)-f(z),\ y-z \rangle, được gọi là ổn định-B, nếu điều kiện này dẫn đến \|y_{n+1}-z_{n+1}\|\leq\|y_{n}-z_{n}\|.

Cho B, M và Q là ba ma trận s\times s được xác định bởi

: B=\operatorname{diag}(b_1,b_2,\ldots,b_s),\, M=BA+A^TB-bb^T,\, Q=BA^{-1}+A^{-T}B-A^{-T}bb^TA^{-1}.

Một phương pháp Runge-Kutta được gọi là ổn định về mặt đại số nếu các ma trận B and M đều xác định không âm. Điều kiện đủ cho độ ổn định-B là: B and Q là xác định không âm.

Nguồn gốc của phương pháp Runge-Kutta bậc bốn
Nói chung một phương pháp Runge-Kutta bậc 4 có thể được viết như sau:
:y_{t + h} = y_t + h \cdot \sum_{i=1}^s a_i k_i +\mathcal{O}(h^{s+1}),
trong đó:
:k_i = f\left(y_t + h \cdot \sum_{j = 1}^s \beta_{ij} k_j,\ t_n + \alpha_i h \right)
là các số gia đạt được đánh giá các đạo hàm của y_t tại bậc thứ i.

Chúng ta rút ra phương pháp Runge-Kutta bậc bốn bằng cách áp dụng công thức tổng quát với trường hợp s=4, như đã giải thích ở trên, tại điểm khởi đầu, điểm giữa và điểm cuối của bất kỳ khoảng (t,\ t +h), do đó chúng ta chọn:

:
\begin{align}
&\alpha_i & &\beta_{ij} \\
\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 & &\\
\end{align}

và mặt khác \beta_{ij} = 0. Chúng ta bắt đầu bằng cách định nghĩa các đại lượng sau đây:
:\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}
trong đó y^1_{t+h/2} = \dfrac{y_t + y^1_{t+h}}{2} và y^2_{t+h/2} = \dfrac{y_t + y^2_{t+h}}{2}
Nếu chúng ta định nghĩa:
:\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}
và với các mối quan hệ trước đó chúng ta có thể thấy rằng các đẳng thức sau đúng đến bậc \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}
trong đó:
:\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
là đạo hàm toàn phần của f theo thời gian.

Nếu bây giờ chúng ta biểu diễn công thức tổng quát bằng những gì chúng ta vừa rút ra, chúng ta có được:
:\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}

và so sánh với chuỗi Taylor của y_{t+h} xung quanh 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}

chúng ta có được một hệ phương trình của các hệ số:

:
\begin{cases}
&a + b + c + d = 1 \\[6pt]
& \frac{1}{2} b + \frac{1}{2} c + d = \frac{1}{2} \\[6pt]
& \frac{1}{4} c + \frac{1}{2} d = \frac{1}{6} \\[6pt]
& \frac{1}{4} d = \frac{1}{24}
\end{cases}
giải hệ trên ta có a = \frac{1}{6}, b = \frac{1}{3}, c = \frac{1}{3}, d = \frac{1}{6} như đã nói trên.

Xem thêm
Phương pháp Euler
Danh sách các phương pháp Runge–Kutta
Phân tích phương trình vi phân thường bằng giải tích số
PottersWheel – Hiệu chỉnh tham số trong các hệ phương trình vi phân thường (ODE) bằng cách sử dụng tích phân Runge-Kutta ẩn
Phương pháp Runge–Kutta (SDE)
Các phương pháp tuyến tính tổng quát

Tham khảo
Chú thích

Sách
.
.
.
.
.
.
.
(see Chapter 6).
.
.
.

.
.
. Also, Section 17.2. Adaptive Stepsize Control for Runge-Kutta .
.
.
.

Liên kết ngoài

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).