scipy中的复杂ODE系统

Question

我在解光学布洛赫方程时遇到困难，这是一个具有复值的一阶常微分方程系统。我发现scipy可以解决这样的系统，但是他们的网页提供的信息太少，我很难理解。

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

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

然后complex_ode(derv）

我的问题是：

1. 我的y不是一个列表，而是一个矩阵，我如何才能给出相应的输出 适合复杂的代码吗？
2. complex_ode()需要jacobian，我不知道如何开始构造jacobian 它应该是什么类型的？
3. 我应该把初始条件放在哪里，比如在正常的ode和 时间紧吗？

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

有谁能给我提供更多的信息以便我能学到更多。

Answer 1

我想我们至少可以给你指出正确的方向。光学的 布洛赫方程是科学界很清楚的一个问题 社区，虽然不是由我来做的：-），所以互联网上已经有了解决方案 对于这个特殊的问题。

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

但是，为了满足您的需要，您提到了使用复杂代码，我想是这样的 很好，但是我认为简单的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，但它可能允许计算继续进行 更快。如果你不知道如何计算，你可能想咨询一下 维基百科或微积分课本。这很简单，但可能很费时。

至于初始条件，你应该已经知道 无论是复杂的还是真实的价值。只要选择的值是 在合理的范围内，这应该不重要。