龙格-库塔方法
M阶Runge-Kutta方法用于求解方程，
d y d x = f ( x , y ) , y ( a ) = y 0 ( 1 ) \frac{\boldsymbol{dy}}{\boldsymbol{dx}}=\boldsymbol{f}\left( \boldsymbol{x},\boldsymbol{y} \right) ,\boldsymbol{y}\left( \boldsymbol{a} \right) =\boldsymbol{y}_0\,\, \left( 1 \right) dxdy=f(x,y),y(a)=y0(1)
其中 y : R → R n , f : R × R ⋉ → R n , x ∈ [ a , b ] \boldsymbol{y}:\mathbb{R}\rightarrow \mathbb{R}^{\boldsymbol{n}},\boldsymbol{f}:\mathbb{R}\times \mathbb{R}^{\ltimes}\rightarrow \mathbb{R}^{\boldsymbol{n}},\boldsymbol{x}\in \left[ \boldsymbol{a},\boldsymbol{b} \right] y:R→Rn,f:R×R⋉→Rn,x∈[a,b] 。在这里，我们假设函数 f j ( x , y ) \boldsymbol{f}_{\boldsymbol{j}}\left( \boldsymbol{x},\boldsymbol{y} \right) fj(x,y) 对于所有 j = 1 , ⋯ , n \boldsymbol{j}=1,\cdots ,\boldsymbol{n} j=1,⋯,n在 [ a , b ] \left[ \boldsymbol{a},\boldsymbol{b} \right] [a,b] 中都是利普希茨连续的，并且至少一个函数 f k ( x , y ) \boldsymbol{f}_{\boldsymbol{k}}\left( \boldsymbol{x},\boldsymbol{y} \right) fk(x,y) 对于某些 k ∈ { 1 , 2 , ⋯ , n } \boldsymbol{k}\in \left\{ 1,2,\cdots ,\boldsymbol{n} \right\} k∈{ 1,2,⋯,n} 是非线性的。
等式（1）是任何以三重态 A = ( a ) i j ∈ R M × M , b = ( b ) i ∈ R M , 和 c = ( c ) i ∈ R M , i , j = 1 , ⋯ , M \boldsymbol{A}=\left( \boldsymbol{a} \right) _{\boldsymbol{ij}}\in \boldsymbol{R}^{\boldsymbol{M}\times \boldsymbol{M}},\boldsymbol{b}=\left( \boldsymbol{b} \right) _{\boldsymbol{i}}\in \boldsymbol{R}^{\boldsymbol{M}},和\boldsymbol{c}=\left( \boldsymbol{c} \right) _{\boldsymbol{i}}\in \boldsymbol{R}^{\boldsymbol{M}},\boldsymbol{i},\boldsymbol{j}=1,\cdots ,\boldsymbol{M} A=(a)ij∈RM×M,b=(b)i∈RM,和c=(c)i∈RM,i,j=1,⋯,M为特征的数值方法 。 三元组 ( A , b , c ) \left( \boldsymbol{A},\boldsymbol{b},\boldsymbol{c} \right) (A,b,c) 满足Butcher阵列
参数 c j , j = 1 , ⋯ , M \boldsymbol{c}_{\boldsymbol{j}},\boldsymbol{j}=1,\cdots ,\boldsymbol{M} cj,j=1,⋯,M使得 0 ⩽ c 1 ⩽ c 2 ⩽ ⋯ ⩽ c M ⩽ 1 0\leqslant \boldsymbol{c}_1\leqslant \boldsymbol{c}_2\leqslant \cdots \leqslant \boldsymbol{c}_{\boldsymbol{M}}\leqslant 1 0⩽c1⩽c2⩽⋯⩽cM⩽1。参数 a i j , i , j = 1 , ⋯ , M \boldsymbol{a}_{\boldsymbol{ij}},\boldsymbol{i},\boldsymbol{j}=1,\cdots ,\boldsymbol{M} aij,i,j=1,⋯,M满足
∑ j = 1 M a i j = c i \sum_{\boldsymbol{j}=1}^{\boldsymbol{M}}{\boldsymbol{a}_{\boldsymbol{ij}}=\boldsymbol{c}_{\boldsymbol{i}}} j=1∑Maij=ci
参数 b i , i = 1 , ⋯ , M \boldsymbol{b}_{\boldsymbol{i}},\boldsymbol{i}=1,\cdots ,\boldsymbol{M} bi,i=1,⋯,M满足
∑ i = 1 M b i = 1 \sum_{\boldsymbol{i}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{i}}=1} i=1∑Mbi=1
下面描述用于求解上述等式（1）的M阶段Runge-Kutta方法。
将时间间隔 [ a , b ] \left[ \boldsymbol{a},\boldsymbol{b} \right] [a,b] 分为 N \boldsymbol{N} N个子间隔 [ x i , x i + 1 ] , i = 0 , ⋯ , N \left[ \boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{x}_{\boldsymbol{i}+1} \right] ,\boldsymbol{i}=0,\cdots ,\boldsymbol{N} [xi,xi+1],i=0,⋯,N，每个子长度为 h i = x i + 1 − x i \boldsymbol{h}_{\boldsymbol{i}}=\boldsymbol{x}_{\boldsymbol{i}+1}-\boldsymbol{x}_{\boldsymbol{i}} hi=xi+1−xi。 然后， a = x 0 < x 1 < ⋯ < x N = b \boldsymbol{a}=\boldsymbol{x}_0<\boldsymbol{x}_1<\cdots <\boldsymbol{x}_{\boldsymbol{N}}=\boldsymbol{b} a=x0<x1<⋯<xN=b。
给定 y ( x i ) \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}} \right) y(xi) ，我们通过在 [ x i , x i + 1 ] \left[ \boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{x}_{\boldsymbol{i}+1} \right] [xi,xi+1] 上积分方程（1）来计算 y ( x i + 1 ) \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}+1} \right) y(xi+1) ，如下
y ( x i + 1 ) = y ( x i ) + ∫ x i x i + 1 f ( x , y ( x ) ) d x , i = 0 , ⋯ , N − 1 ( 2 ) \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}+1} \right) =\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}} \right) +\int_{\boldsymbol{x}_{\boldsymbol{i}}}^{\boldsymbol{x}_{\boldsymbol{i}+1}}{\boldsymbol{f}\left( \boldsymbol{x},\boldsymbol{y}\left( \boldsymbol{x} \right) \right) \boldsymbol{dx}},\boldsymbol{i}=0,\cdots ,\boldsymbol{N}-1 \left( 2 \right) y(xi+1)=y(xi)+∫xixi+1f(x,y(x))dx,i=0,⋯,N−1(2)
在每个子间隔 [ x i , x i + 1 ] \left[ \boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{x}_{\boldsymbol{i}+1} \right] [xi,xi+1] 上，令 x i ( j ) = x i + c j h i \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)}=\boldsymbol{x}_{\boldsymbol{i}}+\boldsymbol{c}_{\boldsymbol{j}}\boldsymbol{h}_{\boldsymbol{i}} xi(j)=xi+cjhi，其中 j = 1 , ⋯ , M \boldsymbol{j}=1,\cdots ,\boldsymbol{M} j=1,⋯,M。 然后将斜率评估为 f ( x i ( j ) , y ( x i ( j ) ) ) \boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) f(xi(j),y(xi(j))) ，在子区间 [ x i ( 1 ) , x i ( j ) ] ; j = 1 , ⋯ , M \left[ \boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)},\boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right] ;\boldsymbol{j}=1,\cdots ,\boldsymbol{M} [xi(1),xi(j)];j=1,⋯,M上的斜率加权总和为
f ( x i ( j ) , y ( x i ( j ) ) ) ≈ ∑ l = 1 M a j l f ( x i ( l ) , y ( x i ( l ) ) ) \boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) \approx \sum_{\boldsymbol{l}=1}^{\boldsymbol{M}}{\boldsymbol{a}_{\boldsymbol{jl}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)} \right) \right)} f(xi(j),y(xi(j)))≈l=1∑Majlf(xi(l),y(xi(l)))
其中 ∑ l = 1 M a j l = c j \sum\nolimits_{\boldsymbol{l}=1}^{\boldsymbol{M}}{\boldsymbol{a}_{\boldsymbol{jl}}=\boldsymbol{c}_{\boldsymbol{j}}} ∑l=1Majl=cj，然后，
y ( x i ( j ) ) = y ( x i ( 1 ) ) + ∫ x i ( 1 ) x i ( j ) f ( x , y ( x ) ) d x ≈ y ( x i ( 1 ) ) + h i ∑ l = 1 M a j l f ( x i ( l ) , y ( x i ( l ) ) ) ( 3 ) \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) =\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)} \right) +\int_{\boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)}}^{\boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)}}{\boldsymbol{f}\left( \boldsymbol{x},\boldsymbol{y}\left( \boldsymbol{x} \right) \right) \boldsymbol{dx}}\\\approx \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)} \right) +\boldsymbol{h}_{\boldsymbol{i}}\sum_{\boldsymbol{l}=1}^{\boldsymbol{M}}{\boldsymbol{a}_{\boldsymbol{jl}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)} \right) \right)}\,\, \left( 3 \right) y(xi(j))=y(xi(1))+∫xi(1)xi(j)f(x,y(x))dx≈y(xi(1))+hil=1∑Majlf(xi(l),y(xi(l)))(3)
则在区间 [ x i , x i + 1 ] \left[ \boldsymbol{x}_{\boldsymbol{i}},\boldsymbol{x}_{\boldsymbol{i}+1} \right] [xi,xi+1] 上的斜率的加权平均值是斜率 f ( x i ( j ) , y ( x i ( j ) ) ) ; j = 1 , ⋯ , M \boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) ;\boldsymbol{j}=1,\cdots ,\boldsymbol{M} f(xi(j),y(xi(j)));j=1,⋯,M的凸组合，由此得出，
∑ l = 1 M b j f ( x i ( j ) , y ( x i ( j ) ) ) , ∑ i = 1 M b i = 1 \sum_{\boldsymbol{l}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{j}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) ,\sum_{\boldsymbol{i}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{i}}}}=1 l=1∑Mbjf(xi(j),y(xi(j))),i=1∑Mbi=1
等式（2）最终为
y ( x i ( M ) ) ≈ y ( x i ( 1 ) ) + h i ∑ j = 1 M b j f ( x i ( j ) , y ( x i ( j ) ) ) ( 4 ) \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{M} \right)} \right) \approx \boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( 1 \right)} \right) +\boldsymbol{h}_{\boldsymbol{i}}\sum_{\boldsymbol{j}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{j}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) \,\, \left( 4 \right)} y(xi(M))≈y(xi(1))+hij=1∑Mbjf(xi(j),y(xi(j)))(4)
公式（3）和（4）定义了M阶龙格库塔法。 如果所有 l > j \boldsymbol{l}>\boldsymbol{j} l>j的系数 a j l = 0 \boldsymbol{a}_{\boldsymbol{jl}}=0 ajl=0，则得出的Runge-Kutta方法是显式的。 如果对于所有 l > j \boldsymbol{l}>\boldsymbol{j} l>j， a j l = 0 \boldsymbol{a}_{\boldsymbol{jl}}=0 ajl=0，但对于某些 j \boldsymbol{j} j， a j l ≠ 0 \boldsymbol{a}_{\boldsymbol{jl}}\ne 0 ajl=0，则该方法为半显式。 否则，它是隐式的。
ξ i = y ( x i ( l ) ) + h i ∑ j = 1 M b j f ( x i ( j ) , y ( t i ( j ) ) ) − y ( x i ( M ) ) \boldsymbol{\xi }_{\boldsymbol{i}}=\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{l} \right)} \right) +\boldsymbol{h}_{\boldsymbol{i}}\sum\nolimits_{\boldsymbol{j}=1}^{\boldsymbol{M}}{\boldsymbol{b}_{\boldsymbol{j}}\boldsymbol{f}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)},\boldsymbol{y}\left( \boldsymbol{t}_{\boldsymbol{i}}^{\left( \boldsymbol{j} \right)} \right) \right) -\boldsymbol{y}\left( \boldsymbol{x}_{\boldsymbol{i}}^{\left( \boldsymbol{M} \right)} \right)} ξi=y(xi(l))+hi∑j=1Mbjf(xi(j),y(ti(j)))−y(xi(M))定义了局部截断误差M阶龙格库塔法。 如果M阶龙格库塔法的局部截断误差为 O ( h P + 1 ) \mathcal{O}\left( \boldsymbol{h}^{\boldsymbol{P}+1} \right) O(hP+1) ，其中 h = max i = 0 , ⋯ , N − 1 { h i } \boldsymbol{h}=\underset{\boldsymbol{i}=0,\cdots ,\boldsymbol{N}-1}{\max}\left\{ \boldsymbol{h}_{\boldsymbol{i}} \right\} h=i=0,⋯,N−1max{ hi} ，则M级Runge-Kutta方法的应该为 P ∈ N \boldsymbol{P}\in \mathbb{N} P∈N阶 。
通过满足系数 ( a ) i j , b i \left( \boldsymbol{a} \right) _{\boldsymbol{ij}},\boldsymbol{b}_{\boldsymbol{i}} (a)ij,bi和 c i , i , j ∈ { 1 , ⋯ , M } \boldsymbol{c}_{\boldsymbol{i}},\boldsymbol{i},\boldsymbol{j}\in \left\{ 1,\cdots ,\boldsymbol{M} \right\} ci,i,j∈{ 1,⋯,M} 的条件，可以实现M级龙格库塔法的阶数。 Butcher指出隐式M阶龙格库塔法的条件数，如下表1所示，为P阶。
表1
在表2中，我们陈述了M阶龙格库塔法的系数 ( A , b , c ) \left( \boldsymbol{A},\boldsymbol{b},\boldsymbol{c} \right) (A,b,c) 的条件，最高可达四阶。 如果满足从1到P的所有条件，则该方法的阶数可以为P。
表2
显式龙格库塔法
隐式龙格库塔法
MATLAB解算器
Python解算器
详情参阅 – 亚图跨际
