суббота, 17 ноября 2007 г.

Метод прогонки

Метод прогонки для решения систем линейных алгебраических уравнений(СЛАУ).


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











Поделив первую строку матрицы, приведенной выше, на -b1 очевидно, что:






и можно вывести формулу для прямого хода:







Затем необходимо выполнить обратный ход - найти вектор X, из последней строки преобразованной матрицы следует, что xn= Qn.
В тоже время остальные элементы вектора считаются по формуле:




Следует заметить, что метод устойчив если(следует из диагонального преобладания матрицы А):




и корректен, если(иначе формулы прямого хода не имеют смысла):




Ниже представлен код, реализующий метод прогонки, принимает трехдиагональную матрицу a размерности N*N, и вектор правых частей b размерности N, результат возвращается в b:


void sweep(double a[N][N],double b[N])
{
int i;
double znam;

b[0]/=a[0][0];//Q1
a[0][1]/=-a[0][0];//P1

for(i=1;i < N-1;i++)
{
znam=-a[i][i]-a[i][i-1]*a[i-1][i]; //общий знаменатель для формул нахождения Pi, Qi
a[i][i+1]/=znam; //Pi
b[i]=(a[i][i-1]*b[i-1]-b[i])/znam; //Qi

}
//строка ниже для вычисления QN
b[N-1]=(a[N-1][N-2]*b[N-2]-b[N-1])/(-a[N-1][N-1]-a[N-1][N-2]*a[N-2][N-1]);


//обратный ход
for(i=N-2;i > -1;i--)
{
b[i]+=b[i+1]*a[i][i+1];
}

return;
}

22 комментария:

Анонимный комментирует...

Статья гавно.
Во-первых нигде не обьясняется по-человечески как получается рекуррентные формулы для P и Q,а только говорится: "и можно вывести формулу для прямого хода".Что значит можно? Можно бабе залететь,когда ее трахают без резинки.А понять как выводятся коэффициенты без подробного обьяснения "не можно".

Во-вторых код программы не пашет нормально.Есть ошибочки в значениях неизвестных.

Алексей комментирует...

Странно, а у меня все заработало. Спасибо автору статьи.

Жнец комментирует...

Хреновая статья. Многие места просто не понятны. Код программы, а точнее прописанный в ней алгоритм - в целом рабочий, но следовало написать, на чём его запускали. На borland c++ 3.01, к примеру, исконный код выдаст кучу ошибок.

Анонимный комментирует...

Хорошая статья, сразу понятно где a b c диагонали, что такое P q , да и код лучше чем 3 листа задротной инмтрукции

kaf комментирует...

Спасиьо автору за статьи, удивлен, что кто-то ещё недовольство высказывает :)

Анонимный комментирует...

Нормальная статья. Все по-простому. Спасибо автору.

Анонимный комментирует...

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

Анонимный комментирует...

если кому-то интересно, откуда выводятся коэффициенты, можно посмотреть тут http://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%BF%D1%80%D0%BE%D0%B3%D0%BE%D0%BD%D0%BA%D0%B8

Анонимный комментирует...

Мне так кажется, что автор допустил пару ошибок в формулах.
для Q1=d1/b1; Qi=-(aiQi-1-di)/(bi-aiPi-1).

Citizen Minister комментирует...

Обсирают статью мудаки и гандоны!!!!!!
Ископал весь рунет, единственный рабочий код. Автор, СПАСИБО!!!

На Visual Studio 2008 все зашибись работает, тестовый пример на отлично решает...

Замира комментирует...

Реализация не плохая, но не экономная. Можно хранить 3 диагонали в массиве размером 3 на N, а потом оперировать этими элементами. А тут жрут память N на N.

1+1 комментирует...

К сожалению, код однозначно работает неверно. Даже на простых примерах.

Анонимный комментирует...

У меня код однозначно работает правильно. И это подтвердил даже препод.

Анонимный комментирует...

Код нормальный чуть подпаравить и все отлично работает))

Анонимный комментирует...

Глупо передавать в метод огромную матрицу, если можно обойтись 3мя массивами.

Анонимный комментирует...

ну если кому-то непонятен вывод формул в статье, то Самарский "Введение в теорию разностных схем" или Самарский, Николаев "Методы решения сеточных уравнений"...
Лучше скажите, где ошибка в коде (если она есть), чтобы автор исправился))

Андрей комментирует...

Все роботает!
Спосибо!

Анонимный комментирует...

Учите МАТЧАСТЬ епт, потом уже разбирайтесь

Анонимный комментирует...

Хороша стаття.Код працює.

Анонимный комментирует...

Формулы получаются так:
Первуя строку делим на ее первый элемент чтобы первый элемент стал равным 1.
Из второй строки вычитаем первую умноженную на соответвтвующий множитель так чтобы первый элемент второй строки стал равным 0. После этого делим вторую строку на ее второй элемент так чтобы он стал равным 1.
Остальные шаги выполняются аналогично.

Анонимный комментирует...

нихера не работает, плохое объяснение

Анонимный комментирует...

Не знаю как код, но объяснение весьма доходчивое, - все просто и ясно!