Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Алгоритмы > Решение СЛАУ в виде сжатой матрицы

Автор: alexgorbach 13.5.2013, 22:01

Добрый день. Помогите написать алгоритм для задачи.
Есть матрица жесткости, большая, симметричная, ленточная, вот такого вида:
http://firepic.org/
Красная и синяя область симметричны между собой, серые области - нулевые.
Так как большинство элементов матрицы нулевые, храню матрицу в сжатом виде, сдвинув каждую строку до упора влево, причем в силу симметрии храню только верхнюю половину значащих элементов (красную область). Примерно так:
http://firepic.org/
Так вот. Решаю систему методом Гаусса. Приписываю к этой матрице столбец свободных членов и передаю эту матрицу на обработку в мой метод, который должен бы решить систему и на место столбца свободных членов записать вектор в ответом, но! Работает неправильно. Я понимаю, почему, но не могу додуматься, как исправить.
Код метода на C#:

Код

public static bool Solve(double[,] M)
        {
            #region{Прямой ход}
            int rowCount = h;

            for (int i = 0; i + 1 < rowCount; i++)
            {
                for (int j = i + 1; j < Math.Min(i + tapeWidth, rowCount); j++)
                {
                    double x = M[i, j - i] / M[i, 0];
                    for (int p = 0; p < w - j + i - 1; p++)
                    {
                        M[j, p] -= M[i, j - i + p] * x;
                    }
                    M[j, w - 1] -= M[i, w - 1] * x;
                }
            }
            #endregion

            #region{Обратный ход}
            int k = 0; //индекс последнего значащего элемента в строке
            for (int row = rowCount - 1; row >= 0; row--)
            {
                for (int i = 1; i <= k; i++)
                {
                    M[row, w - 1] -= M[row, i] * M[row + i, w - 1];
                }
                M[row, w - 1] /= M[row, 0]; 
                if(k < w - 1)
                    k++;
            }
            return true;
            #endregion
        }

tapeWidth - ширина ленты, то бишь, ширина матрицы (не расширенной),
w, h - ширина и высота расширенной матрицы.

Проблема, как я понимаю, в прямом ходе:

Код

double x = M[i, j - i] / M[i, 0];

Здесь я вычисляю коэффициент, на который нужно умножить строку, чтобы при вычитании получить ноль в одной из последующих строк. Так как я храню только половину матрицы, то в числителе дроби выбираю не элемент, стоящий под главной диагональю, а симметричный ему над главной диагональю (я храню только верхнюю половину ленты).
Проблема в том, что после первого прохода (когда ведущая строка имеет номер 0), разные строки меняются по-разному (от них отнимается первая строка, умноженная на разные коэффициенты), и при последующих проходах элемент, симметричный нужному для вычисления коэффициента x, уже не будет равен тому, который подразумевается (который как бы под диагональю), и матрица преобразуется неправильно. Сложно объяснил, но надеюсь, кто-нибудь поймет.
Спасибо откликнувшимся.

Автор: Pavia 13.5.2013, 22:34

Надо смотреть в сторону ленточных алгоритмов. Я только в полных матрицах специализируюсь.
Думаю Гаус не подойдёт нужно применять ортогональные преобразования как в методе Якоби.
Гаусс симметрию нарушает в отличие от ортогонального преобразования в методе Якоби.

Автор: alexgorbach 13.5.2013, 22:46

Pavia, хм, тогда, быть может, стоит избрать другой способ хранения данных матрицы? Ибо мне советовали применять именно метод Гаусса, потому как система хорошо обусловлена и на главной диагонали стоят максимальные элементы, причем ненулевые. Гаусс отлично подходит для матрицы такого вида, только вот подходит ли он для выбранного способа хранения...

Автор: Фантом 13.5.2013, 23:11

А чем не годится прогонка?

Автор: alexgorbach 13.5.2013, 23:18

Фантом, прогонка чего? Не совсем понял...

Автор: Pavia 14.5.2013, 08:02

Есть такой алгоритм прогонки.
http://ru.wikipedia.org/wiki/Метод_прогонки
По сути тот же метод Гауса применённый к трёхдиагональной матрице.
К ленточной матрицы можно применить преобразование отражения(разновидность ортогонального) и свести её к трёх диагональной.

A*x=y
Разумеется при сведение матрицы к трёх диагональной операция умножения применяют над вектором "y".

Автор: Mirkes 14.5.2013, 10:31

Действительно, у вас на выбор два подхода:
1. Развернуть матрицу и использовать любые методы (Гаусса и тп)
2. Использовать методы сохраняющие симметрию (ортогональные разложения, прогонка и т.п.)
Решите что для Вас лучше и дейстщуйте.