在 PYTHON 中求解矩阵耦合微分方程时如何绘制特征值? [英] How to plot the Eigenvalues when solving matrix coupled differential equations in 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()

我基本上有两个问题:

1. 如何减少我的代码?减少我的意思是，有没有办法通过不单独写下所有组件，而是在解决系统的同时处理矩阵来做到这一点ODE?

1. 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(我的代码的最后 2 行)绘制矩阵的元素，我如何绘制特征值(绝对值)(可以说，矩阵 A 的特征值的 abs 作为 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 来解决您的系统问题.它使用 odeintw 通过将您的三个矩阵 A、B 和 C 堆叠成一个三维数组 M 形状为 (3, 2, 2).

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 中求解矩阵耦合微分方程时如何绘制特征值?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！