其中L被称为Lipschitz常数
定义3如果一个常微分系统的Lipschitz常数L很大（大于20），则它是刚性的。
1.2隐式RK方法的研究意义
在常微分方程及常微分方程组的数值解法中，Runge–Kutta方法是目前应用最为广泛的数值解法之一，同时又具有误差小，精度高的特点。
尽管显式Runge–Kutta方法能够非常准确、快速的给出大部分常微分方程组的数值解。
但是在化学、自动控制电力系统等领域中，会出现一些病态的常微分方程组，也就是刚性方程组。
刚性方程组对于数值解法的稳定性要求苛刻，比如方程组
y1=-0.01y1-99.99y2,y10=2y2=-100y2,y20=1
将其表示为矩阵形式：
y1y2=-0.01-99.990-100y1y2
令
A=-0.01-99.990-100
发现A特征值为：
λ1=-100，λ2=-0.01，刚性比s=λ1λ2=10000≫1。
方程组的解为：
y1x=e-100x+e-0.01xy2x=e-100x
快瞬态慢瞬态
解由快瞬态和慢瞬态两部分构成。
由于慢瞬态的部分，y1x衰减变得十分缓慢。
当自变量变到x=391时，函数值还未下降到初值的1%，求解区间至少取为（0,391）。
另一方面，由于快瞬态的部分，y2x衰减的非常快，因此步长要取得非常小。
从绝对稳定性的方面来看，如果用四阶显式经典RK方法求解，绝对稳定区间要求λhϵ（-2.78,0），则要求h<0.0278。
这样，在（0,391）上就要计算14065步，计算量巨大，因此计算区间（0,391）内的解时，舍入误差积累会特别严重。
例如取求解区间为0,1，用不同步长h来计算y11和y21的值。
利用四阶显式经典RK方法求解如下：
h
y11
y21
0.04
2.9802322e+17
2.9802322e+17
0.02
9.9004983e-01
1.3929556e-24
0.01
9.9004983e-01
2.5300364e-43
0.001
9.9004983e-01
3.7204130e-44
0.0001
9.9004983e-01
3.7200760e-44
真值
9.9004983e-01
3.7200760e-44
很显然，保留八位有效数字的情况下，要保持良好的精度，步长要取得非常小，这就增加了计算量。
而随着步长的加大，误差也会越来越大。
当步长加大到绝对稳定与之外时（即h>0.0278），计算结果就完全不可信了。
对于刚性方程组，显式方法已经远远不能胜任了，一般采用绝对稳定性更好的方法（如隐式Runge–Kutta方法）进行求解，本文采用单步隐式Runge–Kutta方法，而对于隐式方法中的级值得求解，本文采用Newton迭代法。
1.3RK方法的研究现状
研究基于标准模型方程的Runge–Kutta方法的常见形式为：
yi+1=yi+hj=1rαjkjk1=fxi,yikj=f（xi+λih,yi+hl=1j-1μjlkl）（j=2,3,…r）
（显式）
yi+1=yi+hj=1rαjkjkj=f（xi+λih,yi+hl=1rμjlkl）（j=1,2,3,…r）
（隐式）
因为Runge–Kutta方法是比较成熟的常微分方程数值解法。
所以如今主要是对于经典的Runge–Kutta方法进行完善和扩充，在一定的条件下，提高级数以提高精度。
或者是将Runge–Kutta方法与某些领域结合使用。
在1994年，费景高[1]给出了一种显式Runge–Kutta并行方法，从而实现Runge–Kutta方法在多处理机上的应用。
1997年，Enenkel和Jackson[2]实现了Runge–Kutta方法的对角隐式并行改进。
1999年，廖文远和李庆扬[5]给出了一类求解刚性常微分方程的半隐式多步Runge–Kutta方法。
2000年，张诚坚和余洪兵[3]针对非线性延迟系统构造了一类并行预校算法，给出其算法的局部误差估计，数值实验表明该算法是有效的，且具一定的可比性。
2003年，李爱雄[4]通过对传统单支方法的计算格式进行改造，得到了解DDEs的两类单支并行算法单支并行预校算法和单支并行算法，并对方法的收敛性和稳定性做出了分析。
2008年，庞丽君和朱永忠[6]给出了一类随机微分方程Runge–Kutta方法的指数稳定性。
2.单步RK方法的收敛性和稳定性
2.1单步RK方法的收敛性
对于常微分初值问题
y=fx,ya≤x≤bya=η
（1）
的单步显式公式
yi+1=yi+hφ（xi,yi,h）i=0,1,⋯,n-1y0=η
（2）
局部截断误差可以表示为
Ri+1=yxi+1-yxi+hφxi,yi,hi=0,1,⋯,n-13
定理2[16]：
设y（x）为
（1）的解，yii=0n为
（2）的解。
如果：
（1）存在常数c，使得
Ri+1≤chp+1（i=0,1,2,⋯,n-1）
（2）存在a>0，使得
maxx,y∈Dσ0≤a≤h∂φ（x,y,h）∂y≤L
其中Dσ=x,y|a≤x≤b,yx-σ≤y≤yx+σ
记c0=cL[eLb-a-1]，则当h≤mina,pσc0时，有
E（h）≤c0hp
证:
由（3）得
yxi+1=yxi+hφxi,y（xi）,h+Ri+1（i=0,1,⋯,n-1）（4）
将（4）与
（2）相减
yxi+1-yxi=yxi-yi+hφxi,yxi,h-φxi,yi,h+
Ri+1i=0,1,⋯,n-1
由yx0-y0=0知道，在i=0时，yxi-yi≤c0hp成立。
现在假设0≤i≤k-1时也是成立的。
在h≤pσc0时有：
yxi-yi≤c0pσc0p=σ（i=0,1,⋯,n-1）
进而
φxi,yxi,h-φxi,yi,h
=∂φxi,ηi,h∂yyxi-yi
≤Lyxi-yi0≤i≤k-1
其中ηi是yxi和yi之间的数。
于是定理结合条件与（4）式，可以得到
yxi+1-yi+1
≤yxi-yi+hφxi,yxi,h-φxi,yi,h+Ri+1
≤yxi-yi+Lhyxi-yi+chp+1
=1+Lhyxi-yi+chp+10≤i≤k-1
从而
y（xk-yk）≤1+Lhkyx0-y0+1+Lhk-11+Lh-1chp+1
因为yx0-y0=0及1+Lhk≤eLkh≤eL（b-a），得到
yxk-yk≤cLeLb-a-1hp=c0hp
所以当i=k时yxi-yi≤c0hp也是成立的。
（证毕）
对于单步显式格式而言，当f（x,y）和∂f（x,y）∂y在Dσ内连续时，那么定理1的条件
（2）总是满足的。
所以满足上述条件的单步显式公式，只要也满足相容性条件
φx,yx,0=f（x,y（x））
那么它一定是收敛的。
2.2单步RK方法的稳定性
定义4对于初值问题
（1）,若yii=0n是由
（2）得到，zii=0n是下面扰动问题的解：
zi+1=zi+hφxi,zi,h+δi+1（i=0,1,2,⋯,n-1）z0=η+δ0
如果存在正常数C，ε0及h0，使所有的ε∈（0,ε0]，h∈（0,h0]，当max0≤i≤nδi≤ε时，有
max0≤i≤nyi-zi≤Cε
就称
（2）是稳定的或者零稳定的。
定理3[16]在定理1的条件下，单步显式公式
（2）是稳定的。
3.三类隐式RK方法
3.1引言
对于初值问题
（1），隐式Runge–Kutta方法的一般形式为
yi+1=yi+hj=1rbjkj（5）kj=fxi+cih,yi+hl=1rajlklj=1,2,3,…r（6）
其中，xi=x0+ih，n=0,1,2,⋯，为数轴上的离散点列；h为步长，yi为解y（xi）的近似值；c1，c2，…，cr称为隐式Runge–Kutta方法的步长；b1，b2，…，br为权系数。
令A=ajl，j=1,2,3,…r。
称A为系数矩阵，因此上式可以简化表示为
c1⋮c2
a11…a1r⋮⋱⋮ar1…arr
b1…br
使用Butcher阵的记号，上表可以表示为
c
A
bT
可以看到，如果A是一个主对角元素均为零的下三角矩阵，那么上表就可以表示一个显式Runge–Kutta方法。
如果A是一个主对角线非零的下三角矩阵，相应的Runge–Kutta方法就是半隐式方法，（5）式等号左右就含有级值k1，k2，⋯，kr，计算ki时要解一个包含r个未知量ki的非线性方程组。
如果A是满足Aij（i≤j）不全为零，则相应的方法就是隐式方法。
若系统是m维的，对于隐式Runge-Kutta方法实现的每一个递推步，均需要求解一个m×r维的非线性方程组，一般用牛顿迭代法求解。
例如
0120
0001200-120
162316
这是Kutta得到的3级3阶显式Runge–Kutta公式。
而
0121
00014140010
162316
与
0121
00052413-124162316
162316
分别是3级4阶的半隐式Runge–Kutta公式和隐式Runge–Kutta公式。
要求出具体的Runge–Kutta公式，一般有两种办法。
一种是将（6）在点xi,yi处进行泰勒展开，并且代入到（5）中，再与y（xi+h）在xi处的泰勒展开式进行对比，从而确定参数c，b，A。
另一种方法就是将微分方程转化为等价的积分方程，用数值积分得到表达式，再进行对比得到参数。
基于后一种方法，可以构造得到以下三种不同的隐式Runge-Kutta方法。
3.2Gauss型隐式RK方法[17]
设c1，c2，…，cs为Ps2c-1=0的根，这里Ps（t）是0,1上的s次Legendre多项式，0求s级的2s阶的Gauss型Runge–Kutta公式的参数c，b，A需要以下步骤:
（1）求出s次的Legendre多项式Ps2c-1=0的s个零点c1，c2，…，cs；
（2）由线性方程组
j=1sbjcjk-1=1kk=1,⋯,s
确定系数bj，j=1,⋯,s；
（3）计算系数矩阵
A=CW-1
在这个基础上，就给出了s=1,2,⋯,5，p=2s的一系列Gauss型隐式Runge–Kutta方法。
定理4[9][10]如果一个数值方法应用于方程y'=λy,其中λ为复常数，则对于所有λ，Reλ<0和对于固定步长h>0，当n→∞时，得到的数值解趋于零，则称这种方法是A稳定的。
定理5[9][10]s级Gauss型隐式Runge–Kutta方法是2s阶的，其稳定函数是s,sPade近似且是A稳定的。
以下列出s=3，p=6的Gauss型隐式Runge–Kutta方法之一：
12-15101212+1510
53629-1515536-1530536+152429536-1524536+153029+1515536
51849518
3.3Radau型隐式RK方法[17]
这种方法的参数c，b，A需要下面几个步骤求得：
（1）求多项式Pst-Ps-1（t）的零点c1，c2，…，cs，并指定c1>0，cs=1；
（2）计算系数bk，
bk=01Pst-Ps-1（t）t-ckPs'ck-Ps-1'（ck）dtk=1,⋯,s
（3）计算系数矩阵A，
A=CW-1
例如s=3，p=5的RadauⅡA型Runge–Kutta方法之一为：
4-6104+6101
88-76360296-16961800-2+36225296+1696180088+76360-2-3622516-63616+63619
16-63616+63619
定理6[31]s级RadauⅠA型Runge–Kutta方法和s级RadauⅡA型Runge–Kutta方法是2s-1阶的，稳定函数是（s-1,s）次对角Pade近似。
这两种都是A稳定的。
3.4Lobatto型隐式RK方法[17]
这种方法的参数c，b，A需要下面几个步骤求得：
（1）求多项式Pst-Ps-2（t）的零点c1，c2，…，cs，并指定c1=0，cs=1；
（2）计算系数bk，
bk=01Pst-Ps-2（t）t-ckPs'ck-Ps-2'（ck）dtk=1,⋯,s
（3）计算系数矩阵A，
A=D-1（WT）-1N-CTD
列出4级6阶RadauⅢA型隐式Runge–Kutta方法
05-5105+5101
011+5120025-512011-512011225+135120512025-1351200-1+512025+5120512-1-5120112
112512512112
定理7[31]s级LobattoⅢA，ⅢB，ⅢC型隐式Runge–Kutta方法是2s-2阶的，LobattoⅢA和ⅢB型Runge–Kutta方法的方法的稳定函数是（s-1,s-1）对角Pade近似，LobattoⅢC型隐式Runge–Kutta方法是（s-2,s）次对角Pade近似。
所以这些方法都是A稳定的。
4隐式RK方法的实现
4.1非线性系统的改进
对于s级隐式隐式Runge–Kutta方法
Yi=yn+hj=1saijfxn+cjh,Yji=1,⋯,s（7）yn+1=yn+hj=1sbjfxn+cjh,Yj（8）
令
zi=Yi-y（9）
于是（7）变为
zi=j=1saijfxn+cjh,yn+zj（10）
当得到的（10）解z1，⋯，z1时，由显式（8）可以很快得到解yn+1。
若隐式Runge–Kutta方法的系数矩阵式非奇异的，那么（10）可以写为
z1⋮zs=Ahf（xn+c1h,yn+z1）⋮hf（xn+csh,yn+zs）（11）
所以（8）可以看成与
yn+1=yn+i=1sdizi
等价的，其中
d1,⋯,ds=b1,…,bsA-1（12）
比如s=3，p=5的RadauⅡA型Runge–Kutta方法中，d=0,0,1。
4.2简化的Newton迭代法
对于一般的非线性微分方程，系统（10）必须用迭代的方法求解。
本文采用Newton迭代法。
Newton迭代法应用于系统（10），每次迭代都需要一个系数矩阵
I-ha11∂f∂yxn+c1h,yn+z1⋯ha1s∂f∂yxn+csh,yn+zs⋮-has1∂f∂yxn+c1h,yn+z1⋯hass∂f∂yxn+csh,yn+zs
的线性方程组。
定义5若s阶矩阵B=bij，m阶矩阵A=aij，且有
B⊗A=b11A⋯b1sA⋮⋱⋮bs1A⋯bssA
则称运算⊗是B与A的Kronecker积。
为了简化（10）的计算，我们用近似值
J≈∂f∂yxn,yn
代替Jacobi矩阵∂f∂yxn+cih,yn+zi的值，于是方程（10）的简化Newton迭代法变为
I-hA⊗J△Zk=-Zk+h（A⊗I）F（Z（k））
Z（k+1）=Z（k）+△Z（k）（13）
这里Z（k+1）=（z1k,⋯,zs（k））T是解的第k次近似，△Zk=△Z1k,⋯,△Zs（k）T是增量，F（Z（k））是FZk=fxn+c1h,yn+z1（k）,⋯,fxn+csh,yn+zs（k）T的缩写。
每一次的迭代需要s次由函数f的求值和一个N∙S维线性方程组的求解。
因为矩阵I-hA⊗J对于所有的迭代是相同的，所以其LU分解在一个积分步内只要进行一次，进而减少了计算时间。
5数值实验与结果分析
问题：
y1'=10-7×y1+sinxy10=1y2'=10-3×y2+cosxy20=1
写成矩阵形式就是：
y1'y2'=10-70010-3y1y2+sinxcosx
很明显该方程组的刚性比R=10-310-7=104≫1,因此是个刚性方程组。
取步长为0.1，在区间[0,1]内分别用四级四阶经典显式Runge–Kutta方法
yi+1=yi+16k1+2k2+2k3+k4k1=fxi,yik2=fxi+12h,yi+
升级会员 