Поиск «правильного» решения при использовании псевдообратного - PullRequest
Новая
       12

Поиск «правильного» решения при использовании псевдообратного

0 голосов
/ 18 октября 2018

У меня есть система обыкновенных дифференциальных уравнений 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)

этот код дает *Напротив, 1031 * 

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])

мы можем вычислить один и тот же объект во временной области, распространяя ODE и интегрируя в конечную точку (код не показан здесь для компактности), что дает: 

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        ]])

очевидно, что существует значительное несоответствие между результатами.

...