Ordinary Differential Equations | XiangSiqi’s Blog

type

Post

status

Published

date

Jan 2, 2026

slug

ODE

summary

1. 龙格-库塔法

平均变化率

对于常微分方程：

将离散化为，这里，倘若已知，我们希望通过数值方法得到。严格的出发点是中值定理：

我们称为，显然如何精确的确定这一项成为中心任务。

欧拉法

普通欧拉法建议的是：

从而迭代公式为：

但这一取法是粗糙的，更为合理的是取和的斜率平均值：

但我们目前不知道，我们利用普通欧拉法得到的作为预报值，从而最终迭代公式为

或者我们采用中点估计：

从普通欧拉法，改进欧拉法到中点欧拉法提示我们通过渐进预报，即从预报到完成高精度逼近。

龙格-库塔法

龙哥-库塔法的表达式为：

其中：

这里表示第个点偏移的量。我们定义龙格库塔的阶数为，即到中取的点的数量。龙格库塔的关键是通过调整利用泰勒展开保证误差精度小于

二阶龙格库塔法

代入公式：

其中：

考虑展开

代入：

而泰勒展开为：

对比得到：

这个方程组有无数组解，显然欧拉法对应着，改进欧拉法对应着，还有中点公式：。

三阶和四阶的龙格库塔

三阶龙格库塔设为：

这里，同样利用泰勒展开得到：

现在我们考虑解这个方程组，首先考虑第三，第四个表达式：

这要求，解：

从而：

显然一个很对称的选择是代入得到，从而，于是：

而对于四阶，结果为：

变步长的龙格-库塔法

方才我们介绍的方法的局部截断误差为，对于阶龙格库塔为。值得注意的是，我们称此处为局部截断误差，局部指的是在假设为真实的情况下，的预测值相对于真实值的误差，即：

这里是步长，对于不同的，误差是不同的，例如缩小一半：

这里系数意味着多走了一倍的数量步长。因此步长减半，截断误差缩小了。得到：

这样我们通过计算：

可以直接估计误差的量级，而无需知道真实值。

当很大时，例如大于一个临界值，我们需要反复折半步长并且计算以至刚好。当很小时，例如小于一个临界值，我们需要反复增倍步长是的刚好大于以加速计算，并且取旧值（即刚好不超过)。这种策略称为变步长的龙格-库塔法，本质是搜索一个使得刚好小于的步长，本例中通过折半增倍完成。

2. 常微分方程组的初值问题

基本理论

我们研究常微分方程组的意义是解高阶常微分方程，例如：

设可以得到：

这就是一个微分方程组。

简而言之，任何阶微分方程都可以化成：

现在我们需要将其应用到龙格库塔法中，考虑矩阵形式：

得到：

现在已知，利用四阶龙格库塔法推出，利用：

这里：

Scipy

我们使用的用于求解常微分方程初值问题的函数为scipy.integrate.solve_ivp，这里ivp的全称是Initial Value Problem.

显然必须要传入的三个参数为，的范围，，值得注意的是的格式为：

为一个数组[y10, y20, ...]，我们得到：

可以看到output对象中中有message, success, status, t, y属性，这里有15个元素，粒度不够，可以通过t_eval参数指定序列。

还有一个技巧是，当func本身带有其他参数时，可以在solve_ivp中通过args=(a, b, c, ..)传入。

刚性问题

刚性问题是自变量缓慢变化时，因变量快速变化的问题，这样的问题需要采取高阶求解器或者极短的步长以获得精确解。

一个经典的问题是火柴燃烧模型：

此处为火球的半径，很容易利用solve_ivp得到：

考察最后50个点的标准差，在量级，这样的高误差是采取4阶龙格-库塔法求解器的结果。通过调整的大小可以改变量级：

或者采取更高阶的求解器，在solve_ivp函数中，利用参数method改变求解器，我们刚才采用的是RK45，它是一个改进版的龙格库塔4，5阶混合求解器。一个更高精度的求解器是Radau，它是隐式5阶龙格库塔变步长求解器。

可以看到，精度迅速提高。

可以获得其解析解为：

这里为Lambert W函数，在scipy.special模块中有该函数。

事件

在常微分方程演化时，可能遇到某些情形使得解进入非物理区域，因此需要设置事件中止演化。

在solve.ivp中，参数event用于控制这一行为，考虑event_func(t, y)，设置属性terminal=True表明当函数等于0时中止演化，设置属性direction=+1表示从负到正时停止演化，direction=-1表示从正到负停止演化。

一个经典的例子是自由落体：

3. 经典物理问题

刚体绕瞬心的转动方程

角动量定理

考虑对于点角动量定义为

求导：

对于第一项，考虑得到：

这里为外力，为总动量。得到：

注意到这个角动量定义中的采取的是绝对系下的动量。因此考虑到（位移定理），代入：

我们总是希望第一项为0，即要么，意味着是质心，或者，即是瞬心。

杆的滑动问题

未离开墙壁

这是经典问题，直接的思路是考虑利用牛顿第二定律对列出两个方程，利用几何关系对列出质心位置两个方程，再利用角动量定理对在点列出一个方程，

但是更佳的思路是寻找瞬心。考虑质心的位置：

因此给出瞬心满足的：

得到，对这个点列角动量定理，考虑到过，因此只有重力力矩：

这里得到：

最终得到：

恰巧离开墙壁

离开墙壁时，因此，得到：

考虑代入得到：

我们希望求出临界的，这在逻辑上可以通过直接演化状态1的方程组，当得到，其解析表达式需要通过能量守恒：

其中，从而得到：

完全离开墙壁

完全离开墙壁以后，几何关系不再成立，恒成立。考虑能量守恒：

有趣的是为什么仍然成立呢，实际上我们只是做出这个假设，如果我们计算出来有小于0的情况，则说明图像有问题。得到：

于是有：

这就是第二阶段的运动方程。

最后得到

圆锥陀螺运动

欧拉运动学方程

考虑转动过程

这个过程中很容易给出：

这里的上角标代表转动序号，我们需要全部写到本体系，即上角标为3。根据几何关系：

全部代入：

从而：

对称陀螺

考虑陀螺动能为：

考虑，代入得到：

显然，这对应着新的本体系为：

根据之前的推导，本体系为：，即只有旋动和章动。这就是说对称陀螺的有效旋转是中心对称轴绕绝对系纵轴的旋转和中心对称轴轴的章动，陀螺绕中心对称轴的自旋不带来新的物理结果。

势能为，因此得到：

圆锥陀螺

对于圆锥陀螺，可以直接计算出转动惯量，考虑高为，底面半径为，得到：

利用拉格朗日方程，得到最终的运动方程：

其中。

4. 边值问题和本征值问题

边值问题

考虑这样的边值问题，电荷分布满足：

电势满足泊松方程：

使用的是高斯单位制。由于球对称，，考虑得到：

现在考虑边界条件，在时，考虑电荷量很小，，在，考虑总电荷量贡献电势，。

直接积分法

直接积分法的思路是，猜测一侧边界，使得演化到另一边界时复合边界条件。在这个问题中，我们需要知道

这样演化的结果，如果边界值猜测不够准确，末态偏差是很大的。

该问题的解析解为，从而，于是

如果猜测的结果是，结果如下：

这就是直接积分法有限之处

两个初值问题

对于以下问题：

可以证明能够转化为：

因此原问题转化为：

本征值问题

本征值问题指的是这样一些微分方程，其中的参数必须取到特定的值才能得到满足边界条件的解。

弦振动方程

弦振动方程如下：

显然可以任意设置，我们需要找到正确的。我们试探性的打点可以得到

因此更精确的值可以通过二分法得到。

一维薛定谔方程的定态解

薛定谔方程为：

考虑势阱，在边界满足，考虑阱的形状为，则，，于是薛定谔方程可以写为：

边界条件满足，考虑使得，显然存在两个点，取其中一个如下图：

现在考虑从左侧出发，积分到得到，从右侧出发，积分到得到。通过选择归一化可以使得，我们要保证一阶导数连续：

考虑用后向：

或许就要求，。偏差为：

调整，让它变为0，考虑弦割法：

可以进行单点迭代：

除打靶法之外，对于薛定谔方程的谐振子势能，我们还可以考虑解析结果告知的第个本征态的波函数恰好有个节点（节点定理），考虑从左向右积分，如果节点数大于，则降低能量，如果节点数小于，则增加能量，从而得到本征值，再利用本征值进行中心对齐。

值得注意的是，如果直接从左向右积分，由于方程是初值敏感的，常常会发散，因此通常再次利用中心对齐微调。

用solve_bvp求边界值问题和本征值问题

solve_bvp的基本语法是solve_bvp(func,bc,x,y,para,S)，这里func就是问题服从的微分方程，bc是问题服从的边界条件，值得注意的是，bvp需要同时传入坐标架x和一个y的初值，例如0向量。因此对于func中的，传入它们的是向量，得到的是个向量，需要将它们用np.vstack()拼接起来。而bc要求传入的是在两个边界满足的条件,即为。

例如对于，代码为：

对于本征值问题，可以通过迭代本征值来观察收敛情况，收敛的好的就是本征值。

我们还关心当方程中有奇点的情形，这就是参数S的用武之地，考虑艾登方程：

显然是这个方程的奇点，有条件保证方程有解。现在设，则：

我们称，严格地讲，对于，且，可用solve_bvp解

最后，对于突变问题，通过参数fun_jac, bc_jac，传入方程的雅可比解析矩阵和边界条件的雅可比矩阵可以提高精度。例如对于问题：

方程雅可比为：

边界条件为：

边界条件雅可比为：