2

我想在 python 中解决大量的双线性 ODE 系统。导数是这样的：

def x_(x, t, growth, connections):
    return x * growth + np.dot(connections, x) * x

我对非常准确的结果不感兴趣，而是对定性行为感兴趣，即组件是否为零。

因为我要解决如此大量的高阶系统，所以我想使用尽可能大的步长。

由于步长较大，ODE 可能会出现在零以下的一个分量中。这应该是不可能的，因为（由于特定 ODE 的结构）每个分量都以零为界。因此 - 为了防止错误的结果 - 我想在每个组件低于时手动将其设置为零。

此外，在我想要解决的系统中，解决方案可能会崩溃。我还想通过设置上限来防止这种情况发生，即如果一个值超过了界限，它会被设置回界限的值。

我希望我可以让我的目标易于理解，并为您提供以下我想要做的伪代码：

for t in range(0, tEnd, dt):
    $ compute x(t) using x(t-dt) $
    x(t) = np.minimum(np.maximum(x(t), 0), upperBound)

我使用 Runge-Kutta 算法实现了这一点。一切正常。只是性能很差。因此，我更喜欢使用像 scipy.integrate.odeint 这样的预实现方法。

因此，我不知道如何设置这样的界限。我尝试过的一个选项是以这种方式操纵 ODE，一旦 x 高于界限，导数变为 0，一旦 x 低于 0，则导数变为 0。此外，为了防止在一个时间步内跳得太高，我还有界导数：

def x_(x, t, growth, connections, bound):
    return (x > 0) * np.minimum((x < bound) * \
           ( x * growth + np.dot(connections, x) * x ), bound) + (x < 0)

虽然这个解决方案（特别是对于零界限）非常难看，但如果它有效就足够了。不幸的是，它不起作用。使用 odeint

x = scipy.integrate.odeint(x_, x0, timesteps, param)

我经常遇到以下两个错误之一：

Repeated convergence failures (perhaps bad Jacobian or tolerances).

Excess work done on this call (perhaps wrong Dfun type).

它们可能是由于我操纵的 ODE 的不连续性。互联网上有很多关于这些错误消息的主题，但它们对我没有帮助。例如，增加允许的步数既不能阻止这个问题，也不是一个好的解决方案，因为我需要使用大步长。此外，通过雅可比也无济于事。

查看解决方案可以看到发生错误时会发生两种奇怪的行为：

  1. 该解决方案在一个时间步长内达到 +-1e250（这应该是不可能的，因为 dx/dt 是有界的）。
  2. 它首先到达边界，但又下降（这应该是不可能的，因为 x 处于边界，因此 x_ 为 0）。

我将不胜感激有关如何解决问题的所有提示-无论它是否有助于

  • 如何防止odeint中的错误
  • 如何正确或打开 ODE 操作
  • 如何编写一个非常快速的 ODE 求解器，我可以直接实现我的需求。

我提前谢谢你！

编辑

我被问到一个最小的例子：

import numpy as np
import random as rd
rd.seed()
import scipy.integrate

def simulate(simParam, dim = 20, connectivity = .8, conRange = 1, threshold = 1E-3, 
             constGrowing=None):
    """
    Creates the random system matrix and starts a simulation
    """

    x0 = np.zeros(dim, dtype='float') + 1
    connections = np.zeros(shape=(dim, dim), dtype='float')

    growth = np.zeros(dim, dtype='float') + 
                    (constGrowing if not constGrowing == None else 0)


    for i in range(dim):        
        for j in range(dim):
            if i != j: 
                connections[i][j] = rd.uniform(-conRange, conRange)

    tend, step = simParam
    return RK4NumPy(x_, (growth, connections), x0, 0, tend, step)


def x_(x, t, growth, connections, bound):
    """
    Derivative of the ODE
    """
    return (x > 0) * np.minimum((x < bound) * 
                (x * growth + np.dot(connections, x) * x), bound) + (x < 0)


def RK4NumPy(x_, param, x0, t0, tend, step, maxV = 1.0E2, silent=True):
    """
    solving method
    """
    param = param + (maxV,)
    timesteps = np.arange(t0 + step, tend, step)

    return scipy.integrate.odeint(x_, x0, timesteps, param)


simulate((300, 0.5))

要查看解决方案，必须绘制 x。使用给定的参数，我经常得到上述错误

Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.
4

0 回答 0