scipy中的复杂ODE系统

Question

我在求解光学布洛赫方程时遇到了问题，这是一个具有复数值的一阶 ODE 系统。我发现 scipy 可以解决这样的系统，但是他们的网页提供的信息太少，我几乎看不懂。

我有 8 个耦合的一阶 ODE，我应该生成如下函数：

def derv(y):
    compute the time dervative of elements in y
    return answers as an array

然后做complex_ode(derv）

我的问题是：

1. 我的 y 不是列表而是矩阵，我怎样才能给出正确的输出 适合 complex_ode() 吗？
2. complex_ode() 需要一个雅可比，我不知道如何开始构建一个 应该是什么类型？
3. 我应该把初始条件放在哪里，就像在正常颂歌和 时间线空间？

这是 scipy 的 complex_ode 链接： http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.complex_ode.html

谁能给我更多的信息，以便我能学到更多。

Answer 1

我认为我们至少可以为您指明正确的方向。光学的 布洛赫方程是一个科学界很好理解的问题 社区，虽然不是我的:-)，所以互联网上已经有解决方案 针对这个特殊问题。

http://massey.dur.ac.uk/jdp/code.html

但是，为了满足您的需求，您谈到了使用 complex_ode，我想这是 很好，但我认为简单的 scipy.integrate.ode 也可以正常工作 根据他们的文档：

 from scipy import eye
 from scipy.integrate import ode

 y0, t0 = [1.0j, 2.0], 0

 def f(t, y, arg1):
     return [1j*arg1*y[0] + y[1], -arg1*y[1]**2]
 def jac(t, y, arg1):
     return [[1j*arg1, 1], [0, -arg1*2*y[1]]]
 r = ode(f, jac).set_integrator('zvode', method='bdf', with_jacobian=True)
 r.set_initial_value(y0, t0).set_f_params(2.0).set_jac_params(2.0)
 t1 = 10
 dt = 1
 while r.successful() and r.t < t1:
     r.integrate(r.t+dt)
     print r.t, r.y

你还有一个额外的好处，那就是更老的更成熟和更好的 记录的功能。我很惊讶你有 8 个而不是 9 个耦合 ODE，但我 确保您比我更了解这一点。是的，您是正确的，您的功能 应该是 ydot = f(t,y) 的形式，你称之为 def derv() 但你是 需要确保你的函数至少有两个参数 喜欢derv(t,y)。如果您的y 在矩阵中，没问题！只需“重塑”它 derv(t,y) 函数如下：

Y = numpy.reshape(y,(num_rows,num_cols));

只要num_rows*num_cols = 8，你的ODE数量应该没问题。然后 在计算中使用矩阵。完成后，请务必返回 一个向量，而不是像这样的矩阵：

out = numpy.reshape(Y,(8,1));

Jacobian 不是必需的，但它可能会允许计算继续进行 快得多。如果您不知道如何计算，您可能需要咨询 维基百科或微积分教科书。这很简单，但可能很耗时。

就初始条件而言，您可能应该已经知道那些应该是什么 是，无论它是复杂的还是真正有价值的。只要您选择的值 在合理的范围内，这应该无关紧要。