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

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


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











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






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







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




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




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




Ниже представлен код, реализующий метод прогонки, принимает трехдиагональную матрицу 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;
}