专题九隐式龙格-库塔法.pdf
龙格-库塔(Runge-Kutta)法
龙格-库塔（Runge-Kutta）法 1.1 龙格-库塔(Runge-Kutta)法的基本思想
Euler公式可改写成
yi1 yi hK1 K1 f ( xi , yi )
则yi+1的表达式y(xi+1)与的Taylor展开式的前两项 完全相同,即局部截断误差为 O(h 2 ) 。
为了进一步提高精度，设除 xi p 外再增加一点
xiq xi qh ( p q 1)
并用三个点 xi ,xi p , xiq 的斜率k1,k2,k3加权平均
得出平均斜率k*的近似值，这时计算格式具有形式:
yi1 yi h(1 )k1 k2 k3
k1 f (xi , yi ) k2 f (xi ph, yi phk1 )
格式。
若取 1 0 ,则 2 法的计算公式为
1,
p
1 2
，此时二阶龙格-库塔
ky1i
1
f
yi hk2 ( xi , yi )
k
2
h
f
(
x
i
1
,
yi
2
2 k1 )
i 0,1,2, n 1
此计算公式称为变形的二阶龙格—库塔法。式中
x 1 i 2
为区间
xi , xi1
的中点。
1.3 三阶龙格-库塔法
拉法,将 xi p 视为 xi1，即可得
k2 f (xi ph, yi phk1 ) 对常微分方程初值问题(7.1)式的解 y=y(x),根据微 分中值定理，存在点 (xi , xi1 ) ，使得
也即
y(xi1 ) y(xi ) y( )( xi1 xi )
y( xi1 ) y( xi ) hK
龙格库塔法推导
于是可考虑用函数f(x,y)在若干点上的函数值的 线性组合来构造近似公式，构造时要求近似公式在 (xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式 的前面几项重合，从而使近似公式达到所需要的阶数。 既避免求高阶导数,又提高了计算方法精度的阶数。 或者说,在[xi,xi+1]这一步内多计算几个点的斜率值， 然后将其进行加权平均作为平均斜率，则可构造出更 高精度的计算格式，这就是龙格—库塔（Runge龙格— 龙格 库塔（ Kutta）法的基本思想 ）法的基本思想。
[
]
K1 =hf(xi, yi)
∴
yi +1 = yi + (c1K1 + c2 K 2 )
= y ( xi ) + c1hf ( xi , yi )
′ ′ + c2 h f ( xi , yi ) + a2 hf x + b21hff y + O ( h 3 )
[
]
= y ( xi ) + (c1 + c2 ) hf ( xi , yi )
存在无穷多个解。所有满足上式的格式统称为2阶 无穷多个解。 无穷多个解 阶 龙格 - 库塔格式。 库塔格式。
2011-12-8 14
1 注意到， 就是二阶龙格 注意到，a 2 = b21 = 1, c1 = c2 = 就是二阶龙格 - 库塔 2 y 公式，也就是改进的欧拉法 改进的欧拉法。 公式，也就是改进的欧拉法。 K2
称为P阶龙格－库塔方法。 其中ai,bij,ci为待定参数，要求上式yi+1在点(xi,yi)处作 Tailor展开，通过相同项的系数确定参数。
2011-12-8 7
龙格-库塔方法-文档资料
c3a 2b32
c3a3 1
6
; 2
O (h4)
常见的2种三阶方法：
库塔三阶方法
h yn1yn6(k14k2k3)
k1
f(xn,yn);
k2
hh f(xn2,yn2k1)
k 3 f(x n h ,y n h k 1 2 h k 2 ) •5
四级方法：N = 4
局部截断误差 O ( h 5 )
可见误差随着 x n 的增加呈指数函数递减
当 f y 0 时，微分方程是不稳定的； 而 f y 0 时，微分方程是稳定的。
上面讨论的稳定性，与数值方法和方程中 f 有关
•21
实验方程： y y C ,R e () 0
D e f 3 对单步法 yn 1ynh(xn,yn,h )应用实验方程，
e n 1 e n h [ ( x n ,y ( x n ) , h ) ( x n ,y n , h ) ] T n 1
•15
因为单步法是 p 阶的：h0,0hh0满足|Tn1|Chp1
|e n 1| |e n| h L |e n| C h p 1|en |
其中 1hL,C hp1
•18
三、绝对稳定性 /*Absolute Stibility*/ 计算过程中产生的舍入误差对计算结果的影响
首先以Euler公式为例，来讨论一下舍入误差的传播:
yn1ynhf(xn,yn)
设实际计算得到的点 x n 的近似函数值为 yn yn n，
其中 y n 为精确值， n 为误差
yn1ynhf(xn,yn)
通过适当选取参数 1,2和 p 的值，使得公式具有 2阶精度！！
•3
由泰勒公式展开，要使公式具有 2 阶精度，只需
龙格-库塔方法(Runge-Kutta)
龙格-库塔方法〔Runge-Kutta〕3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的一般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值方法，若用显式单步法(3.2.2)当，即数值求积用左矩形公式，它就是Euler法(3.1.2)，方法只有一阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数f的值，但方法是二阶精度的.若要得到更高阶的公式，则求积分时必须用更多的f值，根据数值积分公式，可将(3.2.1)右端积分表示为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)一样，用前面已算得的f值表示为(3.2.3)，一般情况可将(3.2.2)的 表示为(3.2.4)其中这里均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K方法.它每步计算r个f值(即由前面(i-1)个已算出的表示，故公式是显式的.例)，而ki如当r=2时，公式可表示为(3.2.5)其中.改进Euler 法(3.1.11)就是一个二级显式R-K 方法.参数取不同的值，可得到不同公式.3.2.2 二、三级显式R-K 方法对r=2的显式R-K 方法(3.2.5)，要求选择参数，使公式的精度阶p 尽量高，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6)令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯一.因r=2，由〔3.2.5〕式可知，于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多，只要都可得到二阶精度R-K方法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)与中点公式(3.2.9)是两个常用的二级R-K方法，注意二级R-K方法只能达到二阶，而不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个方程,加上(3.2.7)的三个方程，共计六个方程求4个待定参数，验证得出是无解的.当然r=2,p=2的R-K方法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对r=3的情形，要计算三个k值，即其中将按二元函数在处按Taylor公式展开，然后代入局部截断误差表达式，可得可得三阶方法，其系数共有8个，所应满足的方程为(3.2.10)这是8个未知数6个方程的方程组，解也是不唯一的，通常.一种常见的三级三阶R-K方法是下面的三级Kutta方法：(3.2.11)附：R-K 的三级Kutta 方法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满足c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K 方法与步长的自动选择利用二元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K 方法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，而33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。
龙格-库塔(Runge-Kutta)方法
——常微分方程及方程组的数值解 常微分方程及方程组的数值解
基本概念 只有一个自变量的微分方程为常微分方程 ODE, 否则称为偏微分方程 PDE。 。 方程中未知函数导数的最高阶数称为方程的阶
[1] [2] C. Runge, "Ueber die numerische Auflsung von Differentialgleichungen" Math. Ann. , 46 (1895) pp. 167–178 W. Kutta, "Beitrag zur naherungsweisen Integration von Differentialgleichungen" Z. Math. und Phys. , 46 (1901) pp. 435–453
2. 递推化 在具有唯一解的条件下， 在具有唯一解的条件下，通过步进法逐步计算出解 在一系列离散点上的值。从而得到原 的数值近似解。 在一系列离散点上的值。从而得到原ODE的数值近似解。 的数值近似解
初值问题解法 讨论一阶ODE（高阶可化为一阶ODEs） （高阶可化为一阶 讨论一阶 ） 的初值问题。初值问题可以一般地写成： 的初值问题。初值问题可以一般地写成：
欧拉（ 欧拉（Euler）方法 ）
Euler方法是求解初值问题的最简单方法， 方法是求解初值问题的最简单方法， 方法是求解初值问题的最简单方法 精度差。然而对理论分析很有用。 精度差。然而对理论分析很有用。Runge-Kutta 法是对Euler法的改进 法是对 法的改进
Euler方法： 方法： ym+1 = ym + hK1 + O(h ) K1 = f ( xm , ym )
K 1 = f ( xm , ym ) K i = f ( xm + hai , ym + h∑ bil K l )
Runge-Kutta法
3½×RK· ¨
1.422
¸ ĽøEuler· ¨
1.42
1.418 0.2965 0.297 0.2975 0.298 0.2985 0.299 0.2995 0.3 0.3005 17
§ 7.2 Runge-Kutta法
龙格-库塔(Runge-Kutta)方法简称R-K法，是一种应用较广的 高精度的单步法。
所谓单步法就是在计算yi时只用到前一步信息yi-1的方法。
本节介绍R-K法的构造原理、常用公式。
1
一、Runge-Kutta方法的构造原理 对于常微分方程的初值问题
y f ( x , y ) y( a ) y0 a xb
另一方面
h2 y( xi 1 ) y( xi ) hy( xi ) y( xi ) o(h3 ) 2
9
式（6）的局部截断误差:
en (h) y ( xi 1 ) yi 1 1 h 1 (1 2 ) y( xi ) h2 22 y( xi ) o(h3 ) 2
1 2 1 1 2 2 2
有无穷多组解，从而可以得到许多具体的二阶R -K公式 如：
1 2 1，2 1
1 1 0， 2 1， 2 2
10
用类似的方法也可以构造三阶R-K公式
h yi yi 1 ( K1 4 K2 K3 ) 6 K1 f ( xi 1, yi 1 ) h h K 2 f ( xi 1 , yi 1 K1 ) 2 2 K3 f ( xi 1 h, yi 1 h(2 K2 K1 ))
因而方法(8)有4阶精度
12
例1. 使用高阶R-K方法计算初值问题
第17节 龙格-库塔方法
10
15
20
Van der Pol 方程
1000
2 1.5 1 0.5 0 -0.5 -1
1500 1000 500 0 -500
-1.5 -2 -2.5 0 500 1000 1500 2000 2500 3000
-1000 -1500 0
500
1000
1500
2000
2500
3000
刚性问题
2 3 h h 2 2 2 y( n ) h( y( n ) tn 1) ( y( n ) tn 1 2tn ) ( y( n ) tn 1 2tn ) 2 6
龙格-库塔方法
目标：求解一阶方程组：y ' f ( x, y ), y ( x(0) ) 0 根据牛顿 - 莱布尼兹公式，有： y ( x( n 1) ) - y ( x( n ) )
高阶微分方程组求解
将高阶微分方程组转换为多元微分方程组。
对于二阶微分方程：x '' f (t , x ', x) 可通过引入中间变量 : y x '， 将其转化为一阶微分方程组进行求解： y ' f (t , y, x) x ' y
Van der Pol 方程
电 容 C 电感 L
微分方程组的求解
x ' f (t , x, y) 二元微分方程组： y ' g (t , x, y) x(t0 ) x0 y(t0 ) y0
n元微分方程组：x ' f (t , x), x(t0 ) x0 , 其中, x, f为n维向量、函数。
n元线性微分方程组可表示为矩阵形式： x ' Ax, x(t0 ) x 0 , 对于此类问题，刚性问题可采用如下方式判定： max(re(i ) S min( re(i ) S 较大，则问题为刚性的。
刚性微分方程组隐式龙格库塔方法
毕业设计题目：刚性系统的隐式RK方法学院：数理学院专业名称：信息与计算科学学号： ************学生姓名：**指导教师：***2016年05月15日摘要本文主要介绍单步隐式Runge – Kutta方法,简要的介绍了Gauss型隐式Runge – Kutta方法、Radau型隐式Runge – Kutta方法和Lobatto型隐式Runge – Kutta方法。
并利用这些基本的隐式Runge – Kutta方法来对刚性方程组进行数值求解，并将隐式Runge – Kutta方法与显式经典Runge – Kutta方法求解的结果进行对比，说明两种数值解法的优缺点。
关键词：刚性系统隐式Runge – Kutta方法单步方法Newton迭代法AbstractThis paper mainly introduces the Implicit Runge-Kutta Methods and a simple description of Gauss implicit Runge-Kutta method ,Radau implicit Runge-Kutta method and Lobatto implicit Runge-Kutta method . These basic Implicit Runge-Kutta methods are used to solve the stiff equations. These implicit Runge-Kutta methods iare compared with the classical explicit Runge-Kutta method. This paper explain the advantages and disadvantages of the two kind of numerical methods.Keywords: Stiff system Implicit Runge-Kutta method One step method Newton iterative method目录1.绪论 (1)1.1刚性方程 (1)1.2隐式RK方法的研究意义 (2)1.3 RK方法的研究现状 (3)2.单步RK方法的收敛性和稳定性 (5)2.1单步RK方法的收敛性 (5)2.2单步RK方法的稳定性 (6)3.三类隐式RK方法 (8)3.1引言 (8)3.2 Gauss型隐式RK方法 (9)3.3 Radau型隐式RK方法 (10)3.4 Lobatto型隐式RK方法 (11)4隐式RK方法的实现 (13)4.1非线性系统的改进 (13)4.2简化的Newton迭代法 (13)5数值实验与结果分析 (15)参考文献 (18)附录 (21)1.绪论1.1刚性方程对于一般的线性常系数系统y′=Ay+φ(t) A为m×m的矩阵，特征值为λi(i=1,2,⋯,m)。
龙格库塔方法ppt课件
z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(
i1,6)=y;
i1=i1+1;
end
11
end
例题4
例4 求(3阶、4阶R ห้องสมุดไป่ตู้公式)解初值问题(步长h 0.1)
y
y
2x y
y(0) 1
(0 x 1)
将步长折半，即取
h 2
为步长从xn跨两步到xn+1，求得一个近似值yn( h21)
每跨一步截断误差是C
(
h 2
)5
,
有y(xn+1
)-y(
h 2
)
n1
2C(
h )5 2
16
步长折半后，误差大约减少为原来的 1 ，即有 16
(h)
y( xn1) y( xn1)
y2 n1
y(h) n1
误差 0.45e-4 0.17e-4 0.15e-4 0.48e-4 0.25e-4 0.55e-4 0.14e-4 0.21e-4 0.54e-4 0.06e-4
y(xn) 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125 1.6733 1.7321
y=y0;z=zeros(c,6);
%z生成c行，6列的零矩阵存放结果；
%每行存放c次迭代结果，每列分别存放
k1~k4及y的结果
10
不断迭代运算： for x=a:h:b
if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4);
(完整版)第二节龙格-库塔方法
因为单步法是 p阶的：h0 , 0 h h0 满足| Tn1 | Chp1
| en1 || en | hL | en | Ch p1 | en |
其中 1 hL, Ch p1
en O(hp )
即取
K * 1K1 2 K2 yi1 yi h(1K1 2 K2 )
其中1 和 2 是待定常数。若取 K1 f ( xi , yi ) ，则
问题在于如何确定 xi p 处的斜率 K2 和常数 1 和 2 。
仿照改进的欧拉方法，用欧拉方法预测 y( xi p ) 的值，
yi p yi phK1
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2：用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x (0,1)
解：经典的四阶龙格-库塔公式： y(0) 1
E(h) 1 h
绝对稳定域： 1 h 1
当 R 时， 1 h 1 2 h 0
绝对稳定区间：(2, 0)
❖经典的R-K公式：yn1
yn
h 6
(k1
2k2
2k3
k4 )
k1 f ( xn , yn ) yn
k2
f
(
yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
第3章03龙格-库塔方法
这样有计算公式：
yn1 yn h[(1 )k1 k2 )
k1 f (xn , yn ), k2 f (xn ph, yn phk1)
y(h) n1
1 16
y( xn1)
(h)
y2 n1
1 15
(h)
y2 n1
y(h) n1
误差事后估计
可以通过检查步长折半前后两次计算结果的偏差
y y ( h ) 2
(h)
n1
n1
来判断所选取的步长是否合适.
1) 对于给定的精度ε,如果δ> ε,则反复将步长折半 进行计算,直到δ< ε为止,这时取步长折半后的 新值作为结果;
k4
k3 f ( xn 2h, yn 21hk1 22hk2 ) f ( xn 3h, yn 31hk1 32hk2 33hk3)
四阶龙格－库塔公式每步要四次计算函数值，具有四阶精 度，局部截断误差是O(h5) .
四阶龙格－库塔格式
下列经典格式是其中常用的一种:
yn1 k2
yk
1.000000 1.005000 1.019025 1.041218 1.070800 1.107076 1.149404 1.197210 1.249975 1.307228 1.368514
y( xk ) yk
0.0 1.6×10－4 2.9×10－4 4.0×10－4 4.8×10－4 5.5×10－4 5.9×10－4 6.2×10－4 6.5×10－4 6.6×10－4 6.6×10－4
调用 2 次函数。有 2 阶精度。
龙格库塔法
c1
c2
1
0,
1 2
c2
0
即常数c1, c2 , 满足条件
c1 c2 1
c2
1 2
方程组有三个未知数，但只有两个方程，因此可得到
局部截断误差为O(h3 )的计算公式.
如果取c1
c2
1 ,
2
1,递推公式为
y0
k1 f (ti , yi )
k2 f (ti h, yi hk1)
yi
代入上式, 得
yi1 yi h(c1 f c2 f ) c2h( ft ffy ) O(h2 )
在局部截断误差的前提假设yi y(ti )下,得
y(ti1)
yi1
h(c1
c2
1)
f
h2(1 2
c2 )( ft
ffy ) O(h3)
要使局部截断误差y(ti1) yi1 O(h3 ),当且仅当
§9-3 龙格—库塔法
一、高阶泰勒法
假设初值问题
dy f (t, y) a t b dt
(1)
y(a)
的解y(t)及f (t, y)足够光滑.
将y(ti1)在ti处作n阶泰勒展开 , 得
y(ti1)
y(ti )
y(ti )h
y(ti ) h2 2!
y(n) (ti ) hn n!
y(n1) (i ) hn1
f
(ti
1 2
h,
yi
1 2 hk1)
(10)
yi1 yi hk2 )
公式(8)、(9)、(10)三式是三种常见的二阶龙格—库塔公式
局部截断误差为 O(h3).
三、三、四阶龙格—库塔法
三阶龙格—库塔法
数值分析9-3(龙格-库塔方法)
总结词
除了Python和MATLAB，还有许多其他编 程语言可以用于实现龙格-库塔方法。
详细描述
例如C、Java和R等编程语言也提供了相应 的数值计算库或框架，可以实现龙格-库塔 方法。使用这些语言实现龙格-库塔方法需 要一定的编程基础和对相应语言的数值计算 库的了解。
龙格-库塔方法可以用于求解偏微分方程的数值解，通过将偏微分方程转化为常微分方程组，利用龙格 -库塔方法进行迭代求解，能够得到较为精确的结果。
积分方程的数值解
积分方程是描述函数与积分之间的关 系的数学模型，常见于物理、工程等 领域。
VS
龙格-库塔方法也可以用于求解积分 方程的数值解，通过将积分方程转化 为常微分方程组，利用龙格-库塔方 法进行迭代求解，能够得到较为精确 的结果。
重要性及应用领域
龙格-库塔方法是数值分析中非常重要的内容， 它为解决常微分方程提供了一种有效的数值方 法。
在科学、工程和经济学等领域中，许多问题都 可以转化为求解常微分方程的问题，因此龙格库塔方法具有广泛的应用价值。
例如，在物理学、化学、生物学、金融学等领 域中，龙格-库塔方法被广泛应用于模拟和预测 各种动态系统的行为。
数值分析9-3：龙格-库塔方法
目录
• 引言 • 龙格-库塔方法概述 • 龙格-库塔方法在数值分析中的应用 • 龙格-库塔方法的实现与编程 • 龙格-库塔方法的改进与优化 • 结论与展望
01 引言
主题简介
龙格-库塔方法是一种用于求解常微 分方程的数值方法。
它通过构造一个离散化的时间序列来 逼近微分方程的解，并利用已知的离 散点来计算新的离散点，逐步逼近微 分方程的真实解。
02 龙格-库塔方法概述
定义与原理
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
∫ ∑ =ຫໍສະໝຸດ h⎡ ⎢⎣1 0
f
(tn
+ hu) du −
s
bi
i =1
f
(tn
+ cih)⎤⎥⎦
∫ ∑ 1 0
f
(tn
+ hu) du
≈
s i =1
bi
f
(tn
+ cih)
(9.8)
具有 q −1次代数精确度，即(9.8)对 f (tn + hu) = (hu)k−1 ， k = 1, 2,L, q 精确成立，
如果当 j ≥ i 时 aij = 0 ，公式(4.2b)中 ki 可由 k1,L, ki−1 得到，此时公式(9.2a)和
(9.2b)为显式方法，若当 j > i 时，aij = 0 而 aii = 0 ，公式称为半隐式 RK 方法，此
时(9.2b)右端 f 含 ki ，因此要求得 ki ，需要解 m 阶非线性方程组（假定 f 关于 y 为
点
c1
,L
cs
就是拉道求积节点，它分别为
Ps*
(
t
)
+
P* s −1
(t
)
的零点或
Ps*
(
t
)
−
Ps
* −1
(
t
)
的
零点．这类方法计算节点及系数 c,b, A 的步骤是：
(1)求
Ps*
(t
)
+
P* s −1
(t
)
=
0
或
Ps*
(t
)
−
P* s −1
(t
)
=
0
的实根，得到
c1 ,L cs
，前者含
c1 = 0 ，后者含 cs = 1．
由以上 3 步求得的 s 级 2s 阶隐式 RK 方法如下： s = 1, p = 2 的隐式 RK 方法公式
11 22
1
yn+1
=
yn
+
hf
⎛ ⎜⎝
tn
+
1 2
h,
1 2
(
yn
+
yn+1
)
⎞ ⎟⎠
s = 2, p = 4 的隐式 RK 方法是
（9.13）
1− 3 26 1+ 3 26
1 4 1+ 3 46
对(9.1)中方程从 tn 到 tn+1 积分可得
∫ ( ) y (tn+1 ) − y (tn ) =
f tn+1
tn
t, y (t ) dt
若用 s 个节点的求积公式计算右端积分，即
∫ ( ) ∑ ( ) ( ) ( ) f tn+1 t, y t tn
s
dt ≈ h bi f tn + cih, y tn + cih
i =1
j =1
(9.12)
这就是隐式 RK 方法的相容性条件． 根据以上讨论可利用高斯型求积公式建立隐式 RK 方法，下面给出 3 种隐式
RK 方法． （I）基于高斯求积公式的隐式 RK 方法
若用 s 点的高斯求积公式．其代数精确度为 2s −1，即 q −1 = 2s −1，故 q = 2s ，
(9.10)
此公式至少具有 s −1次代数精确度，即对 f (tn + hu ) = (hu)k−1 , k = 1, 2,L, s 精确成
立，由此得到
∑s
j =1
aijc j k −1
=
1 k
cik , i
= 1, 2,L, s
(9.11)
对 k = 1, 2,L, s 成立，它表明当求积节点 ci (i = 1, 2,L, s) 确定后，由(9.11)可求出系
由此得到
∑s
i =1
bicik −1
=
1 k
,k
= 1, 2,Lq
(9.9)
当 f (tn + hu) = (hu)q 时，求积公式(9.8)不精确成立，此时
( ) ( ) y tn+1 − yn+1 = O hq+1
若求积节点 ci (i = 1, 2,L, s) 已经给定，对(4.7)从 tn 到 tn,i = tn + cih 求积，则得
于是有
∫ ∫ y (tn+1 ) =
tn+1 tn
f
(t ) dt
=
h
1 0
f
(tn
+ hu) du
s
∑ ≈ h bi f (tn + cih) = yn+1 i=0
若求积公式
( ) ∫ ( ) ∑ ( ) y tn+1
− yn+1 =
f tn+1
tn
t
s
dt − h bi f
i =1
tn + cih
( ) Ps ( x) =
1 ds 2s s! dxs
⎡ ⎢⎣
x2 −1
s⎤ ⎥⎦
令 x = 2t −1，则得 Ps* (t ) = Ps (2t −1) ，当 s = 1, 2,3 时有
P1* (t ) = 2t −1 P2* (t ) = 6t2 − 6t +1 P3* (t ) = 20t3 − 30t2 +12t −1
, j = 1, 2,Ls
对 l = 1, 2,L, s 成立．
（9.16）
这样得到的公式很多，下面只给出 s = 2, p = 3 的两个公式．当 s = 2 时，由
P2* (t ) +
P1*
(t)
=
6t 2
−
4t
=
0
求得
c1
=
0, c2
=
2 3
，再由（9.9）求得
b1
=
1 4
, b2
=
3 4
，
最后由（9.11）求得 a11
1/ 2
1− 3 46
1 4
1/ 2
s = 3, p = 6 的隐式 RK 方法是
（9.14）
5 − 15 10
1 2
5 + 15 10
5 36 1+ 3 46
1/ 2
10 − 3 15 45
2 9
10 + 3 15 45
25 − 6 15 180
10 − 3 15 72
5 / 36
5 /18
1/18
∫ ( ) y (tn + cih) − y (tn ) =
f tn +cih
tn
t, y (t ) dt,i = 1,Ls
若右端积分用求积系数为 aij (i, j = 1, 2L s) 的求积公式，则得
∑ ( ) s
Yi = yn + h aij f tn + cih,Yj , i = 1,L s
y (tn+1 ) 之差 ln+1 ，当 h → 0 时，如果有整数 p > 0 ，使
( ) ( ) ln+1 = y tn+1 − yn+1 = O h p+1
(9.5)
则称方法(9.2)是 p 阶相容的， ln+1 为公式(9.2)的局部截断误差，对 s 级的显式 RK
方法的阶 p ≤ s ，当 s ≤ 4 时有 p = s ．而对隐式 RK 方法的阶 p 最高可达到
=
a12
=
0
及 a21
=
a22
=
1 3
，于是得公式
0
0
0
211 333
1/ 4 3/ 4
(9.17a)
这公式称为 Radau IA 公式，它不是 A − 稳定的，一般不用．另一个 A − 稳定的
Radau IA 公式，其 aij 由（9.16）求出的，公式为
0
1 4
−
1 4
2 15 3 4 12
1/ 4 3/ 4
( ) 这样得到的隐式 RK 方法局部截断误差 ln+1 = ( ) y tn+1 − yn+1 = O h2s+1 ，故方法
的阶 p = 2s ，对 s 级 2s 阶的隐式 RK 方法(4.2)其节点及系数 c,b, A 可由以下 3 步 得到．
(1)求节点 c1, c2,L, cs ，由于[0,1] 上的高斯求积公式节点为区间[0,1] 的 s 阶勒 让德多项式 Ps* (t ) 的 s 个零点，而[−1,1] 上的 s 阶勒让德多项式为
当区间为[0,1] ，其两端均为求积节点，对 s ≥ 2 的洛巴托求积公式，代数精
j =1
为了建立求积系数 bi 与节点系数 ci 之间的关系，可通过对特殊初值问题积分得
到，考虑
y ' = f (t), y (tn ) = 0
(9.7)
由 tn 到 tn+1 积分得
( ) ∫ ( ) y tn+1
=
f tn+1
tn
t dt
令 t = tn + hu ，并用 s 个节点的求积公式得到
Yi = yn + h aij f tn + c jh,Yj ， i = 1, 2,Ls j =1
(9.4b)
计算时，先由(9.4b)求出Y1,LYs ，再由(9.4a)求出 yn+1 ．若假定在 tn 点(9.1)的解 y (tn )
与数值解 yn 相等，即 yn = y (tn ) ．由(9.4)求出 tn+1 = tn + h 的数值解 yn+1 与精确解