imtoken钱包下载官方|ode

通俗易懂的微分方程理论：从ODE到PDE - 知乎

通俗易懂的微分方程理论：从ODE到PDE - 知乎切换模式写文章登录/注册通俗易懂的微分方程理论：从ODE到PDE图通道本文介绍微分方程的基本概念，并以弹簧振子和杆振动两个MATLAB数值算例阐释常微分方程 (ODE) 和 偏微分方程 (PDE) 之间的联系。微分方程（ODE & PDE）方程的形式有很多。代数方程是指含有未知数的等式；而微分方程，是指含有未知函数及其导数的等式。代数方程的解是常数，而微分方程的解是函数。现有未知数x, t,一元未知函数f和二元未知函数u，分别构造出关于x的代数方程，关于f的常微分方程和关于u的偏微分方程如下：当微分方程的未知函数是一元函数时，列出的方程为常微分方程。此处的“常”（Ordinary）表示平常，也就是在一般情况（理想情况）下的微分方程。上图列出的微分方程中f导数的最高阶次为2，微分方程即为2阶。当讨论微分方程的阶或者齐次性时，不关注未知函数的自变量，而只关注未知函数本身。但现实世界是复杂的，如电磁理论和电动力学（麦克斯韦方程组），流体运动（欧拉方程组和纳维-斯托克斯方程组），弹性体变形（平衡微分方程组），量子力学基本规律（薛定锷方程），等等。因此未知函数通常是多元函数（如与时间、位置相关），此时列出的方程为偏微分方程。从定义中，可以看出，常微分方程是偏微分方程的子集。常微分方程（弹簧振子）牛顿第二定律利用微积分，振子的加速度可以表达为位置矢量的二阶导数。即可以把牛顿第二定律表达为如下形式:胡克定律胡克定律表达了振子离开平衡位置的位移与所受弹力成正比。弹力方向始终与位移的方向相反（前提是我们把振子的平衡位置定义为原点，即位移为0的位置）。建立常微分方程常微分方程组（弹簧振子串联）以上的常微分方程中只有一个未知函数u，并且该未知函数只与时间t有关。当有数个弹簧振子串联时，未知函数的个数增加，需用常微分方程组才能描述运动。建模方法与单个振子相同，只是胡克定律部分的表达式稍有改变。之所以需要微分方程组，是因为单个常微分方程不能包含振子的位置信息，因此需要用变量名区分！上述常微分方程组中的每一个方程都具有相似的形式：右端项含有位移的二阶差分。此时可以用微分代替差分，自然地将位置信息引入未知函数u，得到可以描述任意多振子的偏微分方程（无数常微分方程的叠加）。偏微分方程（杆振动）根据上述的推演，我们已经可以得到描述无数个振子串联的运动方程，即可以看成杆的纵向振动方程。这一推导过程可能与常见的力学推导有所不同，读者可以与利用应力-应变关系推导的结果进行比较，结果应该是一致的。弦的横向振动具有相似的运动微分方程。MATLAB仿真本采用有限差分法求解杆的纵/横向振动（一维波动方程数值解）。对于离散系统动画，降低空间离散化数目即可。​% 一维波动方程数值解

% 杆的纵向振动

% 有限差分法

% 编写环境: MATLAB R2019b

function rodvibration_4

close all

c = 3; % 波速

xmin = 0;

xmax = 30;

Nx = 200; % 空间离散化数目

tmin = 0;

tmax = 100*10;

Nt = 2401*10; % 时间离散化数目

bc = 'd'; % 边界条件

x = linspace(xmin, xmax, Nx)';

t = linspace(tmin, tmax, Nt);

dx = (xmax - xmin)/(Nx - 1);

dt = (tmax - tmin)/(Nt - 1);

C = c*dt/dx; C2 = C*C;

u = zeros(Nx, Nt);

u(:,1) = u0(x); % 设置初始波形

u(:,2) = u0(x-c*dt);

D2 = @(v) [0; v(1:end-2) - 2*v(2:end-1) + v(3:end); 0]; % 二阶差分（边界元设为 0）

for n = 2:Nt-1

u(:,n+1) = 2*u(:,n) - u(:,n-1) + C2*D2(u(:,n));

u = bc_set(u, n+1, bc, c, dx, dt);

......仿真结果如下：完整代码下载：编辑于 2022-10-29 00:34微分方程Matlab振动​赞同 69​​3 条评论​分享​喜欢​收藏​申请

ODE笔记：简单方程的求解 - 知乎

ODE笔记：简单方程的求解 - 知乎切换模式写文章登录/注册ODE笔记：简单方程的求解LingoZero一个沉迷于低级趣味的人本文为常微分方程课程笔记，主要内容为：利用级数法研究常系数线性齐次ODE，并得到其求解方法利用Picard序列导出常数变易法，进而得到线性非齐次ODE的求解方法1 常系数线性齐次常微分方程的求解1.1 级数法在ODE基本定理中我们证明了：如下的线性ODE系统\begin{cases} y^{(n)}+a_{n-1}y^{(n-1)}+\cdots+a_0(t)y+b(t)=0\\ y^{(k)}(t)\mid_{t=t_0}=b_k&k=0,1,2,\cdots,n-1 \end{cases}\\ 当 a_k(t) , b(t) 均解析时，方程存在唯一的解，且解是解析的.由于解的解析性，我们可以将解表示为 y(t)=\sum_{i=0}^\infty c_k\frac{(t-t_0)^k}{k!} ,级数法的想法正是将解表示为级数形式，并代入原方程求解对应项系数，进而得到解的级数表达式.【例】作为级数法求解ODE的例子，我们讨论下面的方程y'''=y\\ 设方程的解为 y(t)=\sum_{i=0}^\infty c_k\frac{(t-t_0)^k}{k!}，代入方程有\begin{aligned} y'''&=\sum_{i=3}^\infty c_kk(k-1)(k-2)\frac{(t-t_0)^{k-3}}{k!}=\sum_{i=0}^\infty c_{k+3}\frac{(t-t_0)^k}{k!}\\ &=y=\sum_{i=0}^\infty c_{k}\frac{(t-t_0)^k}{k!} \end{aligned}\\ 即数列满足 c_{k+3}=c_k ，方程的解可以写为\begin{aligned} y(t)&=c_0\sum_{i=0}^\infty \frac{(t-t_0)^{3k}}{(3k)!}+c_1\sum_{i=0}^\infty \frac{(t-t_0)^{3k+1}}{(3k+1)!}+c_2\sum_{i=0}^\infty \frac{(t-t_0)^{3k+2}}{(3k+2)!}\\ &=c_0\varphi_0 (t)+c_1\varphi_1 (t)+c_2\varphi_2 (t) \end{aligned} 考虑\begin{cases} e^{t-t_0}=\sum_{i=0}^\infty \frac{(t-t_0)^k}{k!}\\ e^{\rho(t-t_0)}=\sum_{i=0}^\infty \frac{\rho^k(t-t_0)^k}{k!}\\ e^{\rho^2(t-t_0)}=\sum_{i=0}^\infty \frac{\rho^{2k}(t-t_0)^k}{k!}\\ \end{cases}\\ 其中 \rho=e^{\frac{2\pi i}{3}} (i=0,1,2) （这里引入的目的在1.2中将得到说明）\frac13 (e^{t-t_0}+e^{\rho(t-t_0)}+e^{\rho^2(t-t_0)})=\sum_{i=0}^\infty \frac{1+\rho+\rho^2}3\frac{(t-t_0)^k}{k!}=\sum_{i=0}^\infty \frac{(t-t_0)^{3k}}{(3k)!}=\varphi_0(t) 对上式两边分别求导可依次求得 \varphi_1(t) 和 \varphi_2(t) ，代入即为原方程的解为.【注】：严格来说到这里并没有结束. 还需要说明 \varphi_0(t)，\varphi_1(t) , \varphi_2(t) 构成的解系是线性无关的（级数法的并不总是保证这一点），需要通过朗斯基行列式进行判定.1.2 常系数线性齐次常微分方程简便起见，下面以二阶的情况为例讨论，更高阶的情况是类似的.y''+by'+c=0\\ 设解的形式为 y(t)=\sum_{i=0}^\infty c_k\frac{(t-t_0)^k}{k!}，考虑到y'(t)=\sum_{i=0}^\infty c_{k+1}\frac{(t-t_0)^k}{k!}\\ y'(t)=\sum_{i=0}^\infty c_{k+1}\frac{(t-t_0)^k}{k!} 代入原方程可化简为c_{k+2}+bc_{k+1}+cc_k=0\tag{*}\\ 则问题转化为求解数列 \{c_k\} .(*) 式等价于\begin{bmatrix}c_{k+2}\\c_{k+1}\end{bmatrix}= \begin{bmatrix}-b&-c\\1&0\end{bmatrix} \begin{bmatrix}c_{k+1}\\c_{k}\end{bmatrix}\\ 利用线性代数的相关知识可以求解得：令 \det(r I-A) =0 ，设求出的特征根分别为 r_1 , r_2 （即方程 r^2+br+c=0 的解）（1）当 r_1\not=r_2 时\begin{aligned} c_k&=r_1^{k}\frac{r_2c_0-c_1}{r_2-r_1}+r_2^k\frac{c_1-r_1c_0}{r_2-r_1}\\ y(t)&=\sum_{i=0}^\infty c_k\frac{(t-t_0)^k}{k!}=\frac{r_2c_0-c_1}{r_2-r_1}e^{r_1(t-t_0)}+\frac{c_1-r_1c_0}{r_2-r_1}e^{r_2(t-t_0)} \end{aligned}\\ 如果不考虑参数的具体形式，则解的结构为y(t)=C_1e^{r_1(t-t_0)}+C_2e^{r_2(t-t_0)}\\ （2）当 r_1=r_2=r 时\begin{aligned} c_k&=r^{k-1}kc_1-r^k(k-1)c_0\\ y(t)&=c_0e^{r(t-t_0)}+(c_1-rc_0)(t-t_0)e^{r(t-t_0)} \end{aligned}\\ 类似上面，如果只考虑解的结构，则有y=C_1e^{r(t-t_0)}+C_2te^{r(t-t_0)}\\ 综合上面的讨论有，二阶常系数线性齐次ODE的解为y(t)=\begin{cases} C_1e^{r_1(t-t_0)}+C_2e^{r_2(t-t_0)}&r_1\not=r_2\\ C_1e^{r(t-t_0)}+C_2(t-t_0)e^{r(t-t_0)}&r_1=r_2=r \end{cases}\\ 上面的结论还可以推广到n阶的情况.2 线性非齐次方程的求解2.1 常数变易法常数变易法通过将线性非齐次ODE转化为线性齐次ODE，进而降低求解难度，下面以一阶线性ODE系统为例进行说明：\begin{cases} \frac{{\rm d}y}{{\rm d}t}=a(t)y+g(t)\\ y\mid_{t=t_0}=y_0 \end{cases}\tag{*}\\ 设方程的解具有下面的形式 y(t)=e^{\int a(t){\rm d}t}C(t)\\ 代入原方程并化简可得\frac{{\rm d}C(t)}{{\rm d}t}=e^{-{\int a(t){\rm d}t}}g(t)\\ 这是一个线性齐次ODE，求解该方程并代入有y(t)=e^{\int a(t){\rm d}t}(\int g(t)e^{-\int a{\rm d}t}{\rm d}t+K)\tag{**}\\ 其中 K 为由初值确定的常数. (**) 即为原方程的解.2.2 常数变易法的原理有很多种解释常数变易法的思路，这里选取两种介绍. 一种是基于配凑的（思维难度较小，但是其实显得有些生硬），一种是基于对Picard序列的研究的（与前序知识的连续性较好）.2.2.1 方法1对一阶线性ODE系统\begin{cases} \frac{{\rm d}y}{{\rm d}t}=a(t)y+g(t)\\ y(t)\mid_{t=t_0}=y_0 \end{cases}\tag{*}\\ 我们假设解具有如下形式y(t)=u(t)v(t)\\ 代入 (*) 中的ODE整理可得u\frac{{\rm d}v}{{\rm d}t}+v(\frac{{\rm d}u}{{\rm d}t}-au)=g(t)\\ 考虑到这里的 u(t) 只需要找到一个特定的函数即可，为了简化方程，令\frac{{\rm d}u}{{\rm d}t}-au=0\\ 则方程化为\begin{cases} u\frac{{\rm d}v}{{\rm d}t}=g(t)\\ \frac{{\rm d}u}{{\rm d}t}-au=0\\ \end{cases}\\ 求解这个方程有y(t)=u(t)v(t)=e^{\int a(t){\rm d}t}\cdot(\int e^{-{\int a(t){\rm d}t}}g(t){\rm d}t+C)\\ 2.2.2 方法2构造 (*) 系统的Picard序列\begin{cases} y_0=y_0\\ y_1=y_0+y_0\int_{t_0}^ta(s)ds+\int_{t_0}^tg(s)ds\\ y_2=y_0+y_0\int_{t_0}^ta(s)ds+y_0\int_{t_0}^t(a(s)\int_{t_0}^sa(s_1)ds_1)ds+\int_{t_0}^tg(s)ds+\int_{t_0}^t(a(s)\int_{t_0}^sg(s_1)ds_1)ds\\\cdots \end{cases}\\ y_k 可以分为两部分，即\begin{aligned} y_k&=y_0(1+\int_{t_0}^ta(s)ds+\int_{t_0}^t(a(s)\int_{t_0}^sa(s_1)ds_1)ds+\cdots)\\ &+(\int_{t_0}^tg(s)ds+\int_{t_0}^t(a(s)\int_{t_0}^sg(s_1)ds_1)ds+\cdots)\end{aligned}\\ 一方面，考虑下面的ODE系统（即原方程的齐次形式）\begin{cases} \frac{{\rm d}y}{{\rm d}t}=a(t)y\\ y(t)\mid_{t=t_0}=y_0 \end{cases}\\ 这个系统对应的Picard序列为 y_{hk}=y_0(1+\int_{t_0}^ta(s)ds+\int_{t_0}^t(a(s)\int_{t_0}^sa(s_1)ds_1)ds+\cdots)\\ 另一方面，容易验证系统 (*) 有特解y_p=\int_{t_0}^tg(s)ds+\int_{t_0}^t(a(s)\int_{t_0}^sg(s_1)ds_1)ds+\cdots\\ (严格来说，这里还需要证明上面两个式子右边的级数是收敛的)则 (*) 的解可以表示为y(t)=y_{hk}(t)+y_p(t)\\ 上面的式子其实利用代数知识也可以直接得到，但是真正有意思的是下面对 y_p 的讨论将 g(x) 设为级数形式 g(x)=\sum_{i=0}^\infty b_k\frac{(t-t_0)^k}{k!} ，则有\begin{aligned} y_p&=\int_{t_0}^tg(s)ds+\int_{t_0}^t(a(s)\int_{t_0}^sg(s_1)ds_1)ds+\cdots\\ &=\int_{t_0}^t\sum_{i=0}^\infty b_k\frac{(s-t_0)^k}{k!}ds+\int_{t_0}^t(a(s)\int_{t_0}^s\sum_{i=0}^\infty b_k\frac{(s_1-t_0)^k}{k!}ds_1)ds+\cdots\\ &= b_0((t-t_0)+\int_{t_0}^ta(s)(t-t_0)ds+\cdots)\\ &+b_1(\frac{(t-t_0)^2}{2!}+\int_{t_0}^ta(s)\frac{(s-t_0)^2}{2!}ds+\cdots)\\ &+\cdots \end{aligned}\\ 考虑到y_{hk}=y_0e^{\int_{t_0}^ta(s)ds}=y_0\sum_{i=0}^\infty\frac{(\int_{t_0}^ta(s)ds)^k}{k!}\\ 在不考虑具体函数形式的情况下，y_p总能化为下面的结构y_p(t)=h(t)y_{hk}(t)\\ 则 (*) 的解可以表示为\begin{aligned} y(t)&=y_{hk}(t)+y_p(t)\\ &=y_{hk}+h(t)y_{hk}\\ &=y_{hk}(1+h(t)) \end{aligned}\\ 如果记 C(t)=1+h(t) ，那么有y(t)=C(t)y_{hk}(t)\\ 则问题转化为求解 (*) 对应的齐次方程的解和一个特定函数 C(t) ，也就是常数变易法的思想. 编辑于 2020-12-18 12:44常微分方程高等数学 (大学课程)​赞同 18​​添加评论​分享​喜欢​收藏​申请

百度知道 - 信息提示

百度知道 - 信息提示

百度首页

商城

注册

登录

网页

资讯

视频

图片

知道

文库

贴吧采购

地图更多

搜索答案

我要提问

百度知道>提示信息

知道宝贝找不到问题了>_

该问题可能已经失效。返回首页

15秒以后自动返回

帮助

 | 意见反馈

 | 投诉举报

京ICP证030173号-1   京网文【2023】1034-029号     ©2024Baidu  使用百度前必读 | 知道协议 

常微分方程（ODE）求解方法总结-CSDN博客

常微分方程（ODE）求解方法总结-CSDN博客

常微分方程（ODE）求解方法总结

一只特立独行的蚂蚁

已于 2023-06-27 20:38:14 修改

阅读量1w

收藏

23

点赞数

6

分类专栏：

数值分析算法

文章标签：

算法

于 2023-06-03 21:50:49 首次发布

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

本文链接：https://blog.csdn.net/luolei188/article/details/131025102

版权

数值分析算法

专栏收录该内容

10 篇文章

17 订阅

订阅专栏

常微分（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−mc​v (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−mc​v)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​+lg​sin(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​−ti​v(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​−ti​v(ti+1​)−v(ti​)​=g−mc​v(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−mc​v(ti​)](ti+1​−ti​) （2.3）

注意1：在式（2.3）的中括号内部分

g

−

c

m

v

(

t

i

)

g-\frac{c}{m}v(t_i)

g−mc​v(ti​)其实就是式（1.1）常微分方程的右边

g

−

c

m

v

g-\frac{c}{m}v

g−mc​v，是不过这个是一个估计值，是拿

t

i

t_i

ti​对刻的速度来估计

t

i

t_i

ti​ ~

t

i

+

1

t_{i+1}

ti+1​这段时间的速度，但只要两个时间点距离越近，这个估值就准确。

根据式（2.3），如果给定了在某一对刻

t

i

t_i

ti​的初始值速度

v

i

v_i

vi​，就很容易计算出下一时刻

t

i

+

1

t_{i+1}

ti+1​的速度

v

i

+

1

v_{i+1}

vi+1​。接着又可以用

t

i

+

1

t_{i+1}

ti+1​时刻的新速度值来计算

t

i

+

2

t_{i+2}

ti+2​的速度，然后依次继续下去。这样，用这种方式就可以计算出任意时刻的速度。

新值

=

旧值

+

斜率

×

步长

新值 = 旧值 + 斜率 \times 步长

新值=旧值+斜率×步长 用数学语言表述为:

v

i

+

1

=

v

i

+

ϕ

×

h

v_{i+1} = v_{i} + \phi \times h

vi+1​=vi​+ϕ×h（2.4）

这种方法的就命名为欧拉法(Euler’s method)(或欧拉一柯西法(Euler-Cauchy method)，或折线法 (point-slope method)。 欧拉方法利用斜率（等于一阶导数在

x

i

x_i

xi​处的值)，通过线性外推的办法预测出y在前进步长h之后

x

i

+

1

x_{i+1}

xi+1​处的新值： 利用上面的方法，取步长为2（也就是每隔2秒，计算一次速度值），于是得到一系列时间t与速度v的值，将这些计算结果与精确解同时画在一个图中： 图5 速度v与时间t的关系的常微分方程的数值解与解析解的对比

上图中，如果我们将步长减小一半，以1秒的步长计算得到的结果，其误差将会更小，从而直线段的轨迹更趋近于真实解。但是步长每减小一半，可以得到更精确的结果，所需要的计算量也会翻一番。

2.1.2 欧拉方法的误差分析

常微分方程数值解的误差包括两种类型： (1)截断(Truncation)误差或离散误差，是由逼近Y值的方法本身引起的。 (2)舍入(Round-off误差，是由于计算机能保留的有效数字有限而引起的。 截断误差由两部分组成：第一部分是局部截断误差(local truncation error)，指方法在一步迭代中所产生的误差。第二部分是传播截断误差(propagated truncation error)，它是由前面各步中产生的逼近值引起的，这两部分的和称为整体截断误差或全局截断误差(global truncation error)。

利用泰勒级数展开直接推导欧拉方法，可以同时推出截断误差的大小和属性。假设被积分的微分方程具有一般形式：

y

′

=

f

(

x

,

y

)

y' = f(x, y)

y′=f(x,y)（2.5） 其中，

y

′

=

d

y

/

d

x

y' = dy/dx

y′=dy/dx，x和y分别为自变量和因变量。

如果解(也就是描述y行为的函数)具有连续导数，那么它可以在点(xi, yi)处展开成泰勒级数，即

y

i

+

1

=

y

i

+

y

i

′

∗

h

+

y

i

′

′

2

!

∗

h

2

+

.

.

.

+

y

i

(

n

)

n

!

∗

h

n

+

R

n

y_{i+1}=y_i+y'_i*h+\frac{y''_i}{2!}*h^2+...+\frac{y^{(n)}_i}{n!}*h^n + R_n

yi+1​=yi​+yi′​∗h+2!yi′′​​∗h2+...+n!yi(n)​​∗hn+Rn​（2.6） 其中，

h

=

(

x

i

+

1

−

x

i

)

h = (x_{i+1}-x_i)

h=(xi+1​−xi​)，

R

n

R_n

Rn​为余项，定义为：

R

n

=

y

(

n

+

1

)

(

ξ

)

(

n

+

1

)

!

∗

(

h

)

(

n

+

1

)

R_n = \frac{y^{(n+1)}(\xi)}{(n+1)!}*(h)^{(n+1)}

Rn​=(n+1)!y(n+1)(ξ)​∗(h)(n+1)（2.7） 其中，

ξ

∈

[

x

i

,

x

i

+

1

]

\xi \in [x_{i}, x_{i+1}]

ξ∈[xi​,xi+1​]。

由于

y

′

=

f

(

x

,

y

)

y' = f(x, y)

y′=f(x,y)，所以

y

′

′

=

f

′

(

x

,

y

)

y'' = f'(x, y)

y′′=f′(x,y)，

y

′

′

′

=

f

′

′

(

x

,

y

)

y''' = f''(x, y)

y′′′=f′′(x,y)，依此类推，代入式（2.6）右边，同时也加入式（2.7），得到：

y

i

+

1

=

y

i

+

f

(

x

i

,

y

i

)

∗

h

+

f

′

(

x

i

,

y

i

)

2

!

∗

h

2

+

.

.

.

+

f

(

n

−

1

)

(

x

i

,

y

i

)

n

!

∗

h

n

+

O

(

h

n

+

1

)

y_{i+1}=y_i+ f(x_i, y_i)*h+\frac{ f'(x_i, y_i)}{2!}*h^2+...+\frac{f^{(n-1)}(x_i, y_i)}{n!}*h^n + O(h^{n+1})

yi+1​=yi​+f(xi​,yi​)∗h+2!f′(xi​,yi​)​∗h2+...+n!f(n−1)(xi​,yi​)​∗hn+O(hn+1)（2.8） 此处

O

(

h

n

+

1

)

O(h^{n+1})

O(hn+1)表示局部截断误差与步长的n+1次幂成正比。

因为我们仅使用了泰勒级数的有限项来逼近真解，所以就产生了截断误差。这样，我们就截断或者略去了一部分真实解。例如，在前面的计算中，我们采用

y

i

+

1

=

y

i

+

f

(

x

i

,

y

i

)

×

h

y_{i+1} = y_{i} + f(x_i, y_i) \times h

yi+1​=yi​+f(xi​,yi​)×h来计算下一个点的值，有部分泰勒级数展开项没有被包括在式中，这些剩余项就导致了欧拉方法的截断误差：

E

t

=

f

′

(

x

i

,

y

i

)

2

!

∗

h

2

+

.

.

.

+

f

(

n

−

1

)

(

x

i

,

y

i

)

n

!

∗

(

h

)

n

+

O

(

h

n

+

1

)

E_t = \frac{ f'(x_i, y_i)}{2!}*h^2+...+\frac{f^{(n-1)}(x_i, y_i)}{n!}*(h)^n + O(h^{n+1})

Et​=2!f′(xi​,yi​)​∗h2+...+n!f(n−1)(xi​,yi​)​∗(h)n+O(hn+1)（2.9）

E

t

E_t

Et​为真实的局部截断误差。

当h足够小时，式(2.9)中各项的误差通常随着阶数的增加而减少。于是常常可以表示为：

E

a

=

f

′

(

x

i

,

y

i

)

2

!

∗

h

2

E_a = \frac{ f'(x_i, y_i)}{2!}*h^2

Ea​=2!f′(xi​,yi​)​∗h2（2.10） 或者：

E

a

=

O

(

h

2

)

E_a =O(h^2)

Ea​=O(h2)（2.11）

E

t

E_t

Et​为近似的局部截断误差。

泰勒级数给出了一种量化欧拉方程误差的方法。不过，这种方法使用时存在局限: (1)泰勒级数仅能估计局部截断误差一即方法在一步迭代中所产生的误差。它不能估量传播截断误差，因此也不能估量全局截断误差。 (2)实际问题中处理的函数一般比简单多项式复杂。因此，计算泰勃级数展开时，导数值有时很难得到。

尽管这些限制使泰勒级数无法分析人多数实际问题的精确误差，但是它仍然有助于了解欧拉方法的性态。按照式(2.11)可发现局部误差与步长的平方和微分方程的一阶导数成正比。此外已经证明，全局截断误差为O(h)，全局截断误差与步长成正比(Carnahan等人，1969)。

结论: (1)缩短步长可以减少误差。 (2)因为直线的二阶导数为常数，所以如果基本函数（即微分方程的解)是线性的，那么方法能给出无误差的预测值。

由于欧拉方法使用直线段来逼近解，因此，欧拉方法称为一阶方法(first-order method)。 思考：如果使用二次来逼近解，也就是二阶方法，就能给出基本函数（即微分方程的解)是二次的问题的精确。同理，如果基本函数是n次多项式，则使用n阶方法可得出精确解。

2.1.3 欧拉方法的改进思路1——添加高阶项

减小欧拉方法的误差的一种方法是，在公式中加入泰勒级数展开的高阶项。例如，添加二阶项后，得到：

y

i

+

1

=

y

i

+

f

(

x

i

,

y

i

)

∗

h

+

f

′

(

x

i

,

y

i

)

2

!

∗

h

2

y_{i+1}=y_i+ f(x_i, y_i)*h+\frac{ f'(x_i, y_i)}{2!}*h^2

yi+1​=yi​+f(xi​,yi​)∗h+2!f′(xi​,yi​)​∗h2（2.12） 于是局部截断误差变为：

E

a

=

f

′

′

(

x

i

,

y

i

)

3

!

∗

h

3

E_a = \frac{ f''(x_i, y_i)}{3!}*h^3

Ea​=3!f′′(xi​,yi​)​∗h3（2.13）

尽管添加高阶项对于处理多项式来说十分简单，不过有时候会使常微分方程变得更加复杂。特别地，如果常微分方程同时为自变量和因变量的函数，就需要通过链式法则来计算微分。例如，f(x,y)的一阶导数是：

f

′

(

x

,

y

)

=

∂

f

(

x

i

,

y

i

)

∂

x

+

∂

f

(

x

i

,

y

i

)

∂

y

d

y

d

x

f'(x,y) = \frac{ \partial{f(x_i, y_i)}}{\partial{x}} + \frac{ \partial{f(x_i, y_i)}}{\partial{y}}\frac{dy}{dx}

f′(x,y)=∂x∂f(xi​,yi​)​+∂y∂f(xi​,yi​)​dxdy​ 二阶导数是：

f

′

′

(

x

,

y

)

=

∂

[

∂

f

/

∂

x

+

(

∂

f

/

∂

y

)

(

d

y

/

d

x

)

]

∂

x

+

∂

[

∂

f

/

∂

x

+

(

∂

f

/

∂

y

)

(

d

y

/

d

x

)

]

∂

y

d

y

d

x

f''(x,y) = \frac{ \partial[\partial{f}/\partial{x} + (\partial{f}/\partial{y})(dy/dx) ]}{\partial{x}} + \frac{ \partial[\partial{f}/\partial{x} + (\partial{f}/\partial{y})(dy/dx) ]}{\partial{y}}\frac{dy}{dx}

f′′(x,y)=∂x∂[∂f/∂x+(∂f/∂y)(dy/dx)]​+∂y∂[∂f/∂x+(∂f/∂y)(dy/dx)]​dxdy​

高阶导数会变得越来越复杂。而综合其计算量，其在性能上其实与一阶泰勒级数方法差不多，而一阶泰勒展开的欧拉方法只需要计算一阶导数，实现起来会简洁很多。因此，添加高阶项的方法其实并不是很实用。

2.1.4 欧拉方法的改进思路2——对斜率采用更好的估值方法

在“注意1”中我们提到，欧拉方法是将区间左端点处的斜率当作斜率在整个区间上平均值的逼近。但估值毕竟是估值，不是精确值，所以这个估值的精确度会严重影响计算的结果。 于是我们思考，有没有什么其他估计斜率的其他方法，能够给出更好的估值，从而得到更精确的预测值。 基于这个想法，下面给出了两种简单的修正方法来克服这个估值不精确的问题：修恩法和中点方法。 (这两种方法其实都属于一个更广泛的类型—龙格一库塔法。)

2.1.4.1 修恩法

改进斜率估计值的一种方法是，在区间上计算两个导数值：一个在区间的起始点处，一个在区间的终点处。然后，将这两个导数的平均值作为斜率在整个区间上的改进估计值。这种方法，称为修恩法(Heun’s method)，其示意图如下图所示。 图5 修恩法的示意图：(a)预估和(b)校正

前面介绍过，在欧拉方法中，区间起始点处的斜率为：

y

i

′

=

f

(

x

i

,

y

i

)

y_i' = f(x_i, y_i)

yi′​=f(xi​,yi​) (2.15) 利用这个值线性外推出

y

i

+

1

y_{i+1}

yi+1​：

y

i

+

1

0

=

y

i

+

f

(

x

i

,

y

i

)

h

y^0_{i+1} = y_i + f(x_i, y_i)h

yi+10​=yi​+f(xi​,yi​)h (2.16) 标准的欧拉方法算到这一步就停止了。然而，在修恩法中，由式(2.16)得出的

y

i

+

1

0

y^0_{i+1}

yi+10​并不是最终结果，而是一个中间的预测值。式(2.16)称为预估表达式(predictor equation)。该式给出了

y

i

+

1

y_{i+1}

yi+1​的一个估计值，利用这个估计值又可以计算区间终点处的斜率（形式如式(2.15)）：

y

i

+

1

′

=

f

(

x

i

+

1

,

y

i

+

1

0

)

y_{i+1}' = f(x_{i+1}, y^0_{i+1})

yi+1′​=f(xi+1​,yi+10​) (2.17) 我们得到两个斜率估计值

y

i

′

y_{i}'

yi′​和

y

i

+

1

′

y_{i+1}'

yi+1′​(式(2.15)和式(2.17))，于是可以计算斜率的区间平均值：

y

i

ˉ

′

=

y

i

′

+

y

i

+

1

′

2

=

f

(

x

i

,

y

i

)

+

f

(

x

i

+

1

,

y

i

+

1

0

)

2

\bar{y_i}' = \frac{y_i' + y_{i+1}'}{2} = \frac{ f(x_i, y_i) + f(x_{i+1}, y^0_{i+1})}{2}

yi​ˉ​′=2yi′​+yi+1′​​=2f(xi​,yi​)+f(xi+1​,yi+10​)​ (2.18)

然后，按照欧拉方法，利用这个平均斜率由

y

i

y_i

yi​线性外推出

y

i

+

1

y_{i+1}

yi+1​:

y

i

+

1

=

y

i

+

f

(

x

i

,

y

i

)

+

f

(

x

i

+

1

,

y

i

+

1

0

)

2

h

y_{i+1} = y_i + \frac{ f(x_i, y_i) + f(x_{i+1}, y^0_{i+1})}{2} h

yi+1​=yi​+2f(xi​,yi​)+f(xi+1​,yi+10​)​h (2.19) 此式称为校正表达式((corrector equation)。 修恩法是一种预估一校正方法{predictor-corrector approach）。

修恩法可简洁地表示为： 预估：

y

i

+

1

0

=

y

i

+

f

(

x

i

,

y

i

)

h

y^0_{i+1} = y_i + f(x_i, y_i)h

yi+10​=yi​+f(xi​,yi​)h (2.16) 校正：

y

i

+

1

=

y

i

+

f

(

x

i

,

y

i

)

+

f

(

x

i

+

1

,

y

i

+

1

0

)

2

h

y_{i+1} = y_i + \frac{ f(x_i, y_i) + f(x_{i+1}, y^0_{i+1})}{2} h

yi+1​=yi​+2f(xi​,yi​)+f(xi+1​,yi+10​)​h (2.19)

迭代修恩法： 因为式(2.19)的等号两边都含有

y

i

+

1

y_{i+1}

yi+1​，所以该式可以按迭代方式计算，可以反复使用旧的估计值计算出改进后的

y

i

+

1

y_{i+1}

yi+1​值。// 需要注意的是，这个迭代过程并不一定能收敛到真实解，不过它将收敛到一个截断误差有限的估计值。

校正过程收敛的终止条件为：

∣

ε

a

∣

=

∣

y

i

+

1

j

−

y

i

+

1

j

−

1

y

i

+

1

j

∣

100

%

|\varepsilon_a| = |\frac{y^j_{i+1} - y^{j-1}_{i+1}}{y^j_{i+1}}|100\%

∣εa​∣=∣yi+1j​yi+1j​−yi+1j−1​​∣100% 其中，

y

i

+

1

j

−

1

y^{j-1}_{i+1}

yi+1j−1​和

y

i

+

1

j

y^j_{i+1}

yi+1j​分别为上一校正步的结果和当前校正步的结果。

// 前面没有迭代校正的修恩法，可称为简单修恩法；经过迭代计算校正的修恩法，可以称为迭代修恩法。

特殊情况： 在前面的例子中，导数

d

y

/

d

x

dy/dx

dy/dx同时为因变量y和自变量x的函数（

d

y

/

d

x

=

f

(

x

,

y

)

dy/dx = f(x,y)

dy/dx=f(x,y)）。对于常微分方程仅与自变量有关的情况（

d

y

/

d

x

=

f

(

x

)

dy/dx = f(x)

dy/dx=f(x)），求解过程中并不需要预估步（式(2.16)），而且每次迭代中也只需一次校正。在这种情况下，方法简化为

y

i

+

1

=

y

i

+

f

(

x

i

)

+

f

(

x

i

+

1

)

2

h

y_{i+1} = y_i + \frac{ f(x_i) + f(x_{i+1})}{2} h

yi+1​=yi​+2f(xi​)+f(xi+1​)​h (2.19)

该式与“数值积分方法的总结”第1.1.1节 梯形法则（2点积分）其实是等价的。相应的截断误差也与梯形法则一样为：

E

t

=

−

f

′

′

(

ξ

)

12

h

3

E_t = - \frac{f''(\xi)}{12}h^3

Et​=−12f′′(ξ)​h3

由于真解是二次多项式时，常微分方程的二阶导数为0，所以方法是二阶精度的。另外，局部和全局误差分别为

O

(

h

3

)

O(h^3)

O(h3)和

O

(

h

2

)

O(h^2)

O(h2)，因此，减小步长时，误差减小的速度比欧拉方法快。

2.1.4.2 中点方法

图5和式(2.16)（

y

i

+

1

0

=

y

i

+

f

(

x

i

,

y

i

)

h

y^0_{i+1} = y_i + f(x_i, y_i)h

yi+10​=yi​+f(xi​,yi​)h ）使用了端点处的值

y

i

+

1

y_{i+1}

yi+1​ 来修正斜率的。

那么是不是也可以考虑用中点处的值来修正斜率？？

于是产生了另一种简单的欧拉方法的修正方法（如图6所示）。该方法称为中点方法(midpoint method)(或改进的多边形法（improved polygon），或修正的欧拉方法(modified Euler))，利用欧拉方法预测区间中点处的y值（图6(a)）：

y

i

+

1

/

2

=

y

i

+

f

(

x

i

,

y

i

)

h

/

2

y_{i+1/2} = y_i + f(x_i, y_i)h/2

yi+1/2​=yi​+f(xi​,yi​)h/2 (2.20) 然后用这个值来计算中点处的斜率：

y

i

+

1

/

2

′

=

f

(

x

i

+

1

/

2

,

y

i

+

1

/

2

)

y_{i+1/2}' = f(x_{i+1/2}, y_{i+1/2})

yi+1/2′​=f(xi+1/2​,yi+1/2​) (2.21) 假设这个值是对整个区间的斜率平均值的有效近似，那么由

x

i

x_i

xi​线性外推出

x

i

+

1

x_{i+1}

xi+1​处的值为:

y

i

+

1

=

y

i

+

f

(

x

i

+

1

/

2

,

y

i

+

1

/

2

)

h

y_{i+1} = y_i + f(x_{i+1/2}, y_{i+1/2}) h

yi+1​=yi​+f(xi+1/2​,yi+1/2​)h (2.22)

// 注意对比式(2.20)、(2.20)、(2.20)与 (2.16)、 (2.17)、 (2.19)，就可以清晰地看出中点法和修恩法的区别。

图6 中点公式的示意图。(a)式(2.20)}和(b)式(2.22)

因为式（2.22）只有等式的一边含

y

i

+

1

y_{i+1}

yi+1​，所以这种方法的结果无法像修恩法一样采用迭代校正的办法来改进。

（本方法也可同牛顿一柯特斯开型公式的中点公式结合起来理解。）

因为中点方法利用区间中点处的斜率估计值，它理论上比端点（起点）处的斜率估计值要更好(//关于数值微分的讨论知，中心有限差分逼近导数的效果比前向或后向差分都要好)，所以中点方法优于欧拉方法。中点方法的局部和全局误差分别为

O

(

h

3

)

O(h^3)

O(h3)和

O

(

h

2

)

O(h^2)

O(h2)。

2.1.5 欧拉方法的改进思路3

思考：上一节提到的中点法，就类似计算积分的时候，取中点计算估值（1点积分）；修恩法，就类似于梯形法计算积分（2点积分）。回想在“数值积分方法的总结”积分计算中，除了梯形法则（2点积分），我们还讲到了3点积分和4点积分，区别就在于用二次曲线和三次曲线来逼近使积分结果更加接近真实值。那么，这个地方，如果我们能计算或者估算出三个点或者四个点，是否也能用类似的方法，使求解的结果更接近真实值？

2.2 龙格-库塔法

龙格一库塔法(Range-Kutta)(RK)无需计算高阶导数，也可以达到与泰勒级数方法同样的精度。龙格一库塔法存在许多不同版本，但它们都可以写成式(25.1)的一般形式

y

i

+

1

=

y

i

+

ϕ

(

x

i

,

y

i

,

h

)

h

y_{i+1} = y_i + \phi(x_i, y_i, h)h

yi+1​=yi​+ϕ(xi​,yi​,h)h (2.23) 其中，

ϕ

(

x

i

,

y

i

,

h

)

\phi(x_i, y_i, h)

ϕ(xi​,yi​,h)称为增量函数((increment function)，它表示整个区间上的斜率。增量函数的一般形式是 k1出现在k2的表达式中，k2又出现在k3的表达式中，等等。因为每个k都是由函数赋值的，所以这种 递归关系提高了龙格-库塔法的计算效率。

将增量函数的项数记为n，n值的不同就对应着不同类型的龙格一库塔法。n=1时，一阶龙格-库塔法，实际上就是欧拉方法。

2.2.1 二阶龙格-库塔法

前面提到，欧拉方法其实是利用了泰勒的一阶展开。那么，如果对应泰勒二阶展开，就是

y

i

+

1

=

y

i

+

f

(

x

i

,

y

i

)

∗

h

+

f

′

(

x

i

,

y

i

)

2

!

∗

h

2

y_{i+1}=y_i+f(x_i,y_i)*h+\frac{f'(x_i,y_i)}{2!}*h^2

yi+1​=yi​+f(xi​,yi​)∗h+2!f′(xi​,yi​)​∗h2（2.24） 通过链式微分法则计算

f

′

(

x

,

y

)

f'(x,y)

f′(x,y)

f

′

(

x

i

,

y

i

)

=

∂

f

(

x

i

,

y

i

)

∂

x

+

∂

f

(

x

i

,

y

i

)

∂

y

d

y

d

x

f'(x_i,y_i) = \frac{ \partial{f(x_i, y_i)}}{\partial{x}} + \frac{ \partial{f(x_i, y_i)}}{\partial{y}}\frac{dy}{dx}

f′(xi​,yi​)=∂x∂f(xi​,yi​)​+∂y∂f(xi​,yi​)​dxdy​（2.25） 将式(2.25)代入式(2.24)有

y

i

+

1

=

y

i

+

f

(

x

i

,

y

i

)

∗

h

+

[

∂

f

(

x

i

,

y

i

)

∂

x

+

∂

f

(

x

i

,

y

i

)

∂

y

d

y

d

x

]

∗

h

2

2

!

y_{i+1}=y_i+f(x_i,y_i)*h+[\frac{ \partial{f(x_i, y_i)}}{\partial{x}} + \frac{ \partial{f(x_i, y_i)}}{\partial{y}}\frac{dy}{dx}]*\frac{h^2}{2!}

yi+1​=yi​+f(xi​,yi​)∗h+[∂x∂f(xi​,yi​)​+∂y∂f(xi​,yi​)​dxdy​]∗2!h2​（2.26）

而龙格-库塔法中，当k取2时，即为二阶龙格-库塔法：

y

i

+

1

=

y

i

+

(

a

1

k

1

+

a

2

k

2

)

∗

h

y_{i+1}=y_i + (a_1k_1 + a_2k_2)*h

yi+1​=yi​+(a1​k1​+a2​k2​)∗h（2.27）

k

1

=

f

(

x

i

,

y

i

)

k_1=f(x_i,y_i)

k1​=f(xi​,yi​)（2.27a）

k

2

=

f

(

x

i

+

p

1

h

,

y

i

+

q

11

k

1

h

)

k_2=f(x_i + p_1h, y_i + q_{11}k_1h)

k2​=f(xi​+p1​h,yi​+q11​k1​h)（2.27b） 这个地方，

k

2

k_2

k2​写成了

f

(

x

i

+

p

1

h

,

y

i

+

q

11

k

1

h

)

f(x_i + p_1h, y_i + q_{11}k_1h)

f(xi​+p1​h,yi​+q11​k1​h)的形式，这很自然联想到也可以用泰勒展开表达该式，于是得到： 将k1、k2都带入式2.27中得到： 想要让二阶龙格-库塔法（即式2.27）等价于二阶泰勒展开（2.26），那么就必须满足：

a

1

+

a

2

=

1

a_1+a_2 = 1

a1​+a2​=1

a

1

p

1

=

1

/

2

a_1p_1=1/2

a1​p1​=1/2

a

2

q

11

=

1

/

2

a_2q_{11}=1/2

a2​q11​=1/2

由于有三个方程，四个未知数，所以上述方程组有无穷多个解。

其中，几种典型的解的意义： 如果令a2 = 1/2，那么a1=1/2，p1=q11=1，这时其实就是单步的修恩法； 如果令a2 = 1，那么a1=0，p1=q11=1/2，这时其实就是中点方法。 如果令a2 = 2/3，那么a1=1/3，p1=q11=3/4，这种方法叫罗森方法。

值得注意的是：罗森和Rabinowitz指出，对于二阶龙格一库塔法，在所有取值中，当a2取2/3时对应的截断误差界最小。此时：

二阶龙格-库塔法局部和全局误差分别为

O

(

h

3

)

O(h^3)

O(h3)和

O

(

h

2

)

O(h^2)

O(h2)。

2.2.2 三阶龙格-库塔法

对于n=3的情况，公式的推导过程与二阶方法类似。推导的结果是含8个未知数的6个方程。因此，必须事先指定两个未知数的值，才能确定剩下的参数。一类常用的三阶龙格-库塔法是：

y

i

+

1

=

y

i

+

1

6

(

k

1

+

4

k

2

+

k

3

)

∗

h

y_{i+1}=y_i + \frac{1}{6} (k_1 + 4k_2 + k_3)*h

yi+1​=yi​+61​(k1​+4k2​+k3​)∗h（2.28）

如果导数仅为x的函数，那么这个三阶方法退化为辛普森1/3法则。

三阶龙格一库塔法，其局部和全局误差分别为

O

(

h

4

)

O(h^4)

O(h4)和

O

(

h

3

)

O(h^3)

O(h3)。并且能精确地求出三次多项式解。

2.2.2+ 牛顿积分公式

这个地方，有必要将牛顿-科特斯积分公式的系数拿出来和龙格-库塔法的系数做一下对比：

2.2.3 四阶龙格-库塔法

龙格-库塔法中，最流行的是四阶方法。与二阶方法一样，四阶方法也有无限多个版木。下面是最常用到的形式，称为经典四阶龙格-库塔法(classical fourth-order RK method)：

y

i

+

1

=

y

i

+

1

6

(

k

1

+

2

k

2

+

2

k

3

+

k

4

)

∗

h

y_{i+1}=y_i + \frac{1}{6} (k_1 + 2k_2 + 2k_3 + k_4)*h

yi+1​=yi​+61​(k1​+2k2​+2k3​+k4​)∗h（2.29） 其中， 四阶龙格一库塔法，其局部和全局误差分别为

O

(

h

5

)

O(h^5)

O(h5)和

O

(

h

4

)

O(h^4)

O(h4)。并且能精确地求出四次多项式解。

注意：对于前面的二阶、三阶龙格一库塔法，其系数都与牛顿-科特斯积分公式的系数相同。但是对于牛顿积分公式，对应的辛普森3/8法则的系数为1/8、3/8、3/8、1/8，而这个地方为1/6、2/6、2/6、1/6，读者可以下去自己尝试一下如果四阶龙格-库塔法系数用1/8、3/8、3/8、1/8，计算会不会有什么不同？？

2.2.4 高阶龙格-库塔法

对于精度要求更高的情况，建议使用巴特切五阶龙格-库塔法。

y

i

+

1

=

y

i

+

1

90

(

7

k

1

+

32

k

3

+

12

k

4

+

32

k

5

+

7

k

6

)

∗

h

y_{i+1}=y_i + \frac{1}{90} (7k_1 + 32k_3 + 12k_4 + 32k_5+ 7k_6)*h

yi+1​=yi​+901​(7k1​+32k3​+12k4​+32k5​+7k6​)∗h（2.30） 其中， 请注意，巴特切五阶方法中需要计算6次函数值。

尽管存在许多高阶龙格-库塔法，但是，一般地，对于四阶以上公式来说，精度的提高会被计算量和复杂度的增加而抵消。

2.2.5 龙格-库塔法分析

从上面的一到五阶龙格-库塔法公式可以看出，当阶数<=4时，需要计算函数值的次数等于方法的阶数。当n为5或者大于5的时候，计算函数值的次数或大于方法的阶数。所以，四阶龙格-库塔法是一个比较经典的阶数。

如果分别利用欧拉方法、非迭代修恩法、三阶龙格-库塔法、经典四阶龙格-库塔法和巴特切五阶龙格一库塔法进行计算。各方法的百分比相对误差的绝对值与计算开销的关系如图25-16所示。这里，计算开销等于求解过程中需要计算的函数值的个数，即

开销

=

n

f

b

−

a

h

开销 = n_f \frac{b-a}{h}

开销=nf​hb−a​（2.31） 其中，

n

f

n_f

nf​为特定的龙格一库塔法中需计算的函数值个数。(b - a)/h为整个积分区间长度除以步长，也就是说，它表示计算过程中需要应用龙格-库塔方法的次数。这样，因为函数值计算一般占去了单步迭代中的大部分时间，所以式2.31给出了计算时间的一个粗略估计。

从上面的分析中，可以得出一些结论：首先，在同样的计算开销下，高阶方法的结果更精确；其次，在某点之后，随着计算开销的增加，精度的提高开始变慢(注意曲线起初下降得很快，然后趋于平缓)。在应用中应该优先选择高阶龙格一库塔法。//不过，在选择求解方法时，还有许多其他因素需要考虑，比如编程的代价和问题的精度要求。

3 一阶常微分方程组

在本文第1.1节中讲到，高阶微分方程可以分解成多个一阶方程组的形式。

许多实际问题往往也需要求解几个联立的常微分方程组成的方程组。

方程组的一般表示形式为：

d

y

1

d

x

=

f

(

x

,

y

1

,

y

2

,

.

.

.

,

y

n

)

\frac{dy_1}{dx} = f(x, y_1, y_2, ..., y_n)

dxdy1​​=f(x,y1​,y2​,...,yn​)

d

y

2

d

x

=

f

(

x

,

y

1

,

y

2

,

.

.

.

,

y

n

)

\frac{dy_2}{dx} = f(x, y_1, y_2, ..., y_n)

dxdy2​​=f(x,y1​,y2​,...,yn​) …

d

y

n

d

x

=

f

(

x

,

y

1

,

y

2

,

.

.

.

,

y

n

)

\frac{dy_n}{dx} = f(x, y_1, y_2, ..., y_n)

dxdyn​​=f(x,y1​,y2​,...,yn​)（3.1）

要求解这样的一个方程组，必须知道x的起始点处的n个初始条件。

3.1 欧拉方法

上一节中所讲的关于单个方程数值解法的讨论，都可以推广到上面给出的方程组。

每个方程组的求解过程为，在开始下一步迭代之前，对当前步中的每个方程应用单步法求解。 重复上述过程，直到x=2。

3.2 龙格-库塔法

上一节中所讲的任何高阶龙格-库塔法都可以应用于方程组的求解。// 方法与欧拉法类似。

4 自适应龙格·库塔法

到目前为止，我们讨论了步长取常数的常微分方程数值解法。

但实际中，方程的曲线在有的地方变化缓慢，有的地方变化剧烈。在变化缓慢的地方可以用大一点的步长，在变化剧烈的地方可以取小一点的步长。于是就有了自动调整步长的算法。

自动调整步长的算法可以避免保证计算精度的同时提高效率。

自适应步长控制((adaptive step-sixe control)的方法在执行过程中需要估计每一步的局部截断误差。然后，根据该误差估计值确定步长应该变大还是缩小。

4.1 步长控制中的误差估计方法

在单步法中加入自适应步长控制主要有两种方式： 第一种方式是，计算同阶龙格-库塔法在两个不同步长F的预测值，将这两个值之差作为误差的估计值； 第二种方式是，算出两个不同阶龙格-库塔法的预测值，将它们的差作为局部截断误差的估计值。

4.1.1 自适应龙格-库塔法或步长——对分法

步长控制的第一种方式，就是先用步长h计算一个

y

h

y_h

yh​，然后再用步长为h/2分两步计算

y

h

y_h

yh​，比较这两个计算结果的差值，将这个差值作为局部截断误差的估计值。 若记

y

1