特征值的曲线图混淆了。。(Python、Matplotlib、Numpy、特征值)

2024-04-16 11:31:42 发布

您现在位置:Python中文网/ 问答频道 /正文

请帮助我理解这里出了什么问题,我已经调查了所有可能的问题,我可以想到,但无法解释这一现象

我正在计算N 2x2矩阵的本征值,包含在一个形状数组(N,M,M)中,其中M=2,每个矩阵产生两个本征值,即结果是一个形状数组(N,M),我想绘制所有N个本征值的M图

也许我误解了一些关于本征值或numpy计算它们的方式,但我最终得到的结果是这样的

attempted plot of eigen-values

正如您可以看到蓝色和橙色轨迹(第一和第二本征值)混合在一起,我如何避免这种情况

这是我的密码

class StateSpaceModel:
    def __init__(self, A=None, B=None, C=None, D=None):
        self.A = A
        self.B = B
        self.C = C
        self.D = D

    def passivity_violations(self):
        A, B, C, D = self.A.real, self.B.real, self.C.real, self.D.real

        δ = 0
        ϵ = 0
        tol = 1e-12

        I = np.identity(len(D))

        AΨ = scipy.linalg.block_diag(A, -A.T)
        BΨ = np.block([[B], [-C.T]])
        CΨ = np.block([[C.T], [B]]).T
        DΨ = D + D.T

        def N(δ):
            δ = np.atleast_1d(δ)
            δ = δ.reshape((*δ.shape, 1, 1))

            DΨδ = DΨ - 2 * δ * I
            Nδ = AΨ - BΨ @ np.linalg.inv(DΨδ) @ CΨ

            return np.squeeze(Nδ)

        N0 = N(δ=δ)

        λ = np.linalg.eigvals(N0)
        λ = λ[λ.imag > 0]

        χ = np.sort(λ[np.abs(λ.real / λ.imag) < tol].imag)

        ξk = (np.hstack(([0], χ)) + np.hstack((χ, 2 * χ[[-1]]))) / 2

        Ψ = lambda s: self(s) + np.moveaxis(self(-s).T, -1, 0)

        _, ax = plt.subplots()
        _f = np.linspace(1, 1e6, 2000)
        _eig = np.linalg.eigvals(Ψ(2j * np.pi * _f).real)
        print(_eig.shape)
        _passive = _eig >= 0
        print(_passive.shape)

        # Plot horizontal threshold-line (0).
        ax.plot(_f, np.full(_f.shape, δ), '--', color='black')

        # Plot passive points.
        ax.plot(_f, _eig, 'o')
        # ax.plot(_f[_passive[:, 1]], _eig[_passive[:, 1], 1], 'o', color='green')

        # Plot passivity violation points.
        ax.plot(_f[~_passive[:, 0]], _eig[~_passive[:, 0], 0], 'o', color='red')
        # ax.plot(_f[~_passive[:, 1]], _eig[~_passive[:, 1], 1], 'o', color='deeppink')

        # Plot vertical "δ-crossing" lines.
        for x in χ:
            ax.plot([x / (2 * np.pi)], [0], 'x', color='black', markersize=20)
            # ax.plot([x / (2 * np.pi)], [0], 'x', color='black', markersize=20)

        # # Plot vertical "δ-crossing" lines.
        # for x in χ:
        #     ax0.plot([x / (2 * np.pi)] * 2, [min0, max0])
        #     ax1.plot([x / (2 * np.pi)] * 2, [min1, max1])

        plt.show()


Tags: selfnoneplotdefnppiaxblock