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

Question

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

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. Вместо построения элементов матриц с относительно т (последние 2 строки моего кода), как я могу сюжет собственных значений (абсолютные значения) (позволяет сказать, абс собственных значений матрицы А, как функция от t)?

Answer 1

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

Вы можете проверить весь пакет, используя мерзавца и установить его с помощью сценария setup.py, или вы можете захватить один важный файл, _odeintw.py, переименуйте его в 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() 

Вот участки, полученные с помощью сценария: