2

我想解决这样的问题：

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
4

4 回答 4

2

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

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 乘以 ty 将在相对较短的 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 并在它们之间进行插值。

于 2013-08-22T15:15:50.103 回答
1

作为对 Krastanov 的回答的补充：

除了 PyDSTool 之外，还有其他包，例如Pysundials和Assimulo，它们提供了从 Sundial 到求解器IDA的绑定。该求解器具有寻根功能。

于 2015-05-24T22:32:33.830 回答
0

用于scipy.integrate.odeint处理您的集成，并在之后分析结果。

import numpy as np
from scipy.integrate import odeint

ts = np.arange(0,3000,1)   # time series - start, stop, step

def rhs(y,t):
    return 0.01*y*(1-y)

y0 = np.array([1])    # initial value

ys = odeint(rhs,y0,ts)

然后分析 numpy 数组ys以找到您的答案（数组ts匹配的维度ys）。（这可能不是第一次工作，因为我是从记忆中构建的）。

这可能涉及使用数组的scipy interpolate函数ys，以便您在 time 获得结果t。

编辑：我看到您希望解决 3D 中的弹簧问题。使用上述方法应该没问题；Odeint网站上scipy有可以解决的系统示例，例如耦合弹簧，这些可以扩展。

于 2013-08-22T14:10:33.277 回答
0

您需要的是具有根查找功能的 ODE 集成器。它们存在并且此类集成器的低级代码由 scipy 提供，但它们尚未包装在 python 绑定中。

有关更多信息，请参阅此邮件列表帖子，其中提供了一些替代方案：http: //mail.scipy.org/pipermail/scipy-user/2010-March/024890.html

您可以使用以下使用回溯的示例实现（因此它不是最佳的，因为它是对自身没有根查找的集成器的附加补充）：https ://github.com/scipy/scipy/拉/4904/文件

于 2015-05-24T15:13:36.847 回答