使用伪逆时寻找“正确”解

2024-10-01 00:34:37 发布

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

我有一个常微分方程组x'(t) = Lx(t),我想在x(t=0)考虑初始条件时找到频率空间行为。我可以用两种方法来做:

  1. 在时域传播方程,然后做(单侧)傅里叶变换;或者
  2. 对运动方程做拉普拉斯变换,得到线性方程(1j * w * ident - L) X(w) = x(t=0)并直接求解(注意,我限制拉普拉斯变量为纯虚数)

这可以通过以下方式轻松完成:

def resolvant(w,L):

    R = np.linalg.inv(1j * ident * w - L)

    return R

这适用于1j * ident * w-L是非单数的所有点。现在我的矩阵L有一个零特征值,其余的特征值是复数。因此,分辨率仅在点w=0处是单数的

此时,我可以用伪逆函数替换逆函数:

def resolvant(w,L):
    if w !=0:
        R = np.linalg.inv(1j * ident * w - L)
    else:
        R = np.linalg.pinv(- L)
    return R

在许多情况下,这是可行的,但是在一些重要的情况下,伪逆方法不同于通过时域方法发现的方法。有没有一种约束伪逆的方法来给出w=0的正确行为

编辑: 最小工作示例

import numpy as np

#construct matrix on for the RHS of equation
L = np.array([[0,-1.j,0,1.j,0,0,0,0,0],[-1.j,-0.5,0,0,1.j,0,0,0,0],[0,0,0,0,0,1.j,0,0,0],[1.j,0,0,-0.5,-1.j,0,0,0,0],[0,1.j,0,-1.j,-1.,0,0,0,0],[ 0,0,1.j,0,0,-0.5,0,0,0],[0,0,0,0,0,0,0,-1.j,0],[ 0,0,0,0,0,0,-1.j,-0.5,0],[0,0,0,0,1,0,0,0,0]])

#initial vector is given as:

x0=np.zeros(L.shape[0])
x0[0] = 1

#calculate the pseudo inverse and resultant vector:
R = np.linalg.pinv(-L)

res = np.dot(R, x0)
print(res)

此代码生成

array([  5.00000000e-01 +1.04874666e-17j,
     1.74736751e-16 -3.33333333e-01j,
     1.23748154e-16 +6.69892995e-17j,
     6.13265296e-17 +3.33333333e-01j,
     3.33333333e-01 +2.70701098e-17j,
     1.24640823e-17 +9.26971489e-17j,
     0.00000000e+00 +0.00000000e+00j,
     0.00000000e+00 +0.00000000e+00j,   0.00000000e+00 +0.00000000e+00j])

相反,我们可以通过传播常微分方程并积分到终点来计算时域中的同一对象(此处未显示紧凑性代码),这将给出:

array([[  1.25024870+0.j        ],
   [  0.00000000-0.49999962j],
   [  0.00000000+0.j        ],
   [  0.00000000+0.49999962j],
   [  1.00000113+0.j        ],
   [  0.00000000+0.j        ],
   [  0.00000000+0.j        ],
   [  0.00000000+0.j        ],
   [ 47.75025018+0.j        ]])

显然,结果之间有相当大的差距


Tags: 方法returndefnparray方程单数特征值