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

Метод прогонки или алгоритм Томаса (англ. Thomas algorithm) используется для решения систем линейных уравнений вида $~Ax=F$, где A — трёхдиагональная матрица.

Описание метода

Система уравнений $~Ax=F$ равносильна соотношению

$~A_{i}x_{i-1}+C_{i}x_{i}+B_{i}x_{i+1} = F_{i}.\qquad\qquad(1)$

Метод прогонки основывается на предположении, что искомые неизвестные связаны рекуррентным соотношением:

$x_i = \alpha_{i+1}x_{i+1} + \beta_{i+1},\,\!$ где $~i=n-1,n-2,\dots,1.\qquad\qquad(2)$

Используя это соотношение, выразим x_i-1 и x_i через x_i+1 и подставим в уравнение (1):

$\left(A_i\alpha_i\alpha_{i+1} + C_i\alpha_{i+1} + B_i\right)x_{i+1} + A_i\alpha_i\beta_{i+1} + A_i\beta_i + C_i\beta_{i+1} - F_i = 0$,

где F_i — правая часть i-го уравнения. Это соотношение будет выполняться независимо от решения, если потребовать

$\begin{cases} A_i\alpha_i\alpha_{i+1} + C_i\alpha_{i+1} + B_i = 0\\ A_i\alpha_i\beta_{i+1} + A_i\beta_i + C_i\beta_{i+1} - F_i = 0 \end{cases}$

Отсюда следует:

$\begin{cases} \alpha_{i+1} = \frac{-B_i}{A_i\alpha_i + C_i} \\ \beta_{i+1} = \frac{F_i - A_i\beta_i}{A_i\alpha_i + C_i}\end{cases}$

Из первого уравнения получим:

$\begin{cases} \alpha_2 = \frac{-B_1}{C_1} \\ \beta_2 = \frac{F_1}{C_1}\end{cases}$

После нахождения прогоночных коэффициентов $\alpha$ и $\beta$, используя уравнение (2), получим решение системы. При этом,

$x_i = {\alpha_{i+1}x_{i+1} + \beta_{i+1}},\!$     $i=n-1..1 \,\!$

$x_n = \frac{F_n-A_n\beta_n}{C_n+A_n\alpha_n}$

Другим способом объяснения существа метода прогонки, более близким к терминологии конечно-разностных методов и объясняющим происхождение его названия, является следующий: преобразуем уравнение (1) к эквивалентному ему уравнению

$~A' x=F'\qquad\qquad(1')$

c надиагональной (наддиагональной) матрицей

$A' = \begin{pmatrix} C_1' & B_1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & C_2' & B_2 & 0 & \cdots & 0 & 0 \\ 0 & 0 & C_3' & B_3 & \cdots & 0 & 0 \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \cdots & \cdots & \cdots & \cdots & 0 & C'_{n-1} & B_{n-1} \\ 0 & 0 & 0 & 0 & \cdots & 0 & C_{n}' \end{pmatrix}$.

Вычисления проводятся в два этапа. На первом этапе вычисляются компоненты матрицы $~C_i'$ и вектора $~F'$, начиная с  $~i=2$  до  $~i=n:$

$~C'_1=C_1;$

$C'_i=C_i-\frac{A_{i-1} B_{i-1}}{C'_{i-1}},$

и

$~F'_1=F_1;$

$F'_i=F_i-\frac{A_{i-1} F'_{i-1}}{C'_{i-1}}.$

На втором этапе, для $i=n, n-1, \dots, 1$ вычисляется решение:

$x_n=\frac{F'_n}{C'_n};$

$x_i=\frac{F'_i-B_i x_{i+1}}{C'_i}.$

Такая схема вычисления объясняет также английский термин этого метода «shuttle».

Для применимости формул метода прогонки достаточно свойства строгого диагонального преобладания у матрицы A.

Пример реализации на Си

Данный код работает при предположении, что a[0] = 0, b[n-1] = 0.

        /**
         * n - число уравнений (строк матрицы)
         * b - диагональ, лежащая над главной (нумеруется: [0;n-2])
         * c - главная диагональ матрицы A (нумеруется: [0;n-1])
         * a - диагональ, лежащая под главной (нумеруется: [1;n-1])
         * f - правая часть (столбец)
         * x - решение, массив x будет содержать ответ
         */
void solveMatrix (int n, double *a, double *c, double *b, double *f, double *x)
{
        double m;
        for (int i = 1; i < n; i++)
        {
                m = a[i]/c[i-1];  
                c[i] = c[i] - m*b[i-1];
                f[i] = f[i] - m*f[i-1];
        }
 
        x[n-1] = f[n-1]/c[n-1]; 
 
        for (int i = n - 2; i >= 0; i--) 
                x[i]=(f[i]-b[i]*x[i+1])/c[i];
 
}

Пример реализации на Java

import java.util.Scanner;
 
public class Prorace_method {
    public static void main(String[] args) {
        int n;
        Scanner sc = new Scanner (System.in);
        System.out.println("Введите число уравнений:");
        n = sc.nextInt();
        double a[][] = new double[n+1][n+1];
        double b[] = new double[n+1];
        double x[] = new double[n+1];
        double k[] = new double[n+1];
        double m[] = new double[n+1];
        double t[] = new double[n+1];
        double p[] = new double[n+1];
        double q[] = new double[n+1];               
        for (int i = 1; i <= n; i++) {
            for (int j = 1; j <= n; j++) {
                System.out.print("a(" + (i-1) + "," + (j-1) + ")=");
                a[i][j] = sc.nextDouble();
            }
            System.out.print("b(" + (i-1) + ")=");
            b[i] = sc.nextDouble();
        }
        for (int i = 1; i <= n; i++) {
            if (i == 1) { 
                k[i] = 0;
            }
            else {
                k[i] = a[i][i - 1];
            }
            m[i] = -a[i][i];            
            if (i == n ) {
                t[i] = 0;
            } else {
                t[i] = a[i][i + 1];
            }
        }
        p[1] = t[1] / m[1];
        q[1] = -b[1] / m[1];
        for (int i = 2; i <= n; i++) {
            p[i] = -t[i] / (k[i] * p[i - 1] - m[i]);
            q[i] = (b[i] - k[i] * q[i - 1]) / (k[i] * p[i - 1] - m[i]);
        }
        x[n] = (b[n] - k[n] * q[n - 1]) / (k[n] * p[n - 1] - m[n]);
        for (int i = n - 1; i >= 1; i--) {
            x[i] = p[i] * x[i + 1] + q[i];
        }
        for (int i = 1; i <=n; i++) {
            System.out.println("x[" + i + "] = " + x[i]);
        }
    }
}

Литература

• Introduction to Numerical Methods in Chemical Engineering 1.1 "Tridiagonal matrix algorithm (TDMA)"  (англ.)
• А. А. Абрамов, В. Б. Андреев. О применении метода прогонки к нахождению периодических решений дифференциальных и разностных уравнений. — Ж. вычисл. матем. и матем. физ., 3:2 (1963), 377–381

Ссылки

• The Thomas Algorithm  (англ.)

В этой статье не хватает ссылок на источники информации. Информация должна быть проверяема, иначе она может быть поставлена под сомнение и удалена. Вы можете отредактировать эту статью, добавив ссылки на авторитетные источники. Эта отметка установлена 9 октября 2011.