Линейное уравнение исключения Гаусса в C++

У меня довольно простой вопрос по математике, но хитрость в том, что он мне нужен на C++. В данный момент я следую псевдокоду, указанному в Википедии. Вот моя попытка:

createMatrixForAllSolutions(*this);
std::cout << equationMatrix.to_string() << endl;
bool solved = false;
int rows = equationMatrix.getRows();
int cols = equationMatrix.getCols();
int i = 0;
int j = 0;
int maxi = 0;
double current = 0;
double eqnValue = 0;
double solValue = 0;
std::vector<char> reversedVars;
int sum = 0;
int tempValue;
int tempRHS;
int newValue;
int neRHS;

while (i < rows && j < cols) {
    maxi = i;
    for (int k = i + 1; k < rows; k++) {
        if (abs(equationMatrix.get_element(k, j)) > abs(equationMatrix.get_element(maxi, j)))
            maxi = k;
    }
    if (equationMatrix.get_element(maxi, j) != 0) {
        current = equationMatrix.get_element(i, j);
        for (int x = 0; x < cols; x++) {
            tempValue = equationMatrix.get_element(i, x);
            newValue = equationMatrix.get_element(maxi, x);
            equationMatrix.set_element(i, x, newValue/current);
            equationMatrix.set_element(maxi, x, tempValue);
        }
        tempRHS = solutionMatrix.get_element(i, 0);
        neRHS = solutionMatrix.get_element(maxi, 0);
        solutionMatrix.set_element(i, 0, neRHS/current);
        solutionMatrix.set_element(maxi, 0, tempRHS);
        //SWAP rows i and maxi
        //SWAP RHS i and maxi
        //DIVIDE each entry in row i by current
        //DIVIDE RHS i by current
        for (int u = i + 1; u < rows; u++) {
            eqnValue = equationMatrix.get_element(u, j) - equationMatrix.get_element(i, j) * equationMatrix.get_element(u, j);
            std::cout << "Equation Value: " << eqnValue << endl;
            equationMatrix.set_element(u, j, eqnValue);
            solValue = solutionMatrix.get_element(u, 0) - solutionMatrix.get_element(i, 0) * solutionMatrix.get_element(u, 0);
            std::cout << "Solution Value: " << solValue << endl;
            solutionMatrix.set_element(u, 0, solValue);
        }
        i++;
    }
    j++;
}

И псевдокод, который я использую, взят из Википедии:

i := 1
j := 1
while (i ≤ m and j ≤ n) do
  Find pivot in column j, starting in row i:
  maxi := i
  for k := i+1 to m do
    if abs(A[k,j]) > abs(A[maxi,j]) then
      maxi := k
    end if
  end for
  if A[maxi,j] ≠ 0 then
    swap rows i and maxi, but do not change the value of i
    Now A[i,j] will contain the old value of A[maxi,j].
    divide each entry in row i by A[i,j]
    Now A[i,j] will have the value 1.
    for u := i+1 to m do
      subtract A[u,j] * row i from row u
      Now A[u,j] will be 0, since A[u,j] - A[i,j] * A[u,j] = A[u,j] - 1 * A[u,j] = 0.
    end for
    i := i + 1
  end if
  j := j + 1
 end while

Я сделал все возможное, чтобы соответствовать этому до сих пор, но если кто-нибудь сможет понять, почему моя лань не работает, это было бы прекрасно. Спасибо!


c++ gaussian linear-algebra
person Brandon    schedule 09.10.2011    source источник
comment
Что происходит с вашим кодом? Он компилируется? Это крах? Дает неверный ответ?   -  person David Alber    schedule 09.10.2011
comment
Он дает неправильные ответы... На самом деле очень большие числа.   -  person Brandon    schedule 09.10.2011

Ответы (3)


arrow_upward
1
arrow_downward

Вот одна проблема: вы объявляете tempValue, tempRHS, newValue и neRHS как целые. Даже если ваша матрица начинается со всех целочисленных значений, они не останутся такими долго, как только вы перейдете к исключению. Все они должны быть объявлены двойными - как целые, вы будете постоянно выбрасывать дробные части.

person dewtell    schedule 09.10.2011

arrow_upward
0
arrow_downward

Я не просматривал ваш код построчно, но наиболее вероятная проблема заключается в разнице в соглашениях об индексации массивов между вашей реализацией C++ и алгоритмом Википедии. C++ использует массивы с отсчетом от 0 (т. е. индекс первого элемента массива равен 0), в то время как алгоритм Википедии отсчитывается от 1 (т. е. индекс первого элемента массива равен 1). Вероятно, вы где-то пропустили какое-то преобразование.

Предполагая, что вы понимаете, что пытается сделать алгоритм, лучше всего будет отказаться от кода и начать заново, исходя из вашего понимания алгоритма и C++.

Если у вас возникли трудности с пониманием алгоритма, вы можете просмотреть копию числовых рецептов в C (глава 2 будет наиболее полезной в этом отношении). Поскольку C и C++ являются языками массивов, основанными на 0, необходимые вам изменения должны быть относительно незначительными по сравнению с использованием версии Википедии в качестве основы для вашей реализации.

person andand    schedule 09.10.2011

arrow_upward
0
arrow_downward

Вы не делите на правильный опорный элемент. На каждом шаге алгоритма опорным элементом является A(maxi,j). Однако в своем коде вы пытаетесь объединить два шага замены и деления на опорный элемент. В результате вы говорите

  current = equationMatrix.get_element(i,j)

Это следует читать

  current = equationMatrix.get_element(maxi,j)

Надеюсь, это поможет. Могут быть еще баги. Это просто первое, что я увидел. При отладке может быть полезно распечатать текущую матрицу на каждом шаге алгоритма. Когда ваше исключение Гаусса работает правильно, ваш алгоритм создаст верхнюю треугольную матрицу с единицами по диагонали. (Подробнее см. на странице Википедии.)

Я надеюсь, что этот код предназначен для образовательных целей. Существует множество прекрасных библиотек линейной алгебры (например, LAPACK). Если вам нужно решить линейные системы, я настоятельно рекомендую использовать одну из этих замечательных библиотек, а не пытаться создать свою собственную.

person codehippo    schedule 09.10.2011