hk算法

A. hK是什么硬度

HK是努氏硬度值单位的表示符号，用于努氏硬度。

努氏硬度是以发明人美国的Knoop命名的。过去又称克氏硬度、克努普硬度、努普硬度。在我国列为小负荷硬度试验，也可称其为显微硬度。它与显微维氏硬度试验一样，使用较小的力以特殊形状的压头进行试验。测量压痕对角线求得硬度值。

(1)hk算法扩展阅读：

努氏硬度特点

1、由于使用压头的特殊形状，试验时能产生长短对角线为7:1的棱形压痕。原理规定只测量长对角线长度，具有较高的测量精度。

2、努氏硬度试验的压痕压入深度只有长对角线长度的1/30，维氏硬度试验的压痕压入深度为对角线长度的1/7，所以努氏硬度试验适用于表层硬度和薄件的硬度测试。

3、同一试样在同一负荷下，努氏硬度压痕对角线长度约为维氏硬度压痕对角线长度的3倍，大大优于维氏测量法。

B. Mallat算法

6.9.1.1 尺度空间的有限分解

MRA框架表明f（t）∈L²（R）可分解为无穷个小波分量的直和，但在实际应用中，仅知道f（t）的近似函数。为不失一般性，可假设原信号是在a=1或j=0的分辨率下测得的，用f⁰（t）表示，它属于子空间V₀。而子空间V₀又可分解成两个子空间。因此，在MRA框架下理解为f⁰（t）∈V₀，这样就有如下的尺度空间的有限分解表现

V₀=V₁⊕W₁=（V₁⊕W₁）⊕W₀=（V₂⊕W₂）⊕W₁

=（V_j⊕W_j）⊕W_j-1⊕W_j-2⊕…⊕W₁

︙

=V_J+W_J⊕W_J-1⊕W_J-2⊕…⊕W₁

其中子空间及分量分别为

地球物理信息处理基础

V₀的有限分解关系仅是MRA无穷分解中的一部分。因此，就子空间而言、就函数分量而言以及就频率范围而言，有限分解的含义都与MRA相同。例如，V₀中的元素f⁰（t）是有限频率范围的，f^J（t）∈V_J是f⁰（t）的最低频表现，W_j中元素δ^j（t）是具有特定带宽的，它们互不重叠，这些频带的总和就是f⁰（t）的频率范围。

为了数字计算和分析处理的目的，需要将f^j（t）和δ^j（t）用离散数据来表示。显然，｛c_j，k｝∈l²是合适的，因为

表明f^j（t）∈V_j和｛c_j，k｝∈l²是一一对应的，或者说利用｛c_j，k｝就一定可恢复f^j（t）。根据同样的理由，可用数据｛d_j，k｝∈l²来表示δ^j（t）∈W_j，这类数据｛c_j，k｝和｛d_j，k｝对应着特定的离散表现，j表示尺度方面的2进制规则离散，k表示整节点标记的时移离散。

6.9.1.2 分解算法

分解算法要实现的目标是：在｛φ（t-k）｝是标准正交基条件下，已知｛c_j-1，k｝、｛h_k｝和｛g_k｝，求出｛c_j，k｝和｛d_j，k｝。

ψ（t）关于φ（t）的两尺度关系式（6-95）提供了一条由尺度函数φ（t）构造母小波ψ（t）的途径。根据尺度函数

地球物理信息处理基础

和小波函数

地球物理信息处理基础

得到

地球物理信息处理基础

用l取代2k+l，上式成为

地球物理信息处理基础

计算f（t）∈L²（R）与上式两端的内积，便得到如下计算小波级数系数的一个公式

地球物理信息处理基础

式中c_j-1，l=＜f^j，φ_j-1，l＞=＜f，φ_j-1，l＞是信号f（t）与φ_j-1，l（t）的内积，即c_j，l=＜f，φ_j，l＞。信号在V_j的正交投影

，将两尺度关系式的φ_j，k（t）代入定义式φ_j，k（t）=2^-j/2φ（2^-jt-k），有

地球物理信息处理基础

由此可以计算f（t）与φ_j，k（t）的内积，即

地球物理信息处理基础

在数学上，f^j（t）是平方可积连续函数，它是L²（R）的元素；c_j，l（l∈Z）是平方可和序列，它是平方可和序列矢量空间l²（Z）中的元素。

说明f^j（t）与c_j，k（k∈Z）有等效对应关系，在数学上称l²（Z）与L²（R）之间存在着同构关系。显然，c_j，k比f^j（t）更适合于在计算机上运算。c_j，k称为在分辨2^j下的f（t）离散逼近，因此，可以把f⁰（t_k）=f（t_k）的采样值作为｛c_j，k｝的初始数据，即｛c_0，k｝=｛f（t_k）｝。

综合式（6-109）和（6-110），得到一般的分解公式

地球物理信息处理基础

式（6-111）的实现过程如图6-26所示，每个尺度j所存储的数据｛c_j，k｝都是按整数编号的。以j尺度层为基础来观察j-1尺度层，j尺度层的采样节点编号k对应着j-1尺度层上编号为2k的采样节点，或者说。j-1尺度的采样节点是在j尺度采样节点基础上均匀加密的结果；若以j-1尺度层为基础来观察j尺度层，则j-1尺度层上隔2取样（“隔1取1”）的节点正好对应着j尺度层上的采样节点。图6-26表明了由j-1尺度向j尺度的变换过程，｛h_k｝可看做滤波器的单位冲激响应（权系数）。假设｛h_k｝仅有6个元素，则式（6-111）所表明的变换过程相当于把｛h_-2，h_-1，h₀，h₁，h₂，h₃，｝作为权值，其中心点h₀对准｛c_j，k｝后再作加权平均，即

地球物理信息处理基础

这一实现过程是非常快捷的，相当于先计算序列c_j-1，k（k∈Z）与

的卷积，然后抽取偶数下标（隔1抽1）的卷积结果，即得到逼近级数的系数c_j，k（k∈Z）。

图6-26 分解算法计算｛c_j，k｝示意图

图6-26虽然仅表明由｛c_j-1，k｝计算｛c_j，k｝的变换过程，它同样也表明了由｛c_j-1，k｝计算｛d_j，k｝的变换过程，只不过要将｛h_k｝换成｛g_k｝而已。对正交小波而言，g_k=（-1）^kh_-k+1是由｛h_k｝决定的。

由c_0，k开始，利用式（6-111）进行迭代运算，陆续计算出c_1，k、c_2，k等等，与此同时，利用c_0，k、c_1，k、c_2，k等值，同样不断计算出d_1，k、d_2，k等小波级数系数值。

式（6-111）所表明的计算过程由算子表示会更简单些。记A^j=｛c_j，k｝，D^j=｛d_j，k｝。记算子H：l²➝l²，其运算意义如式（6-111）的第1式所示，即

地球物理信息处理基础

同样，记算子G：l²➝l²，其运算意义如式（6-111）的第2式所示，即

地球物理信息处理基础

采用算子表示后，式（6-111）所表明的分解算法结构如图6-27所示。它表明

A^j=H^jA⁰；D^j=GH^j-1A⁰

即A⁰经H算子j次作用后即可获得A^j，A⁰经H算子j-1次作用后再经G算子作用1次即可获得D^j。

图6-27 分解算法结构示意图

图6-27还表明，只要在细密采样间隔的尺度层次上给定A⁰，就可利用分解算法快速地获得较粗采样间隔尺度层上的有关数据A^j和D^j（0＜j≤J）。假设实际问题中A⁰有N个数据，则2尺度层上的A¹和D¹各有N/2个数据，依此类推；还假设｛h_k｝和｛g_k｝分别有M个数据；那么，用A⁰计算A¹和D¹共需2MN/2次运算，从2到3尺度层共需2MN/4次运算，依此类推可知，要得到A^j和D^j（0＜j≤J）共需2MN（2^-1+2^-2+…2^-J）次运算。由此可见，分解算法是快速的。

原信号f⁰（t）具有下列唯一分解

地球物理信息处理基础

δ_j是信号在W_j（j=1，2，…，J）各子空间上的正交投影，它们是从一个较精细的逼近变成较粗略的逼近（两个逼近的分辨率相邻近）时所丢失的信息，即

地球物理信息处理基础

可将这一不断降低逼近分辨率的过程看成是“一层又一层地把信号进行剥皮”的过程。当J选得足够大时，“剥”下来的信息总和

地球物理信息处理基础

足够多，将足以精确表示原信号f⁰（t），而最终的逼近信号f^J（t）的分辨率已经非常低，这样反而可以把式（6-113）当成原信号的估计，而把f^J（t）看成是估计误差。这就是说，用小波级数在所有分辨率下的全部系数（j=1，2，…，J）来代替原信号，其误差f^J（t）可以任意小。按照这种解释，式（6-111）算法就是将f⁰（t）的信息（c_0，k是它的离散表示）表示成c_1，k～c_J，k等信息和一个估计误差（实际上它是在分辨率最低即2^J下的逼近）c_J，k。这一过程实际上是在一次又一次地改变着正交基（或子空间）。

6.9.1.3 信号重构算法

重构算法是分解算法的逆过程。此时已知数据｛c_j，k｝和｛d_j，k｝（0≤j≤J），希望利用这些数据快速准确地重构出原始数据。

一般而言，相邻两分辨率下的逼近信号存在着下列关系：

地球物理信息处理基础

计算上式左右两端与φ_j-l，k（t）的内积，得

地球物理信息处理基础

为了便于讨论运算过程，在上式中令m=p-2k，所以信号重构公式为

地球物理信息处理基础

根据式（6-115）计算｛

｝的过程见图6-28所示。为了理解式（6-115）和图6-28所示的计算过程，应注意下面两点：一是j-1尺度层偶数编号的采样点对应着j尺度层上的采样点，二是示意图中假设｛h_k｝的有限数据为｛h_-2，h_-1，h₀，h₁，h₂，h₃，｝。

可以看出，计算｛

｝可分为两步进行：

地球物理信息处理基础

第一步计算j-1尺度层上的偶数编号采样点处的｛

｝，此时可对j尺度层采样节点循序推进，所用数据为｛h_-2，h₀，h₂，｝，按公式

地球物理信息处理基础

计算出j-l尺度层上偶数编号节点处的相应值。

第二步计算j-1尺度层上的奇数编号的采样点处的｛

｝，其循序推进过程仍同于第一步，但需另用｛h₃，h₁，h_-1｝，且按公式

地球物理信息处理基础

图6-28也表明了式（6-115）中｛d_j，k｝计算｛

｝的过程，此时要采用｛g_k｝，计算式应为

地球物理信息处理基础

为了更清楚地表明｛c_j，p｝和｛d_j，p｝，0≤ j≤J-1重构｛c_0，p｝的关系，与前述的分解情况一样，记A^j=｛c_j，p｝，D^j=｛d_j，p｝，采用算子的记号H^*：l²➝l²，其运算意义为

地球物理信息处理基础

记G^*：l²➝l²，其运算意义为

地球物理信息处理基础

在这些重构算子意义下，式（6-115）可表示为

地球物理信息处理基础

于是从尺度层j=J到尺度层j=0的重构算法过程可用图6-29表示。重构算法也是快速的，现估计其计算量。设｛h_m｝和｛g_m｝分别有M个数据，0尺度层｛c_j，0｝有N个数据，则尺度层1上｛c_j，1｝有N/2个数据，用｛c_1，p｝计算｛

｝的偶节点值和奇节点值分别需要（M/2）N/2次运算；由｛c_1，p｝计算｛

｝共需MN/2次运算；｛d_1，p｝计算｛

｝也需MN/2次运算；由1尺度层数据重构0尺度层数据共需MN次运算。可依此类推出其它相邻尺度层之间采用重构算法的运算量，最后，得总运算量为

2MN（2^-1+2^-2+…2^-J）

图6-29 信号重构算法示意图

图6-30所示的是信号分解与重构算法的计算流程示意图，图中的“

”表示“按因子2抽取”或“隔1抽1”运算。符号“

”表示，“按因子2内插”或“两相邻样本间内插一个零样本”。输入是c_0，k，相当于原信号的离散表示；输出是不同尺度下的小波级数系数c_j，k（j∈Z）以及信号在子空间V_J的正交投影c_J，k。这是在做信号分析，其分解公式见（6-111）。与此相反，当输入是d_j，k（j∈Z）和c_J，k，输出是c_0，k时，为信号重建或合成，此时有

地球物理信息处理基础

将合成与分析算法合起来画成图6-30所示的形式，称之为Mallat算法。

图6-30 信号分解与重构算法示意图

6.9.1.4 Mallat算法实现中的一些问题

Mallat算法是一种纯数字的快速递推算法，在使用Mallat算法时，有一些具体问题需引起注意。

（1）对正交尺度函数φ（t）而言，Mallat算法中仅需数据｛c_0，k｝和｛h_k｝可进行快速的分解和重构递推运算。要存储的数据为｛c_J，k｝和｛d_j，k｝，0＜j≤J，这些有用的数据的存储量等于｛c_0，k｝的数据存储量。特别值得强调的是，Mallat算法中隐含着两类关系，一类是关于多分辨分析方面的，例如对0＜j≤J，有

地球物理信息处理基础

f^j-1（t）=f^j（t）+δ^j（t）

地球物理信息处理基础

另一类关系是由于尺度函数φ（t）平移正交性产生的，例如

c_j，k=＜f，φ_j，k＞；d_j，k=＜δ^j，ψ_j，k＞

g_k=（-1）^kh_-k+1，0＜j≤J

Mallat算法正是利用了这些关系，在算法实施过程中不需尺度函数φ（t）和小波函数ψ（t）的具体形式，只要求它们存在并找出｛h_k｝，就可以顺利地进行分解和重构处理了。因此，只要查得正交尺度函数双尺度方程的传递系数｛h_k｝，就可以应用Mallat算法了。

顺便指出，如果尺度函数φ（t）（例如样条函数）不是平移正交的，它虽然可以生成MRA，但由此构造出的样条小波ψ（t）仅关于尺度正交，没有平移正交性。此时

地球物理信息处理基础

所以Mallat算法不再适用，必须另行推导相应的分解和重构算法。

（2）初始数据｛c_0，k｝的选用，在正交小波分解中，φ（t）是正交尺度函数，0尺度层上的展180开系数｛c_0，k｝=＜f（t），φ_0，k（t）＞，用复杂的计算来确定初始数据｛c_0，k｝是不合算的，应采用变通处理办法，即简单采用｛c₀，k｝=｛f（t_k）｝，用细尺度层上的采样值作为初始数据｛c_0，k｝。这种做法似乎有些不严密，但可以证明，虽然｛f（t_k）｝作为｛c_0，k｝的近似值时有微小误差，但数据｛f（t_k）｝同样有效地表现了f⁰（t）的变化波动状况和有效的频率范围，这种替代不会影响对f⁰（t）的时频分析；同时还应看到，用｛f（t_k）｝作初始数值，不仅可使问题简单化，而且也可使Mallat算法准确地分解和重构初始数据。总之，用｛f（t_k）｝替代｛c_0，k｝是实用方便的。

（3）分解层数和采样间隔的关系，这个问题主要从以下几方面考虑便可得出结论。

第一，因为最细的0尺度层的采样间隔T决定了f⁰（t）的频率范围，由取样效应可知，最大的频率范围为｜ω（f⁰）｜≤1/（2T）；同样，最粗的J尺度层的采样间隔为J=2^JT，最低的频率范围为｜ω（f^J）｜≤1/（2T）=2^-J｜ω（f⁰）｜。于是可从需要分辨的最高频率和需要分辨的最低频率这两个指标来决定最细尺度层的采样间隔和数据分解的层数。

第二，在最细的0尺度层上，应取用多少个数据才能满足J个层次的数据分解呢？在Mallat算法中，｛c_1，k｝的数据量仅为｛c_0，k｝数据量的一半，依此类推。同样，在J尺度层至少要取用N_J个数据才能表现低频量，于是推知，在0尺度层，至少取用N=2^JN_J个数据，才能满足J个分解层次的需要。

（4）在Mallat算法的运算中，需用到所存储的数据外面的数据（见图6-31），图中实线框内数据是要存储的。在分解算法中，若｛h_k｝仅有6个数据，由图6-26可知，要用到实线框外的数据c_0，-2和c_0，-1才能计算出c_1，-1，要用到c_0，N0+1、c_0，N0+2和c_0，N0+3才能计算出c_1，N1，其它层次的情形类似。由于细密采样层（图6-31中对应着j=0的尺度层次）中的近似函数f⁰（t）∈L²（R），当t➝∞时，f⁰（t）➝0，实线框内足够多的采样数据｛c_0，k｝，k=0，1，…，N₀，已反应了f⁰（t）的基本特性，f⁰（t）在图6-31实线框外的数据几乎消失，因此在实际分解过程中可简单地令实线框外的数据为零。同样，重构过程也会用到实线框以外的数据（见图6-28），也可以简单地令实线框外的数据为零。

图6-31 Mallat算法涉及没有存储的数据示意图

另一种办法也是常用的。在最细尺度层上较多地取用数据，在计算过程中适当地多存储些数据，如图6-31中虚线框所示。此时应以实线框存储数据作为分解重构算法和进行数字分析的依据。

6.9.1.5 Mallat 算法所表现的频域分解特点

有限尺度空间的正交小波子空间的直和分解关系，例如，

V₀=V₁⊕W₁=V₂⊕W₂⊕W₁=V₃⊕W₃⊕W₂⊕W₁

在Mallat算法中是通过算子H和G来表现的。数据A^j表征f^j（t）∈V_j，数据D^j表征δ^j（t）∈W_j，那么在Mallat算法中，A^j=HA^j-1；D^j=GA^j-1，从而实现了子空间的分解。这种子空间分解和算子H与G之间的关系示意于图6-32中（以V₀分解成3层为例）。

图6-32 Mallat算法所确定的数据分解和子空间分解的对应关系

因为f^j-1（t）和f^j（t）都是有限频率范围的，f^j（t）的频率范围仅是f^j-1（t）的相对低频部分，δ^j（t）的频率范围是f^j-1（t）的相对高频部分，所以δ^j（t）是f^j-1（t）的关于高频带的分量。换成子空间来描述，V_j-1所表现的频率范围被分解为两部分，一部分是由V_j表现的低频部分，另一部分是由W_j表现的频带部分。因此，V₀尺度空间的直和分解，表明把f⁰（t）的频率范围分解为由W₁、W₂、W₃和V₃所表现的频带，且这些频带是互不重叠的。这里所说的W₂所表现的频带宽度，也就是其基函数ψ_2，k（t）（窗函数）的频窗宽度。由aD_ω可知，W₂所表现的频带宽是W₁的一半，尺度指标增加，则小波子空间所表现的频带宽分半减小。由于对任意尺度而言，小波子空间所表现的时频窗面积是恒定不变的常数，所以对时频窗的时窗宽度而言，尺度指标增加，则相应的时窗宽度成倍增加。图6-33表示了正交小波分解用于分析t^*处局部时域信号时，在各个频带的时频窗表现。

图6-33 正交小波分解在各频带的时频窗示意图

C. patricia tree是什么

Patricia Tree 简称PAT tree。它是 trie 结构的一种特殊形式。是目前信息检索领域应用十分成功的索引方
法，它是1992年由Connel根据《PATRICIA——Patrical Algorithm to Retrieve Information Coded in Alphanumeric》算法发展起来的。

PAT tree 在字符串子串匹配 上有这非常优异的表现，这使得它经常成为一种高效的全文检索算法，在自然语言处理领域也有广泛的应用。其算法中最突出的特点就是采用半无限长字串(semi-infinite string 简称 sistring) 作为字符串的查找结构。

采用半无限长字串(sistring): 一种特殊的子串信息存储方式。

比如一个字符串CUHK。它的子串有C、CU、CUH、CUHK、U、UH、UHK、H、HK、K十种。如果有n个字符的串，就会有n(n+1)/2种子串，其中最长的子串长度为n。因此我们不得不开辟 n(n+1)/2个长度为n的数组来存储它们，那么存储的空间复杂度将达到惊人的O(n^3)级别。

但是我们发现这样一个特点：

CUHK —— 完全可以表示 C、CU、CUH、CUHK

UHK —— 完全可以表示 U、UH、UHK

HK —— 完全可以表示 H、HK、

K —— 完全可以表示 K

这样我们就得到了4个sistring: CUHK、UHK、HK和K。

PAT tree的存储结构

如果直接用单个字符作为存储结点，势必构造出一棵多叉树(如果是中文字符的话，那就完蛋了)。检索起来将会相当不便。事实上，PAT tree是一棵压缩存储的二叉树结构。现在我们用“CUHK”来构造出这样一棵PAT tree 。

开始先介绍一下PAT tree的结点结构(看了后面的过程就再来理解这些概念)

* 内部结点：用椭圆形表示，用来存储不同的bit位在整个完整bit sequence中的位置。

* 外部节点(叶子结点)： 用方形表示，用来记录sistring的首字符在完整sistring中的开始位置(字符索引)和sistring出现的频次。

* 左指针：如果 待存储的sistring在 内部结点所存储的bit位置上的数据 是0，则将这个sistring存储在该结点的左子树中。

* 右指针：若数据是1，则存储在右子树中。

(1) 将所有sistring的字符转化成1 bytes的ASCII码值，用二进制位来表示。形成一个bit sequence pattern(没有的空字符我们用0来填充)。

sistring bit sequence

完整sistring -> CUHK 010 00011 01010101 01001000 01001011 <- 完整bit sequence

UHK0 010 10101 01001000 01001011 00000000

HK00 01001000 01001011 00000000 00000000

K000 01001011 00000000 00000000 00000000

(2) 从第一个bit开始我们发现所有sistring的前3个bit位都相同010，那么相同的这些0/1串对于匹配来说就毫无意义了，因此我们接下来发现第4个bit开始有所不同了。UHK 的第4个bit是1，而CUHK、HK、K的第4个bit是0。则先构造一个内部结点iNode.bitSize=4（第4个bit），然后将UHK的字符索引 cIndex=2(UHK的开始字符U在完整的CUHK的第2位置上)构造成叶子结点插入到iNode的左孩子上，而CUHK、HK、K放在iNode右子树中。(如下图2)

(3) 递归执行第2步，将CUHK、HK、K进一步插入到PAT tree中。流程如下图所示。所有sistring都插入以后结束。

注意：既然PAT tree 是二叉查找树，那么一定要满足二叉查找树的特点。所以，内部结点中的bit 位就需要满足，左孩子的bit 位< 结点bit 位< 右孩子的bit 位。

PAT tree的检索过程

利用PAT tree可以实现对语料的快速检索，检索过程就是根据查询字串在PAT tree中从根结点寻找路径的过程。当比较完查询字串所有位置后，搜索路径达到PAT tree的某一结点。

若该结点为叶子结点，则判断查询字串是否为叶子结点所指的半无限长字串的前缀，如果判断为真，则查询字串在语料中出现的频次即为叶子结点中记录的频次；否则，该查询字串在语料中不存在。

若该结点为内部结点，则判断查询字串是否为该结点所辖子树中任一叶子结点所指的半无限长字串的前缀。如果判断为真，该子树中所有叶子结点记录的频次之和即为查询字串的出现频次。否则，查询字串在语料中不存在。

这样，通过PAT tree可以检索原文中任意长度的字串及其出现频次，所以，PAT tree也是可变长统计语言模型优良的检索结构。

例如：要查找string= “CU ”(bit sequence=010 00 0 1 1 01010101) 是不是在CUHK 中。

(1) 根据“CUHK ”的PAT tree 结构( 如上图) ，根结点r 的bit position=4 ，那么查找bit sequence 的第4 个bit=0 。然后查找R 的左孩子rc 。

(2) rc 的bit position=5 ，在bit sequence 的第5 个bit=0 。则查找rc 的左孩子rcc 。

(3) rcc= ” CUHK ” 已经是叶子结点了，则确定一下CU 是不是CUHK 的前缀即可。

PAT tree 的效率

特点：PAT tree查找的时间复杂度和树的深度有关，由于树的构造取决于不同bit位上0,1的分布。因此PAT tree有点像二叉查找树 ，最坏情况下是单支树(如上图例子)，此时的时间复杂度是O(n-1)，n为字符串的长度。最好情况下是平衡二叉树 结构，时间复杂度是O(log2(N))。另外，作为压缩的二叉查找树，其存储的空间代价大大减少了。

PAT tree的实际应用

PAT tree在子串匹配上有很好的效率，这一点和Suffix Tree(后缀树)，KMP算法的优点相同。因此PAT tree在信息检索和自然语言处理领域是非常常用的工具。比如：关键字提取，新词发现等NLP领域经常使用这种结构。

D. Runge--kutta算法

龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。

对于一阶精度的欧拉公式有： yi+1=yi+hki

其中h为步长，则yi+1的表达式与y(xi+1)的Taylor展开式的前两项完全相同，即局部截断误差为O(h2)。 当用点xi处的斜率近似值k1与右端点xi+1处的斜率k2的算术平均值作为平均斜率k∗的近似值，那么就会得到二阶精度的改进欧拉公式： yi+1=yi+h(k1+k2)

其中k1=f(xi,yi)，k2=f(xi+h,yi+hk1) 依次类推，如果在区间[xi,xi+1]内多预估几个店上的斜率值k1,k2,…,km，并用他们的加权平均数作为平均斜率k∗的近似值，显然能够构造出具有很高精度的高阶计数公式。 上述两组公式在形式删过的共同点：都是用f(x,y)在某些点上值得线性组合得出y(xi+1)的近似值yi+1，且增加计算的次数，可以提高截断误差的阶，他们的误差估计可以用f(x,y)在xi处的Taylor展开来估计。

于是可考虑用函数f(x,y)在若干点上的函数值的线性组合老构造金斯公式，构造时要求近似公式在f(xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式的前面几项重合，从而使金斯公式达到所需要的阶数。既避免求高阶导数，又提高了计算方法精度的阶数。或者说，在[xi,xi+1]这一步内计算多个点的斜率值，若够将其进行加权平均作为平均斜率，则可构造出更高精度的计算格式，这就是龙格-库塔(Runge-Kutta)方法。
一般的龙格-库塔法的形式为

#include"stdio.h"
#include"stdlib.h"
voidRKT(t,y,n,h,k,z)
intn;/*微分方程组中方程的个数，也是未知函数的个数*/
intk;/*积分的步数(包括起始点这一步)*/
doublet;/*积分的起始点t0*/
doubleh;/*积分的步长*/
doubley[];/*存放n个未知函数在起始点t处的函数值,返回时,其初值在二维数组z的第零列中*/
doublez[];/*二维数组,体积为nxk.返回k个积分点上的n个未知函数值*/
{
externvoidFunc();/*声明要求解的微分方程组*/
inti,j,l;
doublea[4],*b,*d;
b=malloc(n*sizeof(double));/*分配存储空间*/
if(b==NULL)
{
printf("内存分配失败 ");
exit(1);
}
d=malloc(n*sizeof(double));/*分配存储空间*/
if(d==NULL)
{
printf("内存分配失败 ");
exit(1);
}
/*后面应用RK4公式中用到的系数*/
a[0]=h/2.0;
a[1]=h/2.0;
a[2]=h;
a[3]=h;
for(i=0;i<=n-1;i++)
z[i*k]=y[i];/*将初值赋给数组z的相应位置*/
for(l=1;l<=k-1;l++)
{
Func(y,d);
for(i=0;i<=n-1;i++)
b[i]=y[i];
for(j=0;j<=2;j++)
{
for(i=0;i<=n-1;i++)
{
y[i]=z[i*k+l-1]+a[j]*d[i];
b[i]=b[i]+a[j+1]*d[i]/3.0;
}
Func(y,d);
}
for(i=0;i<=n-1;i++)
y[i]=b[i]+h*d[i]/6.0;
for(i=0;i<=n-1;i++)
z[i*k+l]=y[i];
t=t+h;
}
free(b);/*释放存储空间*/
free(d);/*释放存储空间*/
return;
}
main()
{
inti,j;
doublet,h,y[3],z[3][11];
y[0]=-1.0;
y[1]=0.0;
y[2]=1.0;
t=0.0;
h=0.01;
RKT(t,y,3,h,11,z);
printf(" ");
for(i=0;i<=10;i++)/*打印输出结果*/
{
t=i*h;
printf("t=%5.2f ",t);
for(j=0;j<=2;j++)
printf("y(%d)=%e",j,z[j][i]);
printf(" ");
}
}voidFunc(y,d)doubley[],d[];
{
d[0]=y[1];/*y0'=y1*/
d[1]=-y[0];/*y1'=y0*/
d[2]=-y[2];/*y2'=y2*/
return;
}