如何使odeint成功? [英] How to make odeint successful?

问题描述

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

<预><代码>>>>在此调用上完成了过多的工作(可能是错误的 Dfun 类型).使用 full_output = 1 运行以获取定量信息.>>>

因此，我必须更改我的时间步长和最终时间，以使其可集成.为此，我需要尝试不同的组合，这非常痛苦.谁能告诉我如何让 odeint 自动改变时间步长和最终时间以成功集成这个 ode 系统?

这是调用odeint的代码的一部分:

def main(t, init_pop_a, init_pop_b, *args, **kwargs):"""求解给定参数集的 obe"""# 构造初始条件# 最初，rho_ee = 0rho_init = zeros((16,16))*1j ########rho_init[1,1] = init_pop_arho_init[2,2] = init_pop_brho_init[0,0] = 1 - (init_pop_a + init_pop_b)########rho_init_ravel, params = to_1d(rho_init)# 执行整合结果 = odeint(wrapped_bloch3，rho_init_ravel，t，args=args)# BUG:需要通过 kwargs# 重新包装结果返回 from_1d(result, params, prepend=(len(t),))东西 = [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 = 东西args = ( Delta_a, Delta_b, Omega_a, Omega_b )t = arange(0, tfinal + tstep, tstep)数据 = main(t, init_pop_a, init_pop_b, *args)plt.plot(t,abs(data[:,4,4]))

其中，wrapped_bloch3 是计算 dy/dt 的函数.

解决方案

EDIT:我注意到你已经在这里得到了答案:scipy 中的复杂 ODE 系统

odeint 不适用于复数值方程.我得到

from scipy.integrate import odeint将 numpy 导入为 np定义函数(t，y):返回 1 + 1jt = np.linspace(0, 1, 200)y = odeint(func, 0, t)# -> 这输出:## TypeError: 不能将复数转换为浮点数# odepack.error: 函数调用的结果不是正确的浮点数组.

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

from scipy.integrate import ode将 numpy 导入为 npdef myodeint(func, y0, t):y0 = np.array(y0, 复数)func2 = lambda t, y: func(y, t) # odeint 有其他方式:/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)返回 np.array(y)def func(y, t, alpha):返回 1j*alpha*y阿尔法 = 3.3t = np.linspace(0, 1, 200)y = myodeint(lambda y, t: func(y, t, alpha), [1, 0, 0], t)

I am a python beginner, currently using scipy's odeint to compute a coupled ODE system, however, when I run, python shell always tell me that

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

So, I have to change my time step and final time, in order to make it integratable. To do this, I need to try a different combinations, which is quite a pain. Could anyone tell me how can I ask odeint to automatically vary the time step and final time to successfully integrate this ode system?

and here is part of the code which has called 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]))

where wrapped_bloch3 is the function compute dy/dt.

解决方案

EDIT: I note you already got an answer here: complex ODE systems in scipy

odeint does not work with complex-valued equations. I get

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.

You can solve your equation by the other ode solver:

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)

这篇关于如何使odeint成功?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！