сложные системы ODE в scipy - PullRequest
Новая
       13

сложные системы ODE в scipy

3 голосов
/ 26 февраля 2012

У меня проблемы с совмещением оптического уравнения Блоха, представляющего собой систему ОДУ первого порядка со сложными значениями.Я обнаружил, что 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. Куда мне поместить начальные условия, как в обычной оде и времениlinspace?

это ссылка scipy complex_ode: http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.complex_ode.html

Может ли кто-нибудь предоставить мне больше информации, чтобы я мог узнать немного больше.

1 Ответ

5 голосов
/ 26 февраля 2012

Я думаю, мы можем, по крайней мере, указать вам правильное направление. Оптический уравнение Блоха является проблемой, которая хорошо понимается в научной сообщество, хотя и не мной :-), так что уже есть решения в интернете к этой конкретной проблеме.

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

Якобиан не требуется, но он, вероятно, позволит продолжить вычисления намного быстрее Если вы не знаете, как рассчитать это, вы можете проконсультироваться Википедия или учебник по исчислению. Это довольно просто, но может занять много времени.

Что касается начальных условий, вы, вероятно, уже должны знать, что они должны будь то сложный или реальный. Пока вы выбираете значения, которые в пределах разумного это не должно иметь большого значения.

...