为什么 scipy 的 odeint 函数为一个解决方案应该是单调的问题提供一个非单调的解决方案？

计算科学 Python 数值分析颂

2021-12-12 16:10:14

下面 ode 的解看起来是单调递增的：

在此处输入图像描述

然而，仔细观察，我们发现它不是：

在此处输入图像描述

如何确保数值解单调递增？

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def f1(theta):
    return -1*np.cos(theta)/np.tan(theta)
def f2(theta):
    return -1*np.cos(theta)

def f3(y, t) :
    return 1/np.sqrt(f1(y)**2+f2(y)**2)

s = np.linspace(0,30,100)
theta = odeint(f3, -np.pi/2, s)

plt.plot(s, theta)
plt.show()

编辑 1

正如 Geoff Oxberry 在他的回答中提到的那样，ODE 的 RHS 在初始条件下是无限的。所以我想找到另一个点来开始整合，它应该与上述原始案例的轨迹相同。我现在解决 $\frac{ds}{d\theta} = f_3(\theta)$ 代替 $\frac{d\theta}{ds} = \frac{1}{f_3(\theta)}$ . 前者在初始条件下不具有奇异性 $s=0$ 在 $\theta = -\pi/2$ . 然后我整合到任意 $\theta = -\pi/4$ , 并得到 $s$ 那里。然后我使用这一点作为我在原始问题中的初始条件。

def f4(theta):
    return np.sqrt(f1(theta)**2+f2(theta)**2)

# Get another initial condition

theta0 = -np.pi/4
s0 =integrate.quad(f4, -np.pi/2, theta0)

但是我仍然遇到这种情况下的振荡问题，所以我不确定如何使用这种方法。

2个回答

如果您希望数值方法满足某些属性，则必须证明（或以这种方式构造方法）这些属性成立。

例子：

守恒定律的保守有限体积公式
Runge-Kutta 方法保留线性不变量
辛积分器保持辛性，这对哈密顿系统很重要
辛 Runge-Kutta 方法保留二次不变量
总变差递减离散化 (TVD) 是保持单调性的（通常，这些是双曲问题的空间离散化）
时间上的强稳定性保持离散化是那些保持 TVD 属性的离散化，如果空间离散化在使用前向 Euler 进行时间积分时是 TVD
您可以强制某些积分方案（例如：SUNDIALS 中的 BDF 实现）尊重 ODE 解决方案属性的界限（例如：您正在积分 ODE，跟踪物体的质量；该物体的质量必须是非负的）

正如 DavidKetcheson 正确指出的那样，收紧公差不一定会使数值解单调。没有办法保证这一点。总结以下讨论，在实际意义上，如果您正在单调地接近稳态，收紧公差将限制虚假振荡（具体来说，它会限制它们的幅度，但不一定会限制它们的总变化），但不会消除它们对于某些。

您的特定 ODE 是：

\begin{aligned} \dot{θ} = \frac{1}{\sqrt{\frac{\cos^{4} θ}{\sin^{2} θ} + \cos^{2}} θ} \end{aligned}

$\begin{align} \dot{\theta} = \frac{1}{\sqrt{\frac{\cos^{4}{\theta}}{\sin^{2}{\theta}} + \cos^{2}}{\theta}} \end{align}$

经过一些简化和咨询三角恒等式，它变成

\begin{aligned} \dot{θ} = \pm \tan θ, \end{aligned}

$\begin{align} \dot{\theta} = \pm\tan{\theta}, \end{align}$

自从你开始 $\theta(0) = -\pi/2$ ，你的右手边无论如何都是单数的，你的“解决方案”是可疑的。（符号通过歧义的原因 $\pm$ 在余切方面涉及正弦的三角恒等式，符号的选择取决于哪个象限 $\theta$ 是。）您看到解决方案的唯一原因是浮点运算的怪癖；当我将您的右侧评估为零时，我得到 16331239353195370.0，这可能是与系统和硬件相关的结果。

为了 $\theta \in (-\pi/2, 0]$ ，你应该有 $\dot{\theta} = -\tan{\theta}$ ，并且由于 $-\tan{\theta}$ 是非负的 $(-\pi/2, 0]$ 和 $\tan{0} = 0$ , 任何初始条件 $(-\pi/2, 0]$ 应该产生一个单调递增的解决方案。同样，您使用的数值方法不能保证此属性。此外，如果 $\theta \in (0, \pi/2)$ ，你的方程变成 $\dot{\theta} = \tan{\theta}$ 保持右手边的正性（等式的原始形式确保正性），和 $\tan{\pi/2}$ 当然，又是单数。

您的容差在这里的主要功能是确保您获得正确的稳态，并且由于浮点运算是神秘的，您以某种方式得到看起来像“正确”的答案（对于 $\theta(0) = -\pi/2 + \varepsilon$ 对于一些小 $\varepsilon > 0$ ) 尽管初始条件处于奇点，但您会看到与真实稳态的偏差 $\theta = 0$ 是因为就计算解决方案中的估计误差而言，容差对应于“大致为零”。由于浮点运算中的怪癖，您的右侧不一定会计算为零，但由于它相对于您的公差足够小，您应该不会看到幅度大于您的公差的振荡（大约，因为这种推理是启发式的，并不严谨，但它应该成立）。

除了稳定性或必要精度的原因外，我对极小的公差持怀疑态度，除非由于物理或数字原因需要它们。对准确性的过度要求将不必要地限制时间步长并增加模拟执行时间，而不会产生任何明显的结果改进。作为一个反例，燃烧问题因需要对某些物种的小公差来获得近似正确的数值结果而臭名昭著。如果某些物种从未存在，则不会发生某些反应，即使这些物种的浓度是以不可能测量的数量计算的，因此在这种情况下，有必要保证小的绝对和相对公差。对于您的特定问题，我怀疑有足够大的公差， $\theta = \pi/2$ ，所以默认容差可能很好，但我认为没有充分的理由降低这里的容差。

我建议降低 odeint 的容忍度。将“rtol”或“atol”参数设置为小于默认值 1.49012e-8。在此处查看 odeint 的 scipy 文档。

其它你可能感兴趣的问题

上一篇偏微分方程数值解的数学检验方法？下一篇我应该考虑为专为科学计算构建的台式机购买哪种处理器？