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

Матеріал з Вікіпедії — вільної енциклопедії.
Перейти до: навігація, пошук

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

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

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}

та

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{align}
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{align}

Розглянемо таку зміну другого (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 можна нормувати діленням (якщо немає ризику ділення на число занадто близьке до нуля), тепер змінені коефіцієнти позначені рисочкою будуть:

a'_i = 0\,
b'_i = 1\,
c'_1 = \frac{c_1}{b_1}\,
c'_i = \frac{c_i}{b_i - c'_{i - 1} a_i}\,
d'_1 = \frac{d_1}{b_1}\,
d'_i = \frac{d_i - d'_{i - 1} a_i}{b_i - c'_{i - 1} a_i}.\,

це дає нам наступну систему з тими самими невідомими і коефіцієнтами вираженими через початкові:

\begin{array}{lcl}
b'_i x_i + c'_i x_{i + 1} = d'_i \qquad &;& \ i = 1, \ldots, n - 1 \\
b'_n x_n = d'_n \qquad &;& \ i = n. \\
\end{array}
\,

Останнє рівняння містить лише одне невідоме. Розв'язуючи його, приводимо наступне останнє рівняння до одного невідомого. І так далі:

x_n = d'_n/b'_n\,
x_i = (d'_i - c'_i x_{i + 1})/b'_i \qquad ; \ i = n - 1, n - 2, \ldots, 1.

Примітки[ред.ред. код]

  1. Pradip Niyogi (2006). Introduction to Computational Fluid Dynamics. Pearson Education India. с. 76. ISBN 978-81-7758-764-7.