如何让odeint成功？答案

【问题标题】：How to make odeint successful?如何让odeint成功？
【发布时间】：2012-03-16 23:40:42
【问题描述】：

我是python初学者，目前使用scipy的odeint计算耦合ODE系统，但是，当我运行时，python shell总是告诉我

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

所以，我必须更改我的时间步长和最终时间，以使其可积分。为此，我需要尝试不同的组合，这很痛苦。谁能告诉我如何让odeint 自动更改时间步长和最终时间以成功集成此 ode 系统？

这里是调用 odeint 的部分代码：

def main(t, init_pop_a, init_pop_b, *args, **kwargs):
    """
    solve the obe for a given set of parameters
    """
    # construct initial condition
    # initially, rho_ee = 0
    rho_init = zeros((16,16))*1j ########
    rho_init[1,1] = init_pop_a
    rho_init[2,2] = init_pop_b
    rho_init[0,0] = 1 - (init_pop_a + init_pop_b)########
    rho_init_ravel, params = to_1d(rho_init)
    # perform the integration
    result = odeint(wrapped_bloch3, rho_init_ravel, t, args=args)
                        # BUG: need to pass kwargs
    # rewrap the result
    return from_1d(result, params, prepend=(len(t),))

things = [2*pi, 20*pi, 0,0, 0,0, 0.1,100]
Omega_a, Omega_b, Delta_a, Delta_b, \
init_pop_a, init_pop_b, tstep, tfinal = things
args = ( Delta_a, Delta_b, Omega_a, Omega_b )
t = arange(0, tfinal + tstep, tstep)
data = main(t, init_pop_a, init_pop_b, *args)

plt.plot(t,abs(data[:,4,4]))

其中 Wrapped_bloch3 是函数计算 dy/dt。

【问题讨论】：

能否提供更多代码，尤其是对 odeint 的调用？
您将不得不添加比您提供的更多的细节来获得帮助：您正在使用哪种类型的 ODE？他们僵硬吗？您是否向odeint 提供雅可比函数？你确定这是合理的吗？
感谢重播，我已经更新了我的问题:)
@user1233157：最明显的做法是为求解器定义一个雅可比项。这应该会大大加快收敛速度。
@talonmies：谢谢你的建议，但是，我的耦合 ode 系统有 16*16 个方程，这可能需要我很长时间才能写出他们的雅可比行列式。我查看了 odeint 的文档，他们在成功时有一些类似的东西，看起来可能能够解决我的问题，但我就是无法理解。

标签： python scipy

【解决方案1】：

编辑：我注意到你已经在这里得到了答案：complex ODE systems in scipy

odeint 不适用于复值方程。我明白了

from scipy.integrate import odeint
import numpy as np
def func(t, y):
    return 1 + 1j
t = np.linspace(0, 1, 200)
y = odeint(func, 0, t)
# -> This outputs:
#
# TypeError: can't convert complex to float
# odepack.error: Result from function call is not a proper array of floats.

您可以通过其他 ode 求解器求解方程：

from scipy.integrate import ode
import numpy as np

def myodeint(func, y0, t):
    y0 = np.array(y0, complex)
    func2 = lambda t, y: func(y, t)   # odeint has these the other way :/
    sol = ode(func2).set_integrator('zvode').set_initial_value(y0, t=t[0])
    y = [sol.integrate(tp) for tp in t[1:]]
    y.insert(0, y0)
    return np.array(y)

def func(y, t, alpha):
    return 1j*alpha*y

alpha = 3.3
t = np.linspace(0, 1, 200)
y = myodeint(lambda y, t: func(y, t, alpha), [1, 0, 0], t)

【讨论】：