python画微分方程图像,在PYTHON中求解矩阵耦合微分方程时，如何绘制特征值？

Lets say we have three complex matrices and a system of coupled differential equations with these matrices.

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

I have basically two questions:

How to reduce my code? By reduce I meant, is there a way to do this by not writing down all the components separately ,but handling with the matrices while solving the system of ODE?

Instead of plotting elements of the matrices with respect to t (the last 2 lines of my code), how can I plot Eigenvalues (absolute values) (lets say, the abs of eigenvalues of matrix A as a function of t)?

解决方案

Earlier this year I created a wrapper for scipy.integrate.odeint that makes it easy to solve complex array differential equations: https://github.com/WarrenWeckesser/odeintw

You can check out the whole package using git and install it using the script setup.py, or you can grab the one essential file, _odeintw.py, rename it to odeintw.py, and copy it to wherever you need it. The following script uses the function odeintw.odeintw to solve your system. It uses odeintw by stacking your three matrices A, B and C into a three-dimensional array M with shape (3, 2, 2).

You can use numpy.linalg.eigvals to compute the eigenvalues of A(t). The script shows an example and a plot. The eigenvalues are complex, so you might have to experiment a bit to find a nice way to visualize them.

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

Here are the plots generated by the script: