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

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

1 голос
/ 09 октября 2011

У меня довольно простой математический вопрос, но дело в том, что он мне нужен в 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

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

Ответы [ 3 ]

1 голос
/ 09 октября 2011

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

0 голосов
/ 09 октября 2011

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

  current = equationMatrix.get_element(i,j)

Это должно читаться как 

  current = equationMatrix.get_element(maxi,j)

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

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

0 голосов
/ 09 октября 2011

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

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

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

...