你好，我有一个问题，解决以下动态Riccati方程：
这是我写的代码：
import numpy as np
from scipy.integrate import odeint
import networkx as nx
import matplotlib.pyplot as plt
def deriv(sigma, t, A, B, Md):
sigma = np.reshape(sigma, s)
return (-np.matmul(A, sigma) - np.matmul(sigma, A.transpose()) + B - 2 * np.matmul(sigma, Md, sigma)).flatten()
n = 100
alpha = 5
beta = 5
G = nx.barabasi_albert_graph(n=n, m=2, seed=10374196, initial_graph=None)
A = nx.adjacency_matrix(G).toarray() # Convert to dense matrix
B = alpha * np.identity(n)
Md = np.identity(n)
Md[1, 1] = beta
s = (n, n)
sigma0 = np.ones(n**2)
t = np.linspace(0, 100, 101)
sol = odeint(deriv, sigma0, t, args=(A, B, Md))
sol = np.reshape(sol, (len(t),) + s)
Y = np.linalg.eigvals(sol[-1])
Y.sort()
X = [i for i in range(n)]
plt.plot(X, Y)
plt.xlabel("eigenvalue rank")
plt.ylabel("eigenvalue")
plt.show()我在测试稳态解，得到一个零矩阵，尽管我知道这是不可能的，因为这些是，黎卡提方程的期望解，稳态解的特征值。enter image description here有人知道我哪里犯了错误吗？
1条答案
按热度按时间w80xi6nr1#
复制并运行您的确切代码，我无法得到一个解决方案，由于溢出错误。我修改了你的代码，使用更现代的
scipy.integrate.solve_ivpODE求解器，并修改了你的deriv函数，使用一些numpy的简写。通过这些变化，我得到了所有时间的非零解矩阵。我不确定您的问题是否与使用odeint（过时的scipy ODE求解器）或其他有关。