标签: odeint

为了使用SciPy解决简单的ODE,我曾经使用odeint函数,其形式为:

scipy.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)[source]

要集成的简单函数可以包含表单的其他参数:

def dy_dt(t, y, arg1, arg2):
    # processing code here

在SciPy 1.0中,似乎ode和odeint函数已被更新的solve_ivp方法所取代.

scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, **options)

但是,这似乎没有提供args参数,也没有提供文档中关于实现args传递的任何指示.

因此,我想知道是否可以使用新API进行arg传递,或者这是一个尚未添加的功能？(如果故意删除这些功能,对我来说似乎是一种疏忽？)

参考:https: //docs.scipy.org/doc/scipy/reference/integrate.html

我想在多个线程(每个CPU核心一个)中使用scipy.integrate.ode(或scipy.integrate.odeint)实例,以便一次解决多个IVP.然而文档说:" 这个集成器不是可重入的.你不能同时使用"vode"集成器的两个ode实例. "

(如果文档没有说明,odeint会在多次实例化时导致内部错误.)

知道可以做些什么吗？

任何人都可以提供一个integrate.odeint在SciPy中为函数提供jacobian的例子吗？我尝试从SciPy教程odeint示例运行此代码,但似乎永远不会调用Dfunc(渐变).

from numpy import * # added
from scipy.integrate import odeint
from scipy.special import gamma, airy
y1_0 = 1.0/3**(2.0/3.0)/gamma(2.0/3.0)
y0_0 = -1.0/3**(1.0/3.0)/gamma(1.0/3.0)
y0 = [y0_0, y1_0]


def func(y, t):
    return [t*y[1],y[0]]


def gradient(y,t):
    print 'jacobian'  # added
    return [[0,t],[1,0]]


x = arange(0,4.0, 0.01)
t = x
ychk = airy(x)[0]
y = odeint(func, y0, t)
y2 = odeint(func, y0, t, Dfun=gradient)
print y2 # added

使用numba.jit加快了右手端计算odeint从scipy.integrate正常工作：

from scipy.integrate import ode, odeint
from numba import jit

@jit
def rhs(t, X):
    return 1

X = odeint(rhs, 0, np.linspace(0, 1, 11))

但是这样使用integrate.ode：

solver = ode(rhs)
solver.set_initial_value(0, 0)
while solver.successful() and solver.t < 1:
    solver.integrate(solver.t + 0.1)

装饰器产生以下错误@jit：

capi_return is NULL
Call-back cb_f_in_dvode__user__routines failed.
Traceback (most recent call last):
  File "sandbox/numba_cubic.py", line 15, in <module>
    solver.integrate(solver.t + 0.1)
  File "/home/pgermann/Software/anaconda3/lib/python3.4/site-packages/scipy/integrate/_ode.py", line 393, in integrate
    self.f_params, self.jac_params)
  File "/home/pgermann/Software/anaconda3/lib/python3.4/site-packages/scipy/integrate/_ode.py", line …

编辑：这个问题是由于一个错误引起的，已在 Scipy 0.15 中修复

当我开发和测试代码时，我可能会犯一个简单的错误，例如 NameError。当我使用 scipy.integrate.odeint 时，odeint 将打印错误消息，但会继续集成我请求的多个时间步，因此我收到许多相同的错误消息。我认为它具有这种行为，以便在发生算术错误（例如除以零）时它可以继续进行，但这对于编程错误来说是无益的行为。

有没有办法让 scipy 在第一条错误消息后停止？如果我能让它因错误而不是算术异常而停止，那就最好了。

我有两个 numpy 数组：9x9 和 9x1。我想求解离散时间点的微分方程，但在让 ODEInt 工作时遇到问题。我确实不确定我是否在做正确的事。

对于 Mathematica，方程是：

Solution = {A[t]} /. NDSolve[{A'[t] == Ab.A[t] && A[0] == A0}, {A[t]}, {t, 0, .5}, MaxSteps -> \[Infinity]];

time = 0.25;
increment = 0.05;
MA = Table[Solution, {t, 0, time, increment}];

其中 Ab 是 9x9 矩阵，A0 是 9x1 矩阵（初始）。在这里，我解决了时间问题，生活很美好。

在Python实现中，我有以下代码，它给了我错误的答案：

from scipy.integrate import odeint
from numpy import array, dot, pi

def deriv(A, t, Ab):
    return dot(Ab, A)

def MatrixBM3(k12,k21,k13,k31,k23,k32,delta1,delta2,delta3,
              w1, R1, R2):

  K = array([[-k21 -k23, k12, k32, 0., 0., 0., …

我正在尝试使用openmp来并行化我的代码。当我使用恒定的步长时，一切工作都很好，但是当我使用自适应步进器运行相同的代码时，会出现我不理解的错误。

这是代码的基本部分：

using namespace std;
using namespace boost::numeric::odeint;
const int jmax = 10;
typedef  double value_type;
typedef boost::array<value_type ,2*(jmax+1) > state_type;


//The step function
void rhs( const state_type A , state_type &dAdt , const  value_type t )

{

 value_type RHStemp0;
 value_type RHStemp1;
//We will write the RHS of the equations a  big sum
#pragma omp parallel for schedule(runtime)
for(int j = 0; j < jmax+1 ; j++ ) //Real part
{
    RHStemp0 = value_type(0.0);
    RHStemp1 = value_type(0.0);
    for (int k …

我有一个普遍的问题,我将在一个更具体的情况下说.

如果想要找到双摆的动力学,可以在数学上推导出运动方程,重写ODE以具有对数值计算有用的特殊形式,并使用c ++中的odeint来解决ODE(参见此示例)在堆栈溢出/sf/answers/2140791901/).

现在想象一下我们想对n个耦合摆(n在运行时已知)做同样的事情.这需要我们写一个所谓的拉格朗日(动能 - 势能),这个函数的不同衍生物将是我们需要解决的ODE.此外,这些ODE必须以适合odeint的形式重写.一般来说,这很难用手去做.

在像Mathematica和Maple这样的程序中,这实际上非常简单.人们可以明确地从拉格朗日求出所需的ODE,并且ODE求解器不需要我们将方程放在任何特殊形式中(参见mathematica中的示例https://mathematica.stackexchange.com/a/84279).

是否有可能在c ++中做这样的事情而不会经历太多麻烦？

潜在方法:

一种可能的方法是使用c ++包ginac.这可以帮助我们分析地推导出ODE.但我不知道如何将来自ginac的表达式重写为适合odeint中数值计算的形式.有任何想法吗？

下面的脚本集成了封闭路径周围的磁场线,并在使用Python中的Runge-Kutta RK4在某个容差范围内恢复到原始值时停止.我想用SciPy.integrate.odeint,但我看不出当路径近似关闭时我怎么能告诉它停止.

当然odeint可能比在Python中集成要快得多,我可以让它盲目地绕过并在结果中寻找结束,但在将来我会做更大的问题.

有没有办法可以实现" OK,足够接近 - 你现在可以停止! "方法进入odeint？或者我应该整合一段时间,检查,整合更多,检查......

这个讨论似乎很相关,似乎暗示"你不能从SciPy内部"可能是答案.

注意:我通常使用RK45(Runge-Kutta-Fehlberg),它在给定的行程尺寸下更准确,以加快速度,但我在这里保持简单.它还使可变步长成为可能.

更新:但有时我需要固定的步长.我发现它Scipy.integrate.ode确实提供了一种测试/停止方法,ode.solout(t, y)但似乎没有能力在固定点进行评估t.odeint允许在固定点进行评估t,但似乎没有测试/停止方法.

def rk4Bds_stops(x, h, n, F, fclose=0.1):

    h_over_two, h_over_six = h/2.0, h/6.0

    watching = False
    distance_max = 0.0
    distance_old = -1.0

    i = 0

    while i < n and not (watching and greater):

        k1 = F( x[i]                )
        k2 = F( x[i] + k1*h_over_two)
        k3 = F( x[i] + k2*h_over_two)
        k4 = …

SciPy 可以通过 scipy.integrate.odeint 或其他包求解 ode 方程，但它在函数完全求解后给出结果。但是，如果 ode 函数非常复杂，程序将花费大量时间（一两天）才能给出整个结果。那么我如何监控求解方程的步骤（当方程尚未完全求解时打印出结果）？