Сципи неправильно решает собственную задачу матрицы - PullRequest
Новая
       22

Сципи неправильно решает собственную задачу матрицы

0 голосов
/ 15 ноября 2018

Я сталкиваюсь с очень странным поведением SciPy при решении собственной проблемы сингулярной матрицы, а именно, вычисленные собственные значения неверны, если я генерирую матрицу по некоторым функциям (матрица der в коде ниже). Однако, если я вручную введу матрицу в (der2), диагонализация, похоже, даст правильные результаты. Это также можно проверить путем вычитания обеих матриц, что делается в приведенном ниже коде.

код 

import numpy as np
import scipy as sp
from scipy.linalg import eigvals

def cbar(k, n):
    """
    cbar function for coefficients
    """
    if k==0 or k==n-1:
        return np.float128(2.)
    else:
        return np.float128(1.)    

def ChebCollDer(x):
    """
    ChebCollDer  Chebyshev collocation differentiation matrix.
    """
    xx=np.array(x)  
    n=xx.size
    d=np.zeros((n,n))
    for i in range(n):
        for k in range(n):
            if i!=k:
                d[i,k]=cbar(i, n)*np.float128(sp.power(-1., i+k))/(xx[i]-xx[k])/cbar(k,n)
    for i in range(n):
        tmp=-sp.sum(d[i,:])
        d[i,i]=tmp

    return d
nn=5
xx=(np.cos(sp.pi*np.linspace(0,1.0,nn)))/2.
der=ChebCollDer(xx)
print eigvals(der)
der2=[[ 11.0 ,-13.656854,  4.0 ,-2.3431458 , 1.0],
 [ 3.4142136, -1.4142136, -2.8284271,  1.4142136, -0.58578644],
 [-1.0,  2.8284271,  1.110223e-16, -2.8284271,  1.0]
 ,[0.58578644, -1.4142136,  2.8284271,  1.4142136, -3.4142136],
 [-1.0,2.3431458,-4.0 ,13.656854,-11.0]]
print eigvals(der2)

print der-der2

и результаты: Собственные значения матрицы der: 

[ 0.00389434+0.00282825j  0.00389434-0.00282825j -0.00148641+0.00457958j -0.00148641-0.00457958j -0.00481586+0.j        ]

Собственные значения der2: 

[  9.71161644e-02+0.j          -9.71161490e-02+0.j          -3.08158279e-08+0.j 7.69629619e-09+0.09711159j   7.69629619e-09-0.09711159j]

видите, что der2 имеет одно собственное значение, которое численно равно нулю, как и должно быть, поскольку матрица der имеет нулевой собственный вектор, который просто [1,1,1,1,1] Самый большой элемент der-der2 имеет порядок 10E-08. Я подозреваю, что есть какая-то проблема с преобразованием типов, но я не знаю, откуда она.

1 Ответ

0 голосов
/ 25 декабря 2018

Это не полный ответ, но у меня недостаточно репутации, чтобы комментировать.

Возможно, вам больше повезет, если вы разместите на странице математики больше, так как я думаю, что это проблема численного анализа.

Если я сделаю s,v = sp.linalg.eig(der), то все u, похоже, удовлетворяют уравнению собственных значений в том смысле, что der@u - s*u близко к нулю.

Еще одна интересная вещь заключается в том, что собственные значения derвсе имеют одинаковый модуль и имеют аргументы, одинаково расположенные вокруг окружности этого радиуса.Как и масштабированные 5-е корни -1.

Фактически, если вы возьмете и одинаково взвешенную линейную комбинацию собственных векторов из вывода scipy (то есть умножите вправо v на вектор всех единиц), вы получитеблизко к постоянному вектору, соответствующему нулевому собственному значению.

Тем не менее, я не очень уверен, что происходит.

Кроме того, ваша функция ChebCollDer будет выводить регулярный символ массиваfloat вместо float128, поскольку вы не инициализируете его на float128.Это означает, что даже если вы выполняете некоторые вычисления в float128, они преобразуются в плавающее число по умолчанию, прежде чем они будут сохранены в d.

...