081数值计算方法—常微分方程(组)
常微分方程数值解法
用分段的折线逼近函数,此为 “折线法”而非“切线法”， 除第一个点是曲线上的切线，
其它都不是。
2、Euler方法的误差估计
1)局部截断误差。 在一步中产生的误差而非累积误差：
~
T x y y
n1
n1
n1
其中
~
y
是当
y
n
y(
x
)
n
（精确解！）时
n1
由Euler法求出的值，即y 无误差！ n
y y x , y x y 则得： h f
f ,
n1
n2
n
n
n1
n1
同样与Euler法结合，形成迭代算法，对n 0，1，2，
y y x , y 0 hf
n1
n
n
n
y y x , y x y
k 1
推出总体误差与步长的关系。
由微分方程解的存在唯一性,自然假定 （ f x，y）
充分光滑，或满足 Lipschitz条件：
f
x,ny源自xn
f
x
,
n
y
n
L
yxn
y n
第 n 步 的 总 体 截 断 误 差 记 为
en y
xn
y n
则 对 n 1 步：
e x y x y y y T y y ~ ~
用yn1, yn代替y(xn1), y(xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f
(x,
y)dx
hf
(xn ,
yn )
常微分方程组的解法
常微分方程组的解法常微分方程组是由多个关于未知函数及其导数的方程组成的方程组，它是数学中的重要研究对象。
常微分方程组的解法可以分为解析解法和数值解法两种。
解析解法是指通过数学方法求出常微分方程组的解析表达式。
常微分方程组的解析解法主要包括分离变量法、一阶线性方程法、变量代换法、常数变易法、特殊函数法等。
其中，分离变量法是指将常微分方程组中的各个变量分离出来，然后对每个变量分别积分，最后得到常微分方程组的解析解。
一阶线性方程法是指将常微分方程组转化为一阶线性方程，然后通过求解一阶线性方程来得到常微分方程组的解析解。
变量代换法是指通过合适的变量代换将常微分方程组转化为更简单的形式，然后通过求解简化后的方程组得到常微分方程组的解析解。
常数变易法是指将常微分方程组中的常数作为未知量，然后通过求解常数得到常微分方程组的解析解。
特殊函数法是指通过特殊函数的性质求解常微分方程组，如指数函数、三角函数等。
数值解法是指通过计算机数值计算的方法求出常微分方程组的数值解。
常微分方程组的数值解法主要包括欧拉法、龙格-库塔法、变步长法等。
其中，欧拉法是一种简单的数值解法，它的基本思想是将常微分方程组的解曲线上的点离散化为一系列点，然后通过计算机逐步求解得到常微分方程组的数值解。
龙格-库塔法是一种高阶数值解法，它通过计算机采用多个不同的计算公式来逼近常微分方程组的解曲线，从而得到更为准确的数值解。
变步长法是一种自适应数值解法，它通过计算机根据误差大小自动调整步长大小，从而得到更为准确的数值解。
常微分方程组的解法包括解析解法和数值解法两种，每种方法都有其适用的范围和优缺点。
在实际应用中，需要根据具体问题的性质和求解要求选择合适的解法来求解常微分方程组。
常微分方程的数值解法
常微分方程的数值解法…………江南大学信计1203柯恒一、前言对于很多微分方程，我们很难求出解析解，这时我们需要采取数值手段求解。
在数值分析这门课中，老师讲到dy dt =f (t ,y )的数值解法，我们采用了欧拉格式，后退欧拉格式，梯形公式，改进欧拉格式及4阶龙格-库塔方法求解。
而老师没讲关于微分方程组和高阶微分方程的解法。
现在我来粗虐的说一下微分方程组及高阶微分方程的数值解，这里主要是借助matlab 中的相关函数进行计算并仿真出图像。
二、相关问题初值问题问题1，微分方程组问题描述：给定一个微分方程，并且告诉我们初始状态。
我们便可求出整个过程的解。
给定下面方程（L 表示省略号）：11122112112(,,,,)(,,,,)(,,,,)n n n n dy f x y y y dx dy f x y y y dx dy f x y y y dx ⎧=⎪⎪⎪=⎪⎨⎪⎪⎪=⎪⎩而且初始状态y ′i (x 0)已知记成y(0).问题解决：四阶龙格库塔方法：K i1=f i (x n ,y 1n ,y 2n ,……,y nn );K i2=f i (x n +h 2,y 1n +h 2∗k 11,y 2n +h 2∗K 21,……,y nn +h 2∗K n1) K i3=f i (x n +h ,y 1n +h ∗K 12,y 2n +h ∗K 22,……,y nn +h ∗K n2) K i4=f i (x n +h 2,y 1n +h 2∗K 13y 2n +h 2∗K 32,……,y nn +h 2∗K n3),y i,n+1=y i,n+h6∗(K i1+2∗K i2+2∗k i3+2∗K i4)通过龙格库塔方法，我们可以计算后面点的值，在matlab中调用ode45来实现四阶龙格库塔方法的调用。
应用举例：有一个同步地球卫星，现在加速到4km/s进行变轨，试分析变轨后卫星运动的轨迹。
已知地球的质量M=5.97e24,引力常量的值为6.672e-11问题建模：我们已经知道行星的运动是在一个平面上的。
常微分方程中的数值方法
常微分方程中的数值方法常微分方程是数学中的一个重要分支。
它主要研究的对象是随时间变化的函数。
在实际应用中，我们需要求解这些函数的解析解，但通常情况下，解析解并不容易得到，甚至是不可能得到。
因此，我们需要使用数值方法来求解这些函数的数值近似解。
在本文中，我们将介绍常微分方程中的数值方法。
一、欧拉法欧拉法是常微分方程数值解法中最基本的一种方法。
它是根据欧拉公式推导而来的。
具体地，我们可以将一阶常微分方程dy/dt=f(t,y)写成如下形式：y(t+h)=y(t)+hf(t,y(t))其中，h是步长，f(t,y)是t时刻y的导数。
欧拉法就是通过上面的公式进行逐步逼近，然后得到最终的数值解。
欧拉法的计算过程非常简单，但所得到的解可能会出现误差。
这是因为欧拉法忽略了f(t+h,y(t+h))和f(t,y(t))之间的变化。
因此，我们需要使用更为精确的数值方法来解决这个问题。
二、改进欧拉法为了解决欧拉法中的误差问题，我们可以使用改进欧拉法。
改进欧拉法又称作四阶龙格-库塔法。
它的基本思想是对欧拉法公式进行改进，以提高计算精度。
具体地，根据龙格-库塔公式，可将改进欧拉法表示为：y(t+h)=y(t)+1/6(k1+2k2+2k3+k4)其中，k1=h*f(t,y)k2=h*f(t+h/2,y+k1/2)k3=h*f(t+h/2,y+k2/2)k4=h*f(t+h,y+k3)改进欧拉法的计算过程比欧拉法要复杂些，但所得到的数值解比欧拉法更精确。
这种方法适用于一些特殊的问题，但在求解一些更为复杂的问题时，还需要使用其他的数值方法。
三、龙格-库塔法龙格-库塔法是求解常微分方程中数值解的常用方法之一。
它最常用的是四阶龙格-库塔法。
这种方法的基本思想是使用四个不同的斜率来计算数值解。
具体地，我们可以将四阶龙格-库塔法表示为：y(t+h)=y(t)+1/6(k1+2k2+2k3+k4)其中，k1=h*f(t,y)k2=h*f(t+h/2,y+k1/2)k3=h*f(t+h/2,y+k2/2)k4=h*f(t+h,y+k3)与改进欧拉法相比，龙格-库塔法的计算复杂度更高，但所得到的数值解更为精确。
微分方程组的数值求解方法
微分方程组的数值求解方法微分方程组数值求解方法微分方程组是数学中非常重要的一个分支，它描述了许多自然界和社会生活中的现象，例如电路的运行、天体的运行、生命体的生长等等。
我们需要对微分方程组进行求解，才能够得到它们的解析解，从而更好地理解和应用它们。
然而，大多数微分方程组不可能用解析法求解，因此，我们需要采用数值方法来求解微分方程组。
常见的微分方程组数值求解方法包括欧拉法、龙格库塔法和变步长法等。
下面，我们将逐一介绍它们的基本原理和优缺点。
一、欧拉法欧拉法是微分方程组数值求解方法中最简单的一种。
它的基本思想是将微分方程组中的各个变量离散化，然后根据微分方程组的导数计算每一步的值。
具体来讲，欧拉法的数值求解公式为：\begin{aligned} &x_{n+1}=x_n+hf_n(x_n,y_n,z_n),\\&y_{n+1}=y_n+hf_n(x_n,y_n,z_n),\\&z_{n+1}=z_n+hf_n(x_n,y_n,z_n), \end{aligned}其中，$x(t)$，$y(t)$，$z(t)$是微分方程组的解，$f_n(x_n,y_n,z_n)$是微分方程组导数在点$(x_n,y_n,z_n)$处的值，$h$为时间步长。
欧拉法的优点是简单易懂，方便实现，缺点是误差较大，计算不够精确。
因此，在实际应用中，往往需要采用更加精确的数值方法。
二、龙格库塔法龙格库塔法是微分方程组数值求解方法中比较常用的一种。
它的基本思想是通过多次计算微分方程组中的导数，以获得更加精确的数值解。
具体来讲，龙格库塔法的求解公式为：\begin{aligned}&k_{1x}=hf_n(x_n,y_n,z_n),k_{1y}=hf_n(x_n,y_n,z_n),k_{1z}=hf_n (x_n,y_n,z_n),\\&k_{2x}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{1y}}{2},z_n+\frac{k_ {1z}}{2}),k_{2y}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{1y}}{2},z_n+ \frac{k_{1z}}{2}),k_{2z}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{1y}}{ 2},z_n+\frac{k_{1z}}{2}),\\&k_{3x}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{2y}}{2},z_n+\frac{k_ {2z}}{2}),k_{3y}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{2y}}{2},z_n+ \frac{k_{2z}}{2}),k_{3z}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{2y}}{ 2},z_n+\frac{k_{2z}}{2}),\\&k_{4x}=hf_n(x_n+h,y_n+k_{3y},z_n+k_{3z}),k_{4y}=hf_n(x_n+h,y_n+k_{3y},z_n+k_{3z}),k_{4z}=hf_n(x_n+h,y_n+k_{3y},z_n+k_{3 z}),\\&x_{n+1}=x_n+\frac{k_{1x}}{6}+\frac{k_{2x}}{3}+\frac{k_{3x}}{ 3}+\frac{k_{4x}}{6},\\&y_{n+1}=y_n+\frac{k_{1y}}{6}+\frac{k_{2y}}{3}+\frac{k_{3y}}{ 3}+\frac{k_{4y}}{6},\\&z_{n+1}=z_n+\frac{k_{1z}}{6}+\frac{k_{2z}}{3}+\frac{k_{3z}}{ 3}+\frac{k_{4z}}{6}, \end{aligned}其中，$k_{1x}$，$k_{1y}$，$k_{1z}$，$k_{2x}$，$k_{2y}$，$k_{2z}$，$k_{3x}$，$k_{3y}$，$k_{3z}$，$k_{4x}$，$k_{4y}$，$k_{4z}$是微分方程组中导数的值。
常微分方程数值解法
第七章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。
另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。
因此，有必要探讨常微分方程初值问题的数值解法。
本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。
第一节 欧拉法求解常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy（1）的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<<<n x x x x 210 上的近似值 ,,,,,210n y y y y{}n y 称为问题的数值解，数值解所满足的离散方程统称为差分格式，1--=i i ix x h 称为步长，实用中常取定步长。
显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。
定理 设函数()y x f ,在区域+∞≤≤-∞≤≤y b x a D ,:上连续，且在区域D 内满足李普希兹(Lipschitz)条件，即存在正数L ，使得对于R 内任意两点()1,y x 与()2,y x ,恒有()()2121,,y y L y x f y x f -≤-则初值问题(1)的解()x y 存在并且唯一。
一、欧拉法(欧拉折线法)若将函数)xy 在点nx处的导数()n x y '用两点式代替， 即()hx y x y x y n n n )()(1-≈'+，再用n y 近似地代替()n x y ,则初值问题(1)变为⎩⎨⎧==++=+ ,2,1,0),()(001n x y y y x hf y y n n n n(2)(2)式就是著名的欧拉(Euler)公式。
以上方法称 为欧 拉法或欧拉折线法。
常微分方程组的数值解
dx xxi
x xi
由数值微分公式，我们有 向前差商公式
yi1 h
yi
f (xi , yi )
yi1 yi h f (xi , yi )
可以看到，给出初值，就可以用上式求出所有的 {yi}
求的是在一系列离散点列上，求未知函数y在这些点上的值的近 似。称为数值离散方法。
基本步骤如下：
① 对区间作分割：I : a x0 x1 xn b
(1)
y(a) y0
由常微分方程理论可知：只要上式中的函数f(x,y)在区域
G={a≤x≤b,－∞＜y＜∞}内连续，且关于 y 满足Lipschitz条件，即存 在与 x, y 无关的常数L，使
f (x, y1) f (x, y2 ) L y1 y2
对G内任意两个y1, y2 都成立，则方程(1) 的解y y(x)必定存在且唯一。
这类形式的方法也称为差分方法。
当假定 yi为准确值，即在 yi y(xi )的前提下来估计误差
y y ( xi1 )
i 1
，这种截断误差称为局部截断误差。
由（2）、（3）知Euler公式在 xi处的局部截断误差为：
Ri
y(xi1)
yi1
h2 2
y"( i ), (xi
i
xi1), (i
0,1,n 1)
① 收敛性问题 步长充分小时，所得到的数值解能否逼近问题的真解；
② 误差估计
③ 稳定性问题 舍入误差，在以后各步的计算中，是否会无限制扩大；
对于初值问题（1），先将其离散化，即把[a,b]区间n等分， 得各离散节点
xi a ih
(i 0,1,2,n)
其中h b a . n
设y yx为方程的解，则y(xi1 )在(xi , yi )处
常微分方程数值解
p
定义：一个方法的总体截断误差若为O h ， 则称之为 p阶方法。
一般地，方法的总体截断误差阶越高，精度 也越高。
误差分析表
Euler方法 局部截 断误差 向前Euler O(h2) 方法 向后Euler O(h2) 方法 O(h3) 梯形公式 总体截 断误差 O(h)
O(h) O(h2)
微分方程离散化常用方法
用差商代替微商 数值积分 Taylor展开
A 用差商代替微商 y xn1 y xn dy f ( xn , y( xn )) dx x n, y n xn1 xn
yn 1 yn 则 f xn , yn h yn 1 yn hf xn , yn n 0, 1, 2,
迭代收敛 条件
0<hL<1 0<hL<2
(L为Lip常数）
向后Euler 方法收敛条件与截 断误差
局部截断误差
T
O n 1
h
2
，
整体截断误差 O h (当 0 hL 1 时)
y y
n 1
k 1
k
n 1
h f
k 1
n 1
x
,y n 1
k
取左端点的矩形公式
则有Leabharlann xn 1xnf ( x, y )dx h f ( xn , yn )
yn 1 yn hf ( xn , yn ) (n 0,1,)
用yn 1 , yn 代替y ( xn 1 ), y ( xn ), 对右端积分采用 取右端点的矩形公式
则有
xn 1
其中hk=xk+1-xk , 如是等距节点h=(b-a)/n ,
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
科学计算—理论、方法及其基于MATLAB 的程序实现与分析微分方程（组）数值解法§1 常微分方程初值问题的数值解法微分方程（组）是科学研究和工程应用中最常用的数学模型之一。
如揭示质点运动规律的Newton 第二定律：()()()⎪⎩⎪⎨⎧'='==000022x t x x t x t F dt xd m （１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：()()00y y t f ay dtdy=+= （２） 的解为：()()()τττd f e y e t y tt a at⎰-+=00 （３）但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。
一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题：()()000,y t y t t t y t F dtdyf=≤≤= （７）其中()()()()⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=t y t y t y t y n 21 （８）()()()()⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=y t f y t f y t f y t F n ,,,,21（９） 常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。
§１.1 常微分方程（组）的Cauch 问题数值解法概论假设要求在点（时刻）k t ，n k ,,2,1 =处初值问题（７）的解()00,,y t t y y =的（近似）值,如果已求得k t 时刻的值()k t y 或它的近似值k y (如0t 时刻的值()0t y )，那么将式（７）的两端在区间[]1,+k k t t 上积分()()()()dt t y t f t y t y dt dtdyk kk kt t k k t t ⎰⎰++=-=+11,1 (10)可得()()()()dt t y t f t y t y k kt t k k ⎰++=+1,1 (11)或()()()dt t y t f y y t y k kt t k k k ⎰++=≈++1,11 (12)显然,为了利用式（11）或（12）求得()1+k t y 的精确值（近似值1+k y ），必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解, 可以说,不同的式值积分方法将给出不同的Cauch 问题的数值解法。
§1.2 最简单的数值解法——Euler 方法假设要求在点（时刻）kh t t k +=0，nt t h f 0-=，n k ,,2,1 =处初值问题（７）的解()t y y =的近似值。
首先对式（７）的两端积分，得()()()()dt t y t f t y t y dt dtdyk k k kt t k k t t ⎰⎰++=-=+11,1 （13） 对于式（13）的右边，如果用积分下限k t t =处的函数值()()k k t y t f ,代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面 积），则有()()()()()[]k k t tk k t y t hf dt t y t f t y t y k k,,11≈=-⎰++ （14）进而得到下式给出的递推算法—Euler 方法()()001,,2,1,y t y nk y t hf y y k k k k ==+=+ （15）例１ 用Euler 方法解如下初值问题，取3.0=h ，()1030sin 23=≤≤-=y t t y y dtdy解：由（15）得()1010,,2,1sin 6.03.131==-=+y k t y y y kk k k结果如下：如果取1.0=h ，其结果如下图所示：Euler_Method§1.3 改进的Euler 方法对于（15）的右边，如果被积函数用积分限k t t =和1+=k t t 处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有()()()()()()()()[]111,,2,1++++≈=-⎰+k k k k t t k k t y t f t y t f hdt t y t f t y t y k k（16） 进而得到下式给出的递推算法：()()[]()00111,,2,1,,2y t y nk y t f y t f hy y k k k k k k ==++=+++ （17）通常算法（17）比Euler 方法（15）的精度高，但是，按算法（17）求1+k y 时要解（非线性）方程（组），这是算法（17）不如Euler 方法的方面，为了１） 尽可能地保持算法（17）精度高的优点； ２） 尽可能地利用Euler 方法计算简单的长处； 人们采取了如下的称之为改进的Euler 方法的折衷方案：预测 ()k k k k y t hf y y ,01+=+ （18） 修正 ()()[]n k y t f y t f h y y k k k k k k ,,2,1,,2111 =++=+++ （19）()00y t y =例２ Euler 方法与改进的Euler 方法的比较 下图是当3.0=h 时比较的结果：open Improved_Euler_Method.m§1.4 Euler 方法和改进的Euler 方法的误差分析由Taylor 公式()()()()()()()()()()22111,!21h O h t y t f t y t t t y t t t y t y t y k k k k k k k k k k k ++=+-''+-'+=+++（19） 说明Euler 方法的截断误差是()2h O ，类似地，由()()()()()()321!21,h O h t y h t y t f t y t y k k k k k +''++=+ （20） ()()()()()()321111!21,h O h t y h t y t f t y t y k k k k k +''+-=++++ （21） 以及()()()()()()h O t y t y h t y t y t y k k k k k +''=''⇒+'''+''=''++11 (22)让式（20）的两端减式（21）的两端，可得()()()()()()[]()()()()()()()[]()31132111,,241,,2h O t y t f t y t f ht y h O h O h t y t f t y t f ht y t y k k k k k k k k k k k +++=++++=+++++ （23）从上述推导Euler 方法、改进的Euler 方法的过程以及例１、例２容易看出，改进的Euler 方法Euler 方法的精度高，其原因在于： 1 在推导Euler 方法时，我们是用待求解函数()t y y =在一点处的变化 率()()k k t y t f ,代替()t y y =在区间],[1+k k t t 上的平均变化率：()()()()t y t hf dt t y t f k kt t ,,1=⎰+ （24）2 而在推导改进的Euler 方法时，我们是用待求解函数()t y y =在两点处变化率的平均值()()()()[]11,,21+++k k k k t y t f t y t f 代替()t y y =在区间],[1+k k t t 上的平均变化率；显然，通常()()()()[]11,,21+++k k k k t y t f t y t f 比()()k k t y t f ,更接近于()t y y =在区间],[1+k k t t 上的平均变化率。
由此启发人们：适当地选取区间],[1+k k t t 上函数()t y y =若干点处的变化率，用它们加权平均值代替()t y y =在区间],[1+k k t t 上的平均变化率，近似解的精度应更高。
下面将要介绍的Runge —Kutta 法就是基于上述想法得到的。
§2 Runge —Kutta 法Runge —Kutta 法是按选取区间],[1+k k t t 上函数()t y y =变化率()y t f ,的个数的多少和截断误差的阶数()m h O 来区分的一系列方法，如 １ 二阶的Runge —Kutta 法（改进的Euler 方法）()()()⎪⎪⎩⎪⎪⎨⎧++=+==++21111212,,K K hy y hK y t f K y t f K k k k k k k （25） ２ 三阶的Runge —Kutta 法()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧+++=++=⎪⎭⎫ ⎝⎛++==+32112312146,2,2,K K K h y y hK y h t f K K h y h t f K y t f K k k k k k k k k （26） ３ 四阶的Runge —Kutta 法 1) 古典形式()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++++=++=⎪⎭⎫ ⎝⎛++=⎪⎭⎫ ⎝⎛++==+432113423121226,2,22,2,K K K K hy y hK y h t f K K h y h t f K K h y h t f K y t f K kk k k k k k k k k （27） 2) Gill 公式(具有减小舍入误差的优点)()()()()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+++-++=⎪⎪⎭⎫ ⎝⎛++-+=⎪⎪⎭⎫ ⎝⎛-+-++=⎪⎭⎫ ⎝⎛++==+432113242131212222622222,222212,22,2,K K K K hy y hK hK y h t f K hK hK y h t f K K h y h t f K y t f K k k k k k k k k k k （28） 4 Runge —Kutta 法的一般形式()∑=++=+=ni i i k k k k k K a h y h y t h y y 11,,φ (29)其中()h y t k k ,,φ称为增量函数(Increment Function)以及()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧++++=+++=++==-----11,111,1122212123111121,,,,n n n n k n k n k k k k k k hK q hK q y h p t f K hK q hK q y h p t f K hK q y h p t f K y t f K （30）需要特别指出的是：在确定（29）、（30）中的参数i a ，i p 和ij q 时，应该使（29）的右端和适当阶的（如n 阶）Taylor 展式一致，这样，至少对于底阶的R-K 法来说，参与加权的斜率个数n 与方法的阶数是一致的。