Решение систем линейных уравнений с использованием собственных - PullRequest
Новая
       21

Решение систем линейных уравнений с использованием собственных

0 голосов
/ 04 сентября 2018

В настоящее время я работаю над моделированием жидкости в C ++, и часть алгоритма состоит в том, чтобы решить разреженную систему линейных уравнений. Люди рекомендовали использовать для этого библиотеку Eigen. Я решил проверить это с помощью этой короткой программы, которую я написал: 

#include <Eigen/SparseCholesky>

#include <vector>
#include <iostream>

int main() {
    std::vector<Eigen::Triplet<double>> triplets;
    triplets.push_back(Eigen::Triplet<double>(0, 0, 1));
    triplets.push_back(Eigen::Triplet<double>(0, 1, -2));
    triplets.push_back(Eigen::Triplet<double>(1, 0, 3));
    triplets.push_back(Eigen::Triplet<double>(1, 1, -2));

    Eigen::SparseMatrix<double> A(2, 2);
    A.setFromTriplets(triplets.begin(), triplets.end());

    Eigen::VectorXd b(2);
    b[0] = -2;
    b[1] = 2;

    Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(A);
    Eigen::VectorXd x = chol.solve(b);

    std::cout << x[0] << ' ' << x[1] << std::endl;

    system("pause");
}

Это дает ему два уравнения:

x - 2y = -2

3x - 2y = 2

Правильное решение:

x = 2

y = 2

Но проблема в том, что при запуске программы она выдает: 0,181818 -0,727273

Что совершенно неправильно! Я отлаживал это часами, но это очень короткая программа, и я точно следую учебному пособию на сайте Eigen. Кто-нибудь знает, что является причиной этой проблемы?

P.S. Я знаю, что классы, которые я использую, предназначены для разреженных матриц, но единственное различие между ними и обычными классами Matrix заключается в способе хранения элементов.

Ответы [ 2 ]

0 голосов
/ 09 сентября 2018

Вы должны использовать решатель, способный обрабатывать несимметричные разреженные матрицы. Другой возможный подход заключается в поиске решения не исходной системы [A] x = b, а [A] T * [A] x = [A] T * b, где [A] T обозначает транспонирование [A] , Матрица последней системы является симметричной и положительно определенной (при условии, что [A] неособо). Единственным недостатком может быть тот факт, что [A] T [A] может быть довольно плохо обусловленным, если оригинал [A] не является «хорошим» в этом смысле. Просто пример программного обеспечения, предназначенного для решения таких проблем: http://members.ozemail.com.au/~comecau/CMA_LS_Sparse.htm

0 голосов
/ 04 сентября 2018

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

x + 3y = -2
3x -2y = 2

Как вы заметили, для несимметричных квадратичных задач вам нужно использовать прямой решатель на основе LU или BICGSTAB в мире итерационных решателей. Все это обобщено в документ .

...