Как построить собственные значения при решении связанных матричных дифференциальных уравнений в PYTHON?

Допустим, у нас есть три комплексные матрицы и система связанных дифференциальных уравнений с этими матрицами.

import numpy, scipy
from numpy import (real,imag,matrix,linspace,array)
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def system(x,t):

    a1= x[0];a3= x[1];a5= x[2];a7= x[3];
    a2= x[4];a4= x[5];a6= x[6];a8= x[7];
    b1= x[8];b3= x[9];b5= x[10];b7= x[11];
    b2= x[12];b4= x[13];b6= x[14];b8= x[15];
    c1= x[16];c3= x[17];c5= x[18];c7= x[19];
    c2= x[20];c4= x[21];c6= x[22];c8= x[23];

    A= matrix([ [a1+1j*a2,a3+1j*a4],[a5+1j*a6,a7+1j*a8] ])  
    B= matrix([ [b1+1j*b2,b3+1j*b4],[b5+1j*b6,b7+1j*b8] ])
    C= matrix([ [c1+1j*c2,c3+1j*c4],[c5+1j*c6,c7+1j*c8] ])

    dA_dt= A*C+B*C
    dB_dt= B*C
    dC_dt= C

    list_A_real= [dA_dt[0,0].real,dA_dt[0,1].real,dA_dt[1,0].real,dA_dt[1,1].real]
    list_A_imaginary= [dA_dt[0,0].imag,dA_dt[0,1].imag,dA_dt[1,0].imag,dA_dt[1,1].imag]

    list_B_real= [dB_dt[0,0].real,dB_dt[0,1].real,dB_dt[1,0].real,dB_dt[1,1].real]
    list_B_imaginary= [dB_dt[0,0].imag,dB_dt[0,1].imag,dB_dt[1,0].imag,dB_dt[1,1].imag]

    list_C_real= [dC_dt[0,0].real,dC_dt[0,1].real,dC_dt[1,0].real,dC_dt[1,1].real]
    list_C_imaginary= [dC_dt[0,0].imag,dC_dt[0,1].imag,dC_dt[1,0].imag,dC_dt[1,1].imag]

    return list_A_real+list_A_imaginary+list_B_real+list_B_imaginary+list_C_real+list_C_imaginary



t= linspace(0,1.5,1000)
A_initial= [1,2,2.3,4.3,2.1,5.2,2.13,3.43]
B_initial= [7,2.7,1.23,3.3,3.1,5.12,1.13,3]
C_initial= [0.5,0.9,0.63,0.43,0.21,0.5,0.11,0.3]
x_initial= array( A_initial+B_initial+C_initial )
x= odeint(system,x_initial,t)

plt.plot(t,x[:,0])
plt.show()

У меня в основном два вопроса:

  1. Как уменьшить мой код? Под сокращением я имел в виду, есть ли способ сделать это, не записывая все компоненты по отдельности, а работая с матрицами при решении системы ODE?

  2. Вместо того, чтобы строить элементы матриц относительно t (последние 2 строки моего кода), как я могу построить собственные значения (абсолютные значения) (скажем, абс собственных значений матрицы A как функция от t)?


python numpy scipy matrix odeint
person string    schedule 27.10.2014    source источник
comment
Вы уверены, что ваш код правильный? Строка a1= x[0];a2= x[1];a3= x[2];a4= x[3];... в сочетании с A= matrix([ [a1+1j*a2,a3+1j*a4],... говорит о том, что реальная и мнимая части A чередуются в x, но когда вы собираете возвращаемое значение system(x,t), вы перечисляете все реальные части dA_dt, за которыми следуют все мнимые части. Вы хотели их перемежать?   -  person Warren Weckesser    schedule 27.10.2014
comment
дифференциальные уравнения: dA (t) / dt = A (t) * C (t) + B (t) * C (t); дБ (t) / dt = B (t) * C (t) и dC (t) / dt = C (t). Каждый элемент (действительная и мнимая части каждого элемента) любой матрицы является функцией t.   -  person string    schedule 27.10.2014
comment
Кстати, представленный мной код отлично работает для меня, за исключением двух упомянутых мной вопросов.   -  person string    schedule 27.10.2014
comment
Я понимаю уравнения; у меня вопрос как реальный, так и воображаемый. значения в A, B и C хранятся в x. Когда вы распаковываете x в начале system(x,t), ясно, что, например, x[1] - это изображение. часть A[0,0]. Таким образом, второй элемент списка, возвращаемый system(x,t), должен быть производной по времени от мнимой части A[0,0]. То есть должно быть dA_dt[0,0].imag. Вы возвращаете list_A_real + ..., что означает, что вы возвращаете list_A_real[1] как производную по времени от мнимой части A[0,0]. Но list_A_real[1] это dA_dt[0,1].real, а не dA_dt[0,0].imag.   -  person Warren Weckesser    schedule 27.10.2014
comment
@WarrenWeckesser, понятно! Итак, ваша точка зрения заключается в том, что я неправильно поддерживаю порядок, возвращая результат ODE. Но что, если я определю свою новую / возвращающуюся матрицу A как A = matrix ([[x [:, 0] + 1j x [:, 4], x [:, 1] + 1j x [ :, 5]], [x [:, 2] + 1j x [:, 6], x [:, 3] + 1j x [:, 7]]])? Это работает, или я должен придерживаться своего первого приказа о назначении?   -  person string    schedule 27.10.2014
comment
Да, вот в чем суть моего вопроса. Предположим, вы вызываете возвращаемое значение dxdt. Вы можете принять любое соглашение для присвоения компонентов A, B и C x, но вы должны использовать то же соглашение для присвоения производных по времени этих компонентов dxdt.   -  person Warren Weckesser    schedule 27.10.2014
comment
@WarrenWeckesser, приятно это знать, позвольте мне отредактировать свой вопрос, чтобы исправить указанную вами проблему.   -  person string    schedule 27.10.2014

Ответы (1)


arrow_upward
4
arrow_downward

Ранее в этом году я создал оболочку для scipy.integrate.odeint, которая упрощает решение сложных дифференциальных уравнений массива: https://github.com/WarrenWeckesser/odeintw

Вы можете проверить весь пакет с помощью git и установить его с помощью сценария setup.py, или вы можете получить один важный файл, _ 3_, переименуйте его в odeintw.py и скопируйте в нужное место. Следующий сценарий использует функцию odeintw.odeintw для решения вашей системы. Он использует odeintw, складывая ваши три матрицы A, B и C в трехмерный массив M с формой (3, 2, 2).

Вы можете использовать numpy.linalg.eigvals для вычисления собственных значений A(t). В сценарии показан пример и сюжет. Собственные значения сложны, поэтому вам, возможно, придется немного поэкспериментировать, чтобы найти хороший способ их визуализировать.

import numpy as np
from numpy import linspace, array
from odeintw import odeintw
import matplotlib.pyplot as plt


def system(M, t):
    A, B, C = M
    dA_dt = A.dot(C) + B.dot(C)
    dB_dt = B.dot(C)
    dC_dt = C
    return array([dA_dt, dB_dt, dC_dt])


t = np.linspace(0, 1.5, 1000)

#A_initial= [1, 2, 2.3, 4.3, 2.1, 5.2, 2.13, 3.43]
A_initial = np.array([[1 + 2.1j, 2 + 5.2j], [2.3 + 2.13j, 4.3 + 3.43j]])

# B_initial= [7, 2.7, 1.23, 3.3, 3.1, 5.12, 1.13, 3]
B_initial = np.array([[7 + 3.1j, 2.7 + 5.12j], [1.23 + 1.13j, 3.3 + 3j]])

# C_initial= [0.5, 0.9, 0.63, 0.43, 0.21, 0.5, 0.11, 0.3]
C_initial = np.array([[0.5 + 0.21j, 0.9 + 0.5j], [0.63 + 0.11j, 0.43 + 0.3j]])

M_initial = np.array([A_initial, B_initial, C_initial])
sol = odeintw(system, M_initial, t)

A = sol[:, 0, :, :]
B = sol[:, 1, :, :]
C = sol[:, 2, :, :]

plt.figure(1)
plt.plot(t, A[:, 0, 0].real, label='A(t)[0,0].real')
plt.plot(t, A[:, 0, 0].imag, label='A(t)[0,0].imag')
plt.legend(loc='best')
plt.grid(True)
plt.xlabel('t')

A_evals = np.linalg.eigvals(A)

plt.figure(2)
plt.plot(t, A_evals[:,0].real, 'b.', markersize=3, mec='b')
plt.plot(t, A_evals[:,0].imag, 'r.', markersize=3, mec='r')
plt.plot(t, A_evals[:,1].real, 'b.', markersize=3, mec='b')
plt.plot(t, A_evals[:,1].imag, 'r.', markersize=3, mec='r')
plt.ylim(-200, 1200)
plt.grid(True)
plt.title('Real and imaginary parts of the eigenvalues of A(t)')
plt.xlabel('t')
plt.show()

Вот графики, созданные скриптом:

‹code› A [0,0] ‹/code›

собственные значения A

person Warren Weckesser    schedule 27.10.2014
comment
Я скопировал файл и назвал его odeintw.py в ту же папку. Теперь, когда я запускаю приведенный выше код, который вы записали, я получаю следующую ошибку: Traceback (последний вызов последним): File / home / ....., строка 27, в ‹module› sol = odeintw (system, M_initial, t) Файл /home/...../odeintw.py, строка 180, в odeintw y0 = y0.astype (np.complex128, copy = False) TypeError: astype () не принимает аргументов ключевого слова - person string; 27.10.2014
comment
Какую версию numpy вы используете? Я предполагаю, что аргумент copy метода astype был добавлен в более позднюю версию numpy. Вы можете удалить этот аргумент из вызова astype(...), и код все равно будет работать нормально. В файле есть два вызова astype - измените их оба. (Я добавлю make odeintw более терпимым к более старым версиям numpy в свой список дел. :) - person Warren Weckesser; 27.10.2014
comment
Аргумент copy astype был добавлен в numpy 1.7. - person Warren Weckesser; 27.10.2014
comment
@string: я обновил файл _odeintw.py на github для обработки numpy до версии 1.7. У меня это сработало с numpy 1.6.2. - person Warren Weckesser; 27.10.2014
comment
Спасибо за вашу тяжелую работу, но, к сожалению, теперь я получаю еще одну ошибку. Моя версия numpy - 1.6.1, и теперь я использую ваш обновленный пакет _odeint.w. Ошибка: Traceback (последний вызов последним): файл /home/name.py, строка 40, в ‹module› A_evals = np.linalg.eigval (A) Файл /usr/lib/python2.7/dist-packages /numpy/linalg/linalg.py, строка 765, в eigval'ах _assertRank2 (a) Файл /usr/lib/python2.7/dist-packages/numpy/linalg/linalg.py, строка 155, в _assertRank2 двухмерном '% len (a.shape) LinAlgError: задан трехмерный массив. Массив должен быть двухмерным - person string; 27.10.2014
comment
Неважно, все эти ошибки связаны со старой версией numpy, я только что обновился до версии 1.9.0, и все работает без сбоев. Большое спасибо за ваши предложения, я ценю это. Я собираюсь принять ваш ответ. - person string; 28.10.2014
comment
Ах, я забыл, что обобщенное поведение np.linalg.eigvals также является относительно новой функцией numpy. Обновление имеет смысл; Я рад слышать, что это работает. - person Warren Weckesser; 28.10.2014