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

Метод прогонки, також відомий як алгоритм Томаса, дозволяє розв'язувати СЛАР з Тридіагональною матрицею, і є спрощенням методу Гауса для таких обмежень. Працює за лінійний час.

Система має такий вигляд:

$a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},\,\!$

де $a_{1}=0\,$ та $c_{n}=0\,$ . В матричній формі це записується так:

${\begin{bmatrix}{b_{1}}&{c_{1}}&{}&{}&{0}\\{a_{2}}&{b_{2}}&{c_{2}}&{}&{}\\{}&{a_{3}}&{b_{3}}&\ddots &{}\\{}&{}&\ddots &\ddots &{c_{n-1}}\\{0}&{}&{}&{a_{n}}&{b_{n}}\\\end{bmatrix}}{\begin{bmatrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\\vdots \\{x_{n}}\\\end{bmatrix}}={\begin{bmatrix}{d_{1}}\\{d_{2}}\\{d_{3}}\\\vdots \\{d_{n}}\\\end{bmatrix}}.$

В цілому, метод не є числово стійким, але є таким у декількох випадках, таких як діагонально панівна матриця або додатноозначена матриця.^[1]

Розв'язок

Розв'язок проводиться в два кроки, як і в методі Гауса, прямому, та зворотному. В прямому ході ми обчислюємо:

$c'_{1}={\frac {c_{1}}{b_{1}}};\ c'_{i}={\frac {c_{i}}{b_{i}-c'_{i-1}a_{i}}},\,i={\overline {2,n-1}}$

та

$d'_{1}={\frac {d_{1}}{b_{1}}};\ d'_{i}={\frac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}},\,i={\overline {2,n}}$

Тепер розв'язок знаходимо зворотнім ходом:

$x_{n}=d'_{n}\,$

$x_{i}=d'_{i}-c'_{i}x_{i+1};\ i=n-1,n-2,\ldots ,1$

Код на C

/* Розв'язок повертається в x. c та d модифікуються!*/

void TridiagonalSolve (const double *a, const double *b, double *c, double *d, double *x, unsigned int n){
 
	/* Modify the coefficients. */
	c[0] /= b[0];	/* Division by zero risk. */
	d[0] /= b[0];	/* Division by zero would imply a singular matrix. */
	for (unsigned int i = 1; i < n; i++){
		double id = 1 / (b[i] - c[i-1] * a[i]);  /* Division by zero risk. */
		c[i] *= id;	                         /* Last value calculated is redundant. */
		d[i] = (d[i] - d[i-1] * a[i]) * id;
	}
 
	/* Now back substitute. */
	x[n - 1] = d[n - 1];
	for (int i = n - 2; i >= 0; i--)
		x[i] = d[i] - c[i] * x[i + 1];
}

Доведення

Доведення методу вимагає ручного виконання деяких спеціалізованих Гаусових вилучань.

Припустимо, що невідомі - це $x_{1},\ldots ,x_{n}$ , і що рівняння до розв'язання такі:

{\begin{aligned}b_{1}x_{1}+c_{1}x_{2}&=d_{1};&i&=1\\a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}&=d_{i};&i&=2,\ldots ,n-1\\a_{n}x_{n-1}+b_{n}x_{n}&=d_{n};&i&=n.\end{aligned}}

Розглянемо таку зміну другого ( $i=2$ ) рівняння за допомогою першого рівняння:

({\mbox{equation 2}})\cdot b_{1}-({\mbox{equation 1}})\cdot a_{2}

що дасть:

(a_{2}x_{1}+b_{2}x_{2}+c_{2}x_{3})b_{1}-(b_{1}x_{1}+c_{1}x_{2})a_{2}=d_{2}b_{1}-d_{1}a_{2}\,

(b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3}=d_{2}b_{1}-d_{1}a_{2}\,

у висліді маємо, що $x_{1}$ було вилучено з другого рівняння. Використовуючи подібну тактику зі зміненим другим рівнянням, щодо третього маємо:

(a_{3}x_{2}+b_{3}x_{3}+c_{3}x_{4})(b_{2}b_{1}-c_{1}a_{2})-((b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3})a_{3}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}\,

(b_{3}(b_{2}b_{1}-c_{1}a_{2})-c_{2}b_{1}a_{3})x_{3}+c_{3}(b_{2}b_{1}-c_{1}a_{2})x_{4}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}.\,

Цього разу було вилучено $x_{2},$ якщо повторювати цю процедуру до рядка $n$ ; (змінене) рівняння $n$ міститиме лише одну змінну, $x_{n}$ . Таке рівняння ми можемо розв'язати і використати результат для того, щоб розв'язати рівняння $(n-1),$ і так далі допоки всі невідомі не знайдені.

Очевидно, що коефіцієнти у змінених рівняннях ставатимуть все більш заплутаними якщо розписувати їх явно. Але змінені коефіцієнти (з тильдою) можна виразити рекурсивно:

{\tilde {a}}_{i}=0\,

{\tilde {b}}_{1}=b_{1}\,

{\tilde {b}}_{i}=b_{i}{\tilde {b}}_{i-1}-{\tilde {c}}_{i-1}a_{i}\,

{\tilde {c}}_{1}=c_{1}\,

{\tilde {c}}_{i}=c_{i}{\tilde {b}}_{i-1}\,

{\tilde {d}}_{1}=d_{1}\,

{\tilde {d}}_{i}=d_{i}{\tilde {b}}_{i-1}-{\tilde {d}}_{i-1}a_{i}.\,

Для дальшого пришвидшення процесу, ${\tilde {b}}_{i}$ можна нормувати діленням (якщо немає ризику ділення на число занадто близьке до нуля), тепер змінені коефіцієнти позначені рисочкою будуть: