Алгоритм Томаса

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

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

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

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

${\displaystyle {\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}}.}$

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

Розв'язок

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

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

та

${\displaystyle 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}}}$

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

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

${\displaystyle 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]; } 

Доведення

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

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

${\displaystyle {\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}}}$

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

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

що дасть:

${\displaystyle (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}\,}$

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

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

${\displaystyle (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}\,}$

${\displaystyle (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}.\,}$

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

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

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

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

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

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

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

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

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

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

${\displaystyle a'_{i}=0\,}$

${\displaystyle b'_{i}=1\,}$

${\displaystyle c'_{1}={\frac {c_{1}}{b_{1}}}\,}$

${\displaystyle c'_{i}={\frac {c_{i}}{b_{i}-c'_{i-1}a_{i}}}\,}$

${\displaystyle d'_{1}={\frac {d_{1}}{b_{1}}}\,}$

${\displaystyle d'_{i}={\frac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}}.\,}$

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

${\displaystyle {\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}}\,}$

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

${\displaystyle x_{n}=d'_{n}/b'_{n}\,}$

${\displaystyle 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 .