复杂的ODE系统 [英] complex ODE systems in scipy

问题描述

我在求解光学布洛方程时遇到了麻烦，该方程是具有复杂值的一阶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/generation/scipy.integrate.complex_ode.html

任何人都可以向我提供更多信息，以便我可以学到更多.

解决方案

我认为我们至少可以为您指明正确的方向.光学 bloch方程是一个在科学界很容易理解的问题 社区，虽然不是我自己:-)，所以互联网上已经有解决方案 解决这个特定的问题.

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));

不需要雅可比行列式，但是很可能允许进行计算 快得多.如果您不知道如何计算，可能需要咨询 维基百科或微积分教科书.这很简单，但可能很耗时.

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

I am having trouble sovling the optical bloch equation, which is a first order ODE system with complex values. I have found scipy may solve such system, but their webpage offers too little information and I can hardly understand it.

I have 8 coupled first order ODEs, and I should generate a function like:

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

then do complex_ode(derv）

My questions are:

1. my y is not a list but a matrix, how can i give a corrent output fits into complex_ode()?
2. complex_ode() needs a jacobian, I have no idea how to start constructing one and what type it should be?
3. Where should I put the initial conditions like in the normal ode and time linspace?

this is scipy's complex_ode link: http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.complex_ode.html

Could anyone provide me with more infomation so that I can learn a bit more.

解决方案

I think we can at least point you in the right direction. The optical bloch equation is a problem which is well understood in the scientific community, although not by me :-), so there are already solutions on the internet to this particular problem.

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

However, to address your needs, you spoke of using complex_ode, which I suppose is fine, but I think just plain scipy.integrate.ode will work just fine as well according to their documentation:

 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

You also have the added benefit of an older more established and better documented function. I am surprised you have 8 and not 9 coupled ODE's, but I'm sure you understand this better than I. Yes, you are correct, your function should be of the form ydot = f(t,y), which you call def derv() but you're going to need to make sure your function takes at least two parameters like derv(t,y). If your y is in matrix, no problem! Just "reshape" it in the derv(t,y) function like so:

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

As long as num_rows*num_cols = 8, your number of ODE's you should be fine. Then use the matrix in your computations. When you're all done, just be sure to return a vector and not a matrix like:

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

The Jacobian is not required, but it will likely allow the computation to proceed much more quickly. If you do not know how to compute this you may want to consult wikipedia or a calculus text book. It's pretty simple, but can be time consuming.

As far as initial conditions, you should probably already know what those should be, whether it's complex or real valued. As long as you select values that are within reason, it shouldn't matter much.

这篇关于复杂的ODE系统的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！