龙格库塔法
龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。

这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。

由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。

龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。

龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。

其中4 阶龙格库塔法是最常用的一种方法。

因为它相当精确、稳定、容易编程。

在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。

如果需要较高的精度, 可采取减小步长的方法即可。

4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。

1、初值问题
对于一阶常微分方程的初值问题
根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。

2、离散化
取步长h=(b-a)/n，将区间[a , b]分成n个子区间：
a=<=b
在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲
线上至少有一点，它的斜率与整段曲线的平均斜率相同，
得=y’() (0<<1)
其中，=
可以将上式改写成y()=y()+h*K (2.1)
其中K为平均斜率，K=f()
公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。

欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。

3、Euler法
欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙
格库塔法的基础。

首先，令、为y() 及y()的近似值，并且令平均斜
率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1）
4、改进的欧拉法
此种方法是取、两点的斜率的平均值作为平均斜率K，
即K= ,其中、均为y()以及y()的近似值，就得到
改进后的欧拉公式（4.1）
其中、分别为、两点的斜率值，即
= ,=
在上面的（4.1）式中，k2是未知的，采用一种叫预报法的方
法来求解。

即先用欧拉法求得一个初步的近似值，记作 ,称为预报值，预报值的精度不高，但是可以用它代替k2中的再直接计算，便可得到校正后的。

欧拉法求得预报值=，
校正=（）
将带入校正公式可得
改进后的欧拉公式
=（4.2）通常写作
（）
由此可见，为了提高精度，可以多取几点的斜率值作加权平均
来得到平均斜率K，这就是龙格库塔法的基本思想。

5、局部截断误差
初值问题
的单步法可用一般形式表示：
其中，多元函数与有关，如果中含有，那么方法是隐式的，若不含则为显式单步法，所以，显式单步法可以表示为
其中称为增量函数，对欧拉公式有
=f(x,y)
从开始计算，如果考虑每一步产生的误差，直到，则误差
，称该方法为在点的整体截断误差，分析求解整体截断误差是复杂的，为此，仅考虑从到的局部情况，并假设之前的计算没有误差，即，则一般显式单步法的局部截断误差可以定义为
一般，为了求解局部截断误差，我们会将对在
x=点做泰勒展开，进行化简。

6、龙格库塔法
龙格库塔法的基本思想就是用区间上的若干个点的导数
，将它们做加权平均得到平均斜率K，以提高方法的阶数，即提高方法的精度。

一般的显式龙格库塔法可以写成
其中，=
=f
（）
i=2,3......L
其中，，
它的局部截断误差是
其中，与的区别在于，用微分方程精确解代替中所
含的就得到，参数，，待定，确定这些常数的原则：使局
部截断误差误差关于h的阶数尽可能高。

二阶龙格库塔公式
令L=2，则
其局部截断误差为
（）（6.1）利用泰勒公式，将中各项作泰勒展开，特别是在
（，（））点做二元泰勒展开，并利用 ’，，=+y’()=y’’()
则有
+h*++O()
，
将上面各式带入（6.1）式中，整理得
选取、、使局部截断误差尽量小，就是使h和项的系
数为0，于是，得到方程组
这里是含有三个未知量的2个方程，它可以有无穷组解，都统
称为二阶龙格库塔公式。

取=0.5，则，=1，得到中点公式：
以上给出了构造2级龙格库塔公式的方法，它可以推广到L=3，4，5......的情况。

经典的龙格库塔法是一种4阶方法，计算公式是：
在具体计算中，利用已知的h，，，求得、、、
的值，代入公式中，就得到了的值。

7、变步长的龙格库塔法
单从每一步看，步长越小，截断误差就越小，但随着步长的减小，不但引起计算量的增大，而且可能导致舍入误差的严重积累，
因此，同积分的数值计算一样，微分方程的数值解法也要注意合理
的选择步长。

设从出发，以h为步长，经过一步计算，得到的近
似值，记作，若计算方法是p阶的，则有
（7.1）然后，将步长减半，即取为步长，从出发经两步计算，求得的近似值，记作
（7.2）当h很小时，以上两式中的c可看作同一常数，用乘式（7.2）的两端后与式（7.1）相减，得
因此，取=，显然要比
或的精度高。

当p=4时，可取
=
由此，可得下列近似式，
≈ =
上式表明，以作为的近似值，其误差可用先后两
次计算结果的差来表示。

可以用来判断所选取的步长是否合适，具
体可分为以下两种情况：
(1) 对于给定的精度，如果，就将反复将步长减半计算，直至为止，这时取步长折半后的“新值” 作为结果。

(2) 如果，就反复将步长加倍，直至为止，这时取上一次的步长的计算所得到的值作为。

这种通过步长减半或加倍来计算的方法就叫变步长法，表面上看，为了选择步长，每一步的计算量似乎增加了，但从总体上考虑往往是合理的。

8、龙格-库塔法在MATLAB中的用法
ode23:使用二阶龙格-库塔法求微分方程，调用格式为：
[t， ] od 3 ‘ n m ’ tsp n
ode45:使用四阶龙格-库塔法求微分方程，调用格式为：
[t， ] od 4 ‘ n m ’ tsp n
9、隐式龙格库塔公式
（
i=1,2,3......L
隐式龙格库塔法适用于解刚性问题，对于基于数值求积公式的隐式龙格库塔法，可以利用3种方法建立隐式龙格库塔方法，分别是基于高斯求积公式的隐式龙格库塔法，基于拉道求积公式的隐式龙格库塔法，基于洛巴托求积公式的隐式龙格库塔法。

而且一般的隐式龙格库塔法计算比较复杂，为了简化计算，可以考虑一种对角隐式龙格库塔法。

隐式龙格库塔法的稳定性函数
d t
d t ，这是z的有理函数。