问题描述
让我们说我们有三个复矩阵和一个具有这些矩阵的耦合微分方程组.
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()
我基本上有两个问题:
-
如何减少代码?通过减少，我的意思是，有一种方法可以通过不单独写下所有组件，而是在解决系统问题的同时处理矩阵来实现. ODE?
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?
不是相对于 t (代码的最后两行)绘制矩阵元素，而是如何绘制特征值(绝对值)(例如，矩阵A的特征值的绝对值是否是t)的函数?
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)?
推荐答案
今年早些时候，我为scipy.integrate.odeint创建了一个包装器，从而可以轻松地求解复杂的数组微分方程:"> https://github.com/WarrenWeckesser/odeintw
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
您可以使用git检出整个软件包，并使用脚本setup.py安装它，或者可以获取一个基本文件 _odeintw.py ，将其重命名为odeintw.py，然后将其复制到您需要的任何位置.以下脚本使用功能odeintw.odeintw来解决您的系统.它通过将三个矩阵A，B和C堆叠到形状为(3，2，2)的三维数组M中来使用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).
您可以使用numpy.linalg.eigvals计算A(t)的特征值.该脚本显示了一个示例和一个图.特征值很复杂，因此您可能必须进行一些实验才能找到一种可视化它们的好方法.
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()
以下是脚本生成的图:
这篇关于在PYTHON中求解矩阵耦合微分方程时，如何绘制特征值?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！