常微分（ODE）方程求解方法总结
1 常微分方程（ODE）介绍
1.1 微分方程介绍和分类
1.2 常微分方程的非计算机求解方法
1.3 线性微分方程求解的推导过程
2 一阶常微分方程（ODE）求解方法
2.1 欧拉方法
2.1.1 欧拉方法
2.1.2 欧拉方法的误差分析
2.1.3 欧拉方法的改进思路1——添加高阶项
2.1.4 欧拉方法的改进思路2——对斜率采用更好的估值方法
2.1.4.1 修恩法
2.1.4.2 中点方法
2.1.5 欧拉方法的改进思路3
2.2 龙格-库塔法
2.2.1 二阶龙格-库塔法
2.2.2 三阶龙格-库塔法
2.2.2+ 牛顿积分公式
2.2.3 四阶龙格-库塔法
2.2.4 高阶龙格-库塔法
2.2.5 龙格-库塔法分析
3 一阶常微分方程组
3.1 欧拉方法
3.2 龙格-库塔法
4 自适应龙格·库塔法
4.1 步长控制中的误差估计方法
4.1.1 自适应龙格-库塔法或步长——对分法
4.1.2 龙格-库塔-费尔贝格法
4.1.3 其他改进思路
4.2 步长控制
1 常微分方程（ODE）介绍
1.1 微分方程介绍和分类
举例：假设跳伞人的下落速度v于时间有如下关系：
d v d t = g − c m v \frac{dv}{dt} = g-\frac{c}{m}v dtdv=g−mcv (1.1)
其中g为重力常数，m为质量，c为阻力系数。
被微分的量v因变量，与v有关的变量t称为自变量。 如果函数只有一个自变量，那么方程就称为常微分方程（ordinary differential equation，ODE）。 如果函数含有两个或者多个自变量，则成为偏微分方程（partialdifferential equation，PDE）。
此外，微分方程也可以根据阶数来分类：最高阶导数是一阶导数，则称为一阶微分方程（first-order-equation）；最高阶导数是二阶导数，则称为二阶微分方程（second-order-equation）。例如如（1.1）中就是一阶微分方程，下式（1.2）就是一个二阶微分方程。
m d 2 x d t 2 + c d x d t + k x = 0 m\frac{d^2x}{dt^2} +c\frac{dx}{dt} + kx = 0 mdt2d2x+cdtdx+kx=0 (1.2)
高阶微分方程能简化成一阶方程组。考虑上式（1.2），定义新变量y，令 y = d x d t y=\frac{dx}{dt} y=dtdx (1.3)
对上式取微分得： d y d t = d 2 x d t 2 \frac{dy}{dt}=\frac{d^2x}{dt^2} dtdy=dt2d2x (1.4)
将式（1.3）和（1.4）代入式（1.2）中得到： m d y d t + c y + k x = 0 m\frac{dy}{dt} +cy + kx = 0 mdtdy+cy+kx=0 (1.5)
于是原来的二阶微分方程（1.2）可等价于两个一阶方程组（1.3）和（1.5）。
同样的，其他的n阶微分方程也可以用类似的方式简化。
1.2 常微分方程的非计算机求解方法
常微分方程通常采用解析积分得方法来求解。如对于式(1.1)，先乘以dt，在进行积分得到： v = ∫ ( g − c m v ) d t v = \int{(g-\frac{c}{m}v)dt} v=∫(g−mcv)dt (1.6)
对于上式（1.6），是可以精确的推导出该积分得函数表达式的。因为该方程是线性的。
但在实际中，很多方程（是非线性的）精确解是无法求出的。于是提出了一个方法，就是将方程线性化。 （ n n n阶）线性常微分方程的一般形式是： a n ( x ) y ( n ) + . . . + a 2 ( x ) y ( 2 ) + a 1 ( x ) y ′ + + a 0 ( x ) y = f ( x ) a_n(x)y^{(n)}+...+a_2(x)y^{(2)}+a_1(x)y'++a_0(x)y = f(x) an(x)y(n)+...+a2(x)y(2)+a1(x)y′++a0(x)y=f(x) (1.7) 其中， y ( n ) y^{(n)} y(n)是y关于x的n阶导数， a n ( x ) a_n(x) an(x)和 f ( x ) f(x) f(x)都是关于x的函数。因为该方程中未出现因变量y与其导数的乘积，也没有出现非线性函数。所以认为它是线性的。
如下式（1.8）是一个非线性微分方程： d 2 x d t 2 + g l s i n ( x ) = 0 \frac{d^2x}{dt^2} +\frac{g}{l} sin(x)= 0 dt2d2x+lgsin(x)=0 (1.8) 由于含有 s i n ( x ) sin(x) sin(x)为非线性函数，故该微分方程是非线性的。
线性常微分方程是可以通过解析法求解的。但是，大部分非线性方程无法精确求解。
1.3 线性微分方程求解的推导过程
拿一个简单的方程举例。首先给定函数： y = − 0.5 x 4 + 4 x 3 − 10 x 2 + 8.5 x + 1 y=-0.5x^4+4x^3-10x^2+8.5x+1 y=−0.5x4+4x3−10x2+8.5x+1 （1.9） 这是一个四次多项式。对其进行微分，就得到一个常微分方程： d y d x = − 2 x 3 + 12 x 2 − 20 x + 8.5 \frac{dy}{dx}=-2x^3+12x^2-20x+8.5 dxdy=−2x3+12x2−20x+8.5 （1.10）
对式（1.10）乘以dx，在进行积分得到： y = ∫ ( − 2 x 3 + 12 x 2 − 20 x + 8.5 ) d x y=\int{(-2x^3+12x^2-20x+8.5)}dx y=∫(−2x3+12x2−20x+8.5)dx （1.11） 应用积分法则得出解为： y = − 0.5 x 4 + 4 x 3 − 10 x 2 + 8.5 x + C y=-0.5x^4+4x^3-10x^2+8.5x+C y=−0.5x4+4x3−10x2+8.5x+C （1.12） 除了相差一个C外，其余都与原函数相同。这个C称为积分常数（constant of integration）。 出现一个任意常数C表明，积分的结果并不算是唯一的。无限多个常数C对应无限多个可能的函数，都满足微分方程。下图给出了6个满足条件的函数：
为了将解完全确定下来，微分方程通常伴随有辅助条件（auxiliary conditions）。对于一阶常微分方程，有一类被称为**初值（initial value）**的辅助条件，这类条件用于确定常数值，从而使得解是唯一的。例如，给式（1.11）添加初始条件x=0，y=1。带入式（1.12）中，可推导出C=1。于是就得到了唯一解： y = − 0.5 x 4 + 4 x 3 − 10 x 2 + 8.5 x + 1 y=-0.5x^4+4x^3-10x^2+8.5x+1 y=−0.5x4+4x3−10x2+8.5x+1 这个解同时满足常微分方程和指定的初始条件。
当处理n阶微分方程时，就需要n个条件来确定唯一解。如果所有的条件都是在自变量同一值处指定的，那么问题就称为初值问题（initial-value problem）。与之相对的，边值问题（boundary-value problems），就是指在自变量的不同值处指定初始条件。
2 一阶常微分方程（ODE）求解方法
2.1 欧拉方法
2.1.1 欧拉方法
还是拿式（1.1）举例：
d v d t ≈ Δ v Δ t = v ( t i + 1 ) − v ( t i ) t i + 1 − t i \frac{dv}{dt} \approx \frac{\Delta v}{\Delta t} = \frac{v(t_{i+1}) - v(t_{i})}{t_{i+1} - t_{i}} dtdv≈ΔtΔv=ti+1−tiv(ti+1)−v(ti) （2.1） 其中， Δ v \Delta v Δv和 Δ t \Delta t Δt分别为速度与时间的差分， v ( t i ) v(t_{i}) v(ti)为初始时刻 t i t_i ti的速度， v ( t i + 1 ) v(t_{i+1}) v(ti+1)为下一个时刻 t i + 1 t_{i+1} ti+1的速度。注意， d v d t ≈ Δ v Δ t \frac{dv}{dt} \approx \frac{\Delta v}{\Delta t} dtdv≈ΔtΔv 是一个近似计算，因为： d v d t = lim Δ t → 0 Δ v Δ t \frac{dv}{dt} = \lim_{\Delta t \to 0}\frac{\Delta v}{\Delta t} dtdv=limΔt→0ΔtΔv
式(2.1)称为在时刻 t i t_i ti处导数的有限差商(finite divided difference)逼近。将其代入式（1.1）中可得： d v d t ≈ v ( t i + 1 ) − v ( t i ) t i + 1 − t i = g − c m v ( t i ) \frac{dv}{dt} \approx\frac{v(t_{i+1}) - v(t_{i})}{t_{i+1} - t_{i}} = g-\frac{c}{m}v(t_i) dtdv≈ti+1−tiv(ti+1)−v(ti)=g−mcv(ti) （2.2）
对该方程进行整理可得: v ( t i + 1 ) = v ( t i ) + [ g − c m v ( t i ) ] ( t i + 1 − t i ) v(t_{i+1}) = v(t_{i}) + [g-\frac{c}{m}v(t_i)](t_{i+1} - t_{i}) v(ti+1)=v(ti)+[g−mcv(ti)](ti+1−ti) （2.3）
注意1：在式（2.3）的中括号内部分 g − c m v ( t i ) g-\frac{c}{m}v(t_i) g−mcv(ti)其实就是式（1.1）常微分方程的右边 g − c m v g-\frac{c}{m}v g−mcv，是不过这个是一个估计值，是拿 t i t_i ti对刻的速度来估计 t i t_i ti ~ t i + 1 t_{i+1} ti+1这段时间的速度，但只要两个时间点距离越近，这个估值就准确。
根据式（2.3），如果给定了在某一对刻 t i t_i t