python - 关于常微分方程 (ODE) 和优化，在 Python 中

Question

我想解决这样的问题:

dy/dt = 0.01*y*(1-y), find t when y = 0.8 (0<t<3000)

我在 Python 中尝试过 ode 函数，但它只能在给定 t 时计算 y。

那么有没有什么简单的方法可以在 Python 中解决这个问题呢？

PS:这个函数只是一个简单的例子。我真正的问题是如此复杂，无法通过分析来解决。所以我想知道如何数值求解。而且我认为这个问题更像是一个优化问题:

Objective function y(t) = 0.8, Subject to dy/dt = 0.01*y*(1-y), and 0<t<3000 

PPS:我真正的问题是:

objective function: F(t) = 0.85, 
subject to: F(t) = sqrt(x(t)^2+y(t)^2+z(t)^2), 
            x''(t) = (1/F(t)-1)*250*x(t), 
            y''(t) = (1/F(t)-1)*250*y(t), 
            z''(t) = (1/F(t)-1)*250*z(t)-10, 
            x(0) = 0, y(0) = 0, z(0) = 0.7, 
            x'(0) = 0.1, y'(0) = 1.5, z'(0) = 0, 
            0<t<5

Answer 1

这个微分方程可以很容易地解析求解:

dy/dt = 0.01 * y * (1-y)

重新排列以收集相反侧的 y 和 t 项

100 dt = 1/(y * (1-y)) dy

lhs 简单地集成到 100 * t，rhs 稍微复杂一些。我们总是可以将两个商的乘积写为两个商之和 * 一些常数:

1/(y * (1-y)) = A/y + B/(1-y)

A 和 B 的值可以通过将 rhs 放在相同的分母上并比较两边的常数项和一阶 y 项来计算出来。在这种情况下很简单，A=B=1。因此我们必须整合

1/y + 1/(1-y) dy

第一项积分到 ln(y)，第二项可以通过变量 u = 1-y 到 -ln(1-y) 的变化积分。我们的综合方程如下:

100 * t + C = ln(y) - ln(1-y)

不要忘记积分常数(这里写在lhs上很方便)。我们可以结合两个对数项:

100 * t + C = ln( y/(1-y) )

为了求解 t 的精确值 y，我们首先需要计算出 C 的值。我们使用初始条件来完成此操作。很明显，如果 y 从 1 开始，则 dy/dt = 0 并且 y 的值永远不会改变。因此在开头插入 y 和 t 的值

100 * 0 + C = ln( y(0)/(1 - y(0) )

这将为 C 提供一个值(假设 y 不是 0 或 1)，然后使用 y=0.8 来获得 t 的值。请注意，由于对数和因子 100 乘以 t y 将在相对较短的 t 值范围内达到 0.8，除非 y 的初始值非常小。当然，重新排列上面的方程以用 t 表示 y 也很简单，然后您也可以绘制函数。

编辑:数值积分

对于无法解析求解的更复杂 ODE，您将不得不尝试数值计算。最初我们只知道函数在零时间 y(0) 的值(我们至少必须知道它才能唯一地定义函数的轨迹)，以及如何评估梯度。数值积分的想法是我们可以利用我们对梯度的知识(它告诉我们函数如何变化)来计算出函数在我们起点附近的值。最简单的方法是欧拉积分:

y(dt) = y(0) + dy/dt * dt

欧拉积分假设梯度在 t=0 和 t=dt 之间是常数。一旦 y(dt) 已知，也可以在那里计算梯度，进而用于计算 y(2 * dt) 等，逐渐建立函数的完整轨迹。如果您正在寻找特定的目标值，只需等到轨迹超过该值，然后在最后两个位置之间进行插值以获得精确的 t。

欧拉积分(以及所有其他数值积分方法)的问题在于其结果仅在其假设有效时才是准确的。由于时间点对之间的梯度不是恒定的，因此每个积分步骤都会产生一定量的误差，随着时间的推移，这些误差会逐渐累积，直到答案完全不准确。为了提高积分的质量，有必要对梯度使用更复杂的近似值。例如，查看 Runge-Kutta 方法，它们是一系列积分器，以增加计算时间为代价去除误差项的渐进阶数。如果你的函数是可微的，知道二阶甚至三阶导数也可以用来减少积分误差。

当然幸运的是，其他人已经在这里完成了艰苦的工作，您不必太担心解决数值稳定性等问题或对所有细节有深入的了解(尽管大致了解正在发生的事情会很有帮助) )。退房 http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode有关您应该能够直接使用的积分器类的示例。例如

from scipy.integrate import ode
def deriv(t, y):
    return 0.01 * y * (1 - y)
my_integrator = ode(deriv)
my_integrator.set_initial_value(0.5)
t = 0.1  # start with a small value of time
while t < 3000:
    y = my_integrator.integrate(t)
    if y > 0.8:
        print "y(%f) = %f" % (t, y)
        break
    t += 0.1

此代码将在 y 超过 0.8 时打印出第一个 t 值(如果从未达到 0.8，则不打印)。如果您想要更准确的 t 值，请保留前一个 t 的 y 并在它们之间进行插值。