Решение алгебраических уравнений методом Гаусса

[m_ax]

Приветствую)

Есть ли у кого готовая реализация решения системы алг. уравнений методом Гаусса?

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

Буду оч признателен))

Спасибо за внимание)   

[m_ax]

Короче, пришлось писать самому)

Пролистал перед тем как писать одну книжонку "Самоучитель Turbo Pascal" О.А. Меженный, 2003. - 336с. Кстати, эта первая книга, с которой  началось моё знакомство с программированием. Очень даже неплохо написана, между прочим и там собственно есть алгоритм и исходники метода гаусса.

В общем, может кому пригодится:

gauss.h

Код

C++ (Qt)
#ifndef GAUSS_H
#define GAUSS_H
 
#include "array1d.h"
#include "array2d.h"
 
enum Status { NO_ERROR = 1,
              SYSTEM_INCOMPATIBLE = -1,
              INVALID_MATRIX = -2 };
 
Status gauss(Array2D<double> &extMatrix, Array1D<double> &x);
 
#endif // GAUSS_H
 

array2d.h

Код

C++ (Qt)
#ifndef ARRAY2D_H
#define ARRAY2D_H
 
 
template <class T>
class Array2D
{
public:
    Array2D(int rows = 1, int columns = 1, const T &initValue = T());
    Array2D(const Array2D<T> &other);
    ~Array2D();
    T &operator() (int row, int column, bool *isValid = 0);
    const T &operator() (int row, int column, bool *isValid = 0) const;
    Array2D &operator=(const Array2D<T> &array);
    void resize(int rows, int columns);
    int columns() const { return _columns; }
    int rows() const { return _rows; }
private:
    T **_arr;
    T *_trash;
    int _rows;
    int _columns;
};
//-----------------------------------------------------------------------------
 
template <class T>
Array2D<T>::Array2D(int rows, int columns, const T &initValue)
    :_rows(rows), _columns(columns)
{
    if (_rows < 1)
        _rows = 1;
    if (_columns < 1)
        _columns = 1;
 
    _arr = new T*[_rows];
    for (int i = 0; i < _rows; ++i)
        _arr[i] = new T[_columns];
 
    for (int i = 0; i < _rows; ++i) {
        for (int j = 0; j < _columns; ++j) {
            _arr[i][j] = initValue;
        }
    }
 
    _trash = new T;
}
//-----------------------------------------------------------------------------
 
template <class T>
Array2D<T>::Array2D(const Array2D<T> &other)
{
    _rows = other.rows();
    _columns = other.columns();
 
    _arr = new T*[_rows];
    for (int i = 0; i < _rows; ++i)
        _arr[i] = new T[_columns];
 
    for (int i = 0; i < _rows; ++i) {
        for (int j = 0; j < _columns; ++j) {
            _arr[i][j] = other(i, j);
        }
    }
 
    _trash = new T;
}
//-----------------------------------------------------------------------------
 
template <class T>
Array2D<T>::~Array2D()
{
    for (int i = 0; i < _rows; ++i)
        delete []_arr[i];
    delete []_arr;
 
    delete _trash;
}
//-----------------------------------------------------------------------------
 
template <class T>
Array2D<T> &Array2D<T>::operator=(const Array2D<T> &array)
{
    if (this == &array) {
        return *this;
    }
    if (_rows != array.rows() || _columns != array.columns()) {
        for (int i = 0; i < _rows; ++i)
            delete []_arr[i];
        delete []_arr;
 
        _rows = array.rows();
        _columns = array.columns();
 
        _arr = new T*[_rows];
        for (int i = 0; i < _rows; ++i)
            _arr[i] = new T[_columns];
    }
 
    for (int i = 0; i < _rows; ++i) {
        for (int j = 0; j < _columns; ++j) {
            _arr[i][j] = array(i, j);
        }
    }
    return *this;
}
//-----------------------------------------------------------------------------
 
template <class T>
T &Array2D<T>::operator ()(int row, int column, bool *isValid)
{
    if (row >= _rows || row < 0) {
        if (isValid)
            *isValid = false;
        return *_trash;
    }
 
    if (column >= _columns || column < 0) {
        if (isValid)
            *isValid = false;
        return *_trash;
    }
 
    if (isValid)
        *isValid = true;
 
    return _arr[row][column];
}
//-----------------------------------------------------------------------------
 
template <class T>
const T &Array2D<T>::operator ()(int row, int column, bool *isValid) const
{
    if (row >= _rows || row < 0) {
        if (isValid)
            *isValid = false;
        return *_trash;
    }
 
    if (column >= _columns || column < 0) {
        if (isValid)
            *isValid = false;
        return *_trash;
    }
 
    if (isValid)
        *isValid = true;
 
    return _arr[row][column];
}
//-----------------------------------------------------------------------------
 
template <class T>
void Array2D<T>::resize(int rows, int columns)
{
    if (rows < 1 || columns < 1)
        return;
 
    if (rows != _rows || columns != _columns) {
        T **newArr = new T*[rows];
        for (int i = 0; i < rows; ++i)
            newArr[i] = new T[columns];
 
        for (int i = 0; i < rows; ++i) {
            for (int j = 0; j < columns; ++j) {
                newArr[i][j] = T();
            }
        }
 
        int r = (_rows > rows) ? rows : _rows;
        int c = (_columns > columns) ? columns : _columns;
 
        for (int i = 0; i < r; ++i) {
            for (int j = 0; j < c; ++j) {
                newArr[i][j] = _arr[i][j];
            }
        }
 
        for (int i = 0; i < _rows; ++i)
            delete []_arr[i];
        delete []_arr;
        _arr = newArr;
 
        _columns = columns;
        _rows = rows;
    }
}
//-----------------------------------------------------------------------------
 
#endif // ARRAY2D_H
 

array1d.h

Код

C++ (Qt)
#ifndef ARRAY1D_H
#define ARRAY1D_H
 
template <class T>
class Array1D
{
public:
    Array1D(int size = 1, const T &initValue = T());
    Array1D(const Array1D<T> &other);
    ~Array1D();
    T &operator() (int i, bool *isValid = 0);
    const T &operator() (int i, bool *isValid = 0) const;
    Array1D &operator=(const Array1D<T> &array);
    void resize(int size);
    int size() const { return _size; }
private:
    T *_arr;
    T *_trash;
    int _size;
};
 
template <class T>
Array1D<T>::Array1D(int size, const T &initValue)
    :_size(size)
{
    if (_size < 1)
        _size = 1;
 
    _arr = new T[_size];
 
    for (int i = 0; i < _size; ++i) {
        _arr[i] = initValue;
    }
 
    _trash = new T;
}
//-----------------------------------------------------------------------------
 
template <class T>
Array1D<T>::Array1D(const Array1D<T> &other)
{
    _size = other.rows();
    _arr = new T[_size];
 
    for (int i = 0; i < _size; ++i)
        _arr[i] = other(i);
 
    _trash = new T;
}
//-----------------------------------------------------------------------------
 
template <class T>
Array1D<T>::~Array1D()
{
    delete []_arr;
 
    delete _trash;
}
//-----------------------------------------------------------------------------
 
template <class T>
Array1D<T> &Array1D<T>::operator=(const Array1D<T> &array)
{
    if (this == &array) {
        return *this;
    }
    if (_size != array.size()) {
        delete []_arr;
 
        _size = array.size();
 
        _arr = new T[_size];
    }
 
    for (int i = 0; i < _size; ++i) {
        _arr[i] = array(i);
    }
    return *this;
}
//-----------------------------------------------------------------------------
 
template <class T>
T &Array1D<T>::operator ()(int i, bool *isValid)
{
    if (i >= _size || i < 0) {
        if (isValid)
            *isValid = false;
        return *_trash;
    }
 
    if (isValid)
        *isValid = true;
 
    return _arr[i];
}
//-----------------------------------------------------------------------------
 
template <class T>
const T &Array1D<T>::operator ()(int i, bool *isValid) const
{
    if (i >= _size || i < 0) {
        if (isValid)
            *isValid = false;
        return *_trash;
    }
 
    if (isValid)
        *isValid = true;
 
    return _arr[i];
}
//-----------------------------------------------------------------------------
 
template <class T>
void Array1D<T>::resize(int size)
{
    if (size < 1)
        return;
 
    if (size != _size) {
        T *newArr = new T[size];
 
        for (int i = 0; i < size; ++i) {
            newArr[i] = T();
        }
 
        int s = (_size > size) ? size : _size;
 
        for (int i = 0; i < s; ++i) {
            newArr[i] = _arr[i];
        }
 
        delete []_arr;
        _arr = newArr;
 
        _size = size;
    }
}
//-----------------------------------------------------------------------------
 
#endif // ARRAY1D_H
 

gauss.cpp

Код

C++ (Qt)
#include <limits>
#include <cmath>
#include "gauss.h"
 
Status gauss(Array2D<double> &extMatrix, Array1D<double> &x) {
 
    if (extMatrix.rows()+1 != extMatrix.columns())
        return INVALID_MATRIX;
 
    const double epsilon = std::numeric_limits<double>::epsilon();
    const int rows = extMatrix.rows();
    const int columns = extMatrix.columns();
 
    for (int i = 0; i < columns-1; ++i) {
        double maxVal = extMatrix(i, i);
        int index = i;
        for (int j = i+1; j < rows; ++j) {
            if (fabs(extMatrix(j, i)) > maxVal) {
                maxVal = extMatrix(j, i);
                index = j;
            }
        }
        if (fabs(maxVal) < epsilon)
            return SYSTEM_INCOMPATIBLE;
 
        if (index != i) {
            for (int j = 0; j < columns; ++j) {
                double x = extMatrix(index, j);
                extMatrix(index, j) = extMatrix(i, j);
                extMatrix(i, j) = x;
            }
        }
        for (int r = i+1; r < rows; ++r) {
            for (int c = columns-1; c >= i; --c) {
                double s = -extMatrix(r, i)/extMatrix(i, i);
                extMatrix(r, c) += (s*extMatrix(i, c));
            }
        }
    }
 
    x.resize(extMatrix.rows());
    for (int i = rows-1; i >= 0; --i) {
        double s = 0.0;
        for (int j = i+1; j < columns-1; ++j) {
            s += x(j)*extMatrix(i, j);
        }
        x(i) = (extMatrix(i, columns-1)-s)/extMatrix(i, i);
    }
 
    return NO_ERROR;
}
 

Пример использования:
main.cpp

Код

C++ (Qt)
#include <iostream>
#include "gauss.h"
 
using namespace std;
 
int main()
{
   Array2D<double> m(3, 4);
   m(0, 0) = 2.0;
   m(0, 1) = -4.0;
   m(0, 2) = 4.0;
   m(0, 3) = 6.0;
 
   m(1, 0) = -5.0;
   m(1, 1) = 12.0;
   m(1, 2) = -14.0;
   m(1, 3) = -35.0;
 
   m(2, 0) = 6.0;
   m(2, 1) = -7.0;
   m(2, 2) = 5.0;
   m(2, 3) = 10.0;
 
   Array1D<double> x(3);
   Status status = gauss(m, x);
   if (status) {
       for (int i = 0; i < m.rows(); ++i) {
           for (int j = 0; j < m.columns(); ++j)
                cout << m(i, j) << "  ";
           cout << x(i) << endl;
       }
   }
    return 0;
}
 

Несколько тонких мест:
В реализации Array2D<T> для индексации элементов перегружены операторы

Код

C++ (Qt)
T &operator() (int row, int column, bool *isValid = 0);
const T &operator() (int row, int column, bool *isValid = 0) const;
 

хотя привычнее казалось бы нужно делать так:

Код

C++ (Qt)
T *operator[] (int i);
const T *operator[] (int i) const {
    return _arr[i]; // здесь нужно быть осторожным и следить за валидностью i
}
 

тогда мы могли бы обращаться к элементам массива так: m[1][2] = 3.4;
Но, было решено отказаться от ентого и сделать так, как есть сейчас) Почему?

Далее, касательно функции gauss:
первый аргумент - это расширенная матрица системы. (столбцов на один больше чем строк). Она передаётся не как константная ссылка, поскольку она в процессе решения приводится к треугольному виду. Можно было, конечно, внутри gauss создать временный объект:

Код

C++ (Qt)
Array2D<double> m = extMatrix;
 

но это слишком накладно..
Последний аргумент Array1D<double> &x - это результат, в том случае, если функция вернёт статус NO_ERROR

Остальное, думаю прозрачно)   

[Авварон]

передавать bool * ok в оператор индексации смысла 0. Чтобы не было ошибок, надо юзать такой код:

Код:

const T &operator[] (int i) const {
if (i > 0 && i < size())
    return _arr[i];
else
    return T();
}

А вот в ф-ию рассчета гаусса булку передавать надо, и возвращать она должна не эррор, а НОВЫЙ вектор-столбец, который является решением ур-я.

[m_ax]

передавать bool * ok в оператор индексации смысла 0. Чтобы не было ошибок, надо юзать такой код:
Код:

Код

C++ (Qt)
const T &operator[] (int i) const {
if (i > 0 && i < size())
    return _arr[i];
else
    return T();
}
 

Да, Вы правы, если речь идёт о одномерном массиве. Здесь проверка реализуется элементарно.
Кстати, зачем там else?

Код

C++ (Qt)
const T &operator[] (int i) const {
    if (i > 0 && i < size())
        return _arr[i];
    return T(); // а вот здесь компилятор выдаст предупреждение)) Какое? 
}
 

В Array1D<T> я так сделал, ради симетрии)) Симметрия - это круто, неправда ли))
А теперь попробуйте реализовать проверку в случае двумерного массива)

А вот в ф-ию рассчета гаусса булку передавать надо, и возвращать она должна не эррор, а НОВЫЙ вектор-столбец, который является решением ур-я.

Тож посещала меня такая мысля) Но эт накладно. Нужно в функции создавать объект Array1D<double> x(rows);
Затем при присваивании, опять копировать из одного вектора в другой.. Нее..
У меня оч специфичная задача и там большая потребность в скорости.

Спасибо за замечания)

[Авварон]

// а вот здесь компилятор выдаст предупреждение)) Какое?

возврат ссылки на локальный объект? Лень проверять) Это ваше дело, как заставить это работать:) Насчет квадратных скобочек vs. круглых для матриц 2на2 я бы использовал квадратные, а размерностью выше - круглые, но исключительно из соображений удобства.
Используйте шаред поинтер. Удобство использования часто должно превалировать над скоростью (иначе передавайте везде указатели, что уж там:)). И btw у вас и так создается объект)

[Igors]

Код

C++ (Qt)
 x(i) = (extMatrix(i, columns-1)-s)/extMatrix(i, i);
 

Не наблюдаю проверки на ноль

Код

C++ (Qt)
const T &operator() (int row, int column, bool *isValid = 0) const;
 

Воэиться с isVaiid при каждом вызове нереально, по-хорошему надо возбуждать exception,  Если лень, можно обойтись напр так

Код

C++ (Qt)
const T &operator() (int row, int column, bool *isValid = 0) const
{
#if DEBUG_DIAGNOSE
 if (row < 0 || row >= rows || column < 0 || column >= columns)
   DebugMessage("Out of range");
#endif
 ..
}
 

[m_ax]

Код
C++ (Qt)
 x(i) = (extMatrix(i, columns-1)-s)/extMatrix(i, i);
 
Не наблюдаю проверки на ноль

Давайте посмотрим по-внимательнее:

Код

C++ (Qt)
for (int i = 0; i < columns-1; ++i) {
        double maxVal = extMatrix(i, i);
        int index = i;
        for (int j = i+1; j < rows; ++j) {                   // Вот в этом цикле мы ищем максимальное значение
            if (fabs(extMatrix(j, i)) > maxVal) {        // элемента, на который будем, в последствии делить..
                maxVal = extMatrix(j, i);                 // и переставляем ту строку где его нашли на самый верх
                index = j;
            }
        }
        if (fabs(maxVal) < epsilon)                     // А вот здесь внимание!! Если этот максимальный элемент < epsilon
            return SYSTEM_INCOMPATIBLE;     // то выходим со статусом SYSTEM_INCOMPATIBLE это и есть проверка! 
 
        if (index != i) {
            for (int j = 0; j < columns; ++j) {
                double x = extMatrix(index, j);
                extMatrix(index, j) = extMatrix(i, j);
                extMatrix(i, j) = x;
            }
        }
        for (int r = i+1; r < rows; ++r) {
            for (int c = columns-1; c >= i; --c) {
                double s = -extMatrix(r, i)/extMatrix(i, i);       // Здесь деление на ноль уже исключено
                extMatrix(r, c) += (s*extMatrix(i, c));
            }
        }
    }
 

Я Вас убедил?

Код

Код

C++ (Qt)
const T &operator() (int row, int column, bool *isValid = 0) const;
 

Воэиться с isVaiid при каждом вызове нереально, по-хорошему надо возбуждать exception, 

А мы и не возимся с валид) Я вызываю элемент массива просто m(i, j) если уверен, что не выйду за пределы..
Кстати, даже если и выйду, то грохинга программы не произойдёт)) Я буду иметь дело с Tresh-ом.

Мне в одной ситуации (если интересно могу рассказать) очень был удобен такой финт с трэшом и указателем isValid..

Я всё же не вижу причин кидать исключения: с их отловом и обработкой возни ещё больше, на мой взгляд)

Благодарю за обсуждение)   

[Igors]

Давайте посмотрим по-внимательнее:

Код

C++ (Qt)
        double maxVal = extMatrix(i, i);
        int index = i;
        for (int j = i+1; j < rows; ++j) {                   // Вот в этом цикле мы ищем максимальное значение
            if (fabs(extMatrix(j, i)) > maxVal) {        // элемента, на который будем, в последствии делить..
                maxVal = extMatrix(j, i);                 // и переставляем ту строку где его нашли на самый верх
                index = j;
            }
        }
        if (fabs(maxVal) < epsilon)                     // А вот здесь внимание!! Если этот максимальный элемент < epsilon
            return SYSTEM_INCOMPATIBLE;     // то выходим со статусом SYSTEM_INCOMPATIBLE это и есть проверка! 
 

Здесь пробой со знаком maxVal. Лучше так

Код

C++ (Qt)
        int index = i;
        for (int j = i + 1; j < rows; ++j)               
            if (fabs(extMatrix(j, i)) > fabs(extMatrix(index, i))
                index = j;
 
        T maxVal = extMatrix(index, i);
        if (fabs(maxVal) < epsilon)      
            return SYSTEM_INCOMPATIBLE;  
 

[m_ax]

Igors
Да, действительно) Спасибо, что поправили, так бы сам проглядел))

[Admin]

На будущее используйте либу
http://www.gnu.org/software/gsl/

там практически есть все, что есть в современной математике.

[m_ax]

На будущее используйте либу
http://www.gnu.org/software/gsl/

там практически есть все, что есть в современной математике.

Да, спасибо за ссылку) Посмотрел, там даже 9-j символы есть))
Жаль она на чистом C.. Зато исходники открыты)