при определенном решении системы python - PullRequest
Новая
       60

при определенном решении системы python

0 голосов
/ 22 февраля 2019

Цель: вычислить вектор из определенной линейной системы (2x3) Ax = b
Третье уравнение должно быть уравнением единицы (x ^ 2 + y ^ 2 +z ^ 2 = 1).У меня есть правильные матричные коэффициенты, но я не могу получить правильный результат;пытаясь решить Ax = b следующим образом:

Функция возвращает пустое пространство оператора.Затем я устанавливаю матрицу и пытаюсь ее решить. 

from scipy.linalg import qr, null_space, svd
from scipy import transpose, compress

def null(A, eps=1e-17):
    u, s, vh = svd(A)
    padding = max(0,np.shape(A)[1]-np.shape(s)[0])
    null_mask = np.concatenate(((s <= eps), np.ones((padding,),dtype=bool)),axis=0)
    null_space = compress(null_mask, vh, axis=0)
    return transpose(null_space)

У нас есть 3 вершины, которые устанавливают треугольник: 

vh0 = [0., -1., 0.]
vh1 = [-0.03806, -0.98078501, -0.191341]
vh2 = [-0.074658, -0.98078501, -0.18024001]

# normal vector of vh0
normal_vec = [ 0., -0.23760592, 0.]

cap_vec10 = np.subtract(vh1, vh0)
cap_vec20 = np.subtract(vh2, vh0)

a1 = np.array(np.subtract(cap_vec20, cap_vec10))
a2 = np.array(np.dot(-1, capvec10))

# orientation bit of the normal vector
ob = np.sign(np.linalg.det([x_k, x_k1, normal_vec])) 

# normal vector of vertex vh1 that I want to get solving the system
normal_vec1 = [-0.04744975, -0.97674069, -0.209108]

Lm   = np.dot(np.subtract(vh2, vh1), normal_vec1) 
Lm_1 = np.dot(np.subtract(vh0, vh1), normal_vec1) 

# solving under determined system
A = np.array([a1, a2]) 
b = np.array([Lm, Lm_1])
x_lstsq = np.linalg.lstsq(A, b)[0]
wanted_norm = np.sqrt(abs(1 - (np.linalg.norm(x_lstsq)*np.linalg.norm(x_lstsq))))

Z = null(A)*wanted_norm 
new_normal_vec = np.add(Z[:, 0], x_lstsq)

if np.sign(np.linalg.det([x_k, x_k1, Z[:, 0]])) != ob:
    new_normal_vec[list(np.abs(x_lstsq)).index(min(np.abs(x_lstsq)))] *= ob

print("should_be:   {}\ncounted_nv:  {}".format(np.round(normal_vec1, 3), np.round(new_normal_vec, 3)))

normal_vec1 - это нужный мне вектор.И для обоих векторов Z * vector == 1.
Коэффициенты в коде: L_m = , <> - скалярное умножение
Как я понял, два уравнения задают линию, а нормализация дает единицусфера.Таким образом, мое решение - пересечение точек линии и сферы единства.Но также не могу понять, как получить оба решения.

Ответы [ 2 ]

0 голосов
/ 25 февраля 2019

Цель состояла в том, чтобы: в основном решить линейную при определенной системе 2x3 с ограничением, что результирующий вектор должен быть нормализован.
То, что я сделал, дает одно из решений (поскольку в трехмерном пространстве имеется от 0 до 2 векторов):

1.Расчетное решение наименьших квадратов: 

x_lstsq = np.linalg.lstsq(A, b)[0]

2.Вычисляем нулевое пространство (ядро оператора A) 

wanted_norm = np.sqrt(abs(1 - (np.linalg.norm(x_lstsq)*np.linalg.norm(x_lstsq))))
Z = null(A)*wanted_norm

3.Вычислили результирующий вектор 

result = np.add(Z[:, 0], x_lstsq)

Таким образом, я получил один из двух векторов, которые не подходили для моего проекта, но были верны для этой конкретной линейной системы.Поэтому мой вопрос был: как заставить второго делать те же шаги (через пространство нулей).

В процессе поиска решения я понял еще один: в основном решить эту линейную систему вручную, используяуравнение нормализации как третье уравнение.Геометрически оба первых двух уравнения задают две плоскости.Пересечение этих плоскостей дает линию.Такое уравнение: x ^ 2 + y ^ 2 + z ^ 2 = 1 устанавливает единичную сферу.Итак, пересечение этой линии и этой сферы дает нам две точки.Следовательно, решение квадратного уравнения одной из координат дает форму от 0 до 2 корней (от 0 до 2 векторов).

0 голосов
/ 22 февраля 2019

Попробуйте использовать псевдообратную функцию np.linalg.pinv: https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.linalg.pinv.html

...