Алгоритм трехдиагональной матрицы

В числовой линейной алгебре алгоритм трехдиагональной матрицы , также известный как алгоритм Томаса (названный в честь Ллевеллина Томаса ), представляет собой упрощенную форму исключения Гаусса , которую можно использовать для решения трехдиагональных систем уравнений . Трехдиагональную систему для n неизвестных можно записать как

$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}}.$

Для таких систем решение можно получить в $O(n)$ операции вместо $O(n^{3})$ требуется методом исключения Гаусса . Первая проверка устраняет $a_{i}$'s, а затем (сокращенная) обратная замена дает решение. Примеры таких матриц обычно возникают в результате дискретизации одномерного уравнения Пуассона и интерполяции естественным кубическим сплайном .

Алгоритм Томаса нестабилен в целом, но является таковым в некоторых особых случаях, например, когда матрица является диагонально доминирующей (либо по строкам, либо по столбцам) или симметричной положительно определенной ; более точную характеристику устойчивости алгоритма Томаса см. в теореме Хайэма 9.12. Если в общем случае требуется стабильность, метод исключения Гаусса с частичным поворотом (GEPP). вместо этого рекомендуется использовать

Метод

Прямой анализ состоит из вычисления новых коэффициентов следующим образом, обозначая новые коэффициенты штрихами:

$c'_{i}={\begin{cases}{\cfrac {c_{i}}{b_{i}}},&i=1,\\{\cfrac {c_{i}}{b_{i}-a_{i}c'_{i-1}}},&i=2,3,\dots ,n-1\end{cases}}$

и

$d'_{i}={\begin{cases}{\cfrac {d_{i}}{b_{i}}},&i=1,\\{\cfrac {d_{i}-a_{i}d'_{i-1}}{b_{i}-a_{i}c'_{i-1}}},&i=2,3,\dots ,n.\end{cases}}$

Тогда решение получается обратной заменой:

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

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

Приведенный выше метод не изменяет исходные векторы коэффициентов, но также должен отслеживать новые коэффициенты. Если векторы коэффициентов можно изменить, то алгоритм с меньшими затратами на бухгалтерский учет будет выглядеть следующим образом:

Для $i=2,3,\dots ,n,$ делать

$w={\cfrac {a_{i}}{b_{i-1}}},$

$b_{i}:=b_{i}-wc_{i-1},$

$d_{i}:=d_{i}-wd_{i-1},$

с последующей обратной заменой

$x_{n}={\cfrac {d_{n}}{b_{n}}},$

$x_{i}={\cfrac {d_{i}-c_{i}x_{i+1}}{b_{i}}}\quad {\text{for }}i=n-1,n-2,\dots ,1.$

Реализация в виде функции C , которая использует рабочее пространство, чтобы избежать изменения входных данных для переменного тока, что позволяет их повторно использовать:

void thomas(const int X, double x[restrict X],
            const double a[restrict X], const double b[restrict X],
            const double c[restrict X], double scratch[restrict X]) {
    /*
     solves Ax = d, where A is a tridiagonal matrix consisting of vectors a, b, c
     X = number of equations
     x[] = initially contains the input v, and returns x. indexed from [0, ..., X - 1]
     a[] = subdiagonal, indexed from [1, ..., X - 1]
     b[] = main diagonal, indexed from [0, ..., X - 1]
     c[] = superdiagonal, indexed from [0, ..., X - 2]
     scratch[] = scratch space of length X, provided by caller, allowing a, b, c to be const
     not performed in this example: manual expensive common subexpression elimination
     */
    scratch[0] = c[0] / b[0];
    x[0] = x[0] / b[0];

    /* loop from 1 to X - 1 inclusive */
    for (int ix = 1; ix < X; ix++) {
        if (ix < X-1){
        scratch[ix] = c[ix] / (b[ix] - a[ix] * scratch[ix - 1]);
        }
        x[ix] = (x[ix] - a[ix] * x[ix - 1]) / (b[ix] - a[ix] * scratch[ix - 1]);
    }

    /* loop from X - 2 to 0 inclusive */
    for (int ix = X - 2; ix >= 0; ix--)
        x[ix] -= scratch[ix] * x[ix + 1];
}

Вывод

Вывод алгоритма трехдиагональной матрицы является частным случаем исключения Гаусса .

Предположим, что неизвестные $x_{1},\ldots ,x_{n}$, и что уравнения, которые необходимо решить, следующие:

${\begin{alignedat}{4}&&&b_{1}x_{1}&&+c_{1}x_{2}&&=d_{1}\\&a_{i}x_{i-1}&&+b_{i}x_{i}&&+c_{i}x_{i+1}&&=d_{i}\,,\quad i=2,\ldots ,n-1\\&a_{n}x_{n-1}&&+b_{n}x_{n}&&&&=d_{n}\,.\end{alignedat}}$

Рассмотрите возможность изменения второго ( $i=2$) уравнение с первым уравнением следующим образом:

$({\mbox{equation 2}})\cdot b_{1}-({\mbox{equation 1}})\cdot 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}$ было исключено из второго уравнения. Использование аналогичной тактики с модифицированным вторым уравнением третьего уравнения дает:

$(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^{th}$ ряд; (измененный) $n^{th}$ уравнение будет включать только одно неизвестное, $x_{n}$. Это можно решить, а затем использовать для решения $(n-1)^{th}$ уравнение и так далее, пока все неизвестные не будут решены.

Очевидно, что коэффициенты модифицированных уравнений становятся все более и более сложными, если это указано явно. Изучив процедуру, измененные коэффициенты (обозначенные тильдами) вместо этого могут быть определены рекурсивно:

${\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}x_{i}+c'_{i}x_{i+1}=d'_{i}\qquad &;&\ i=1,\ldots ,n-1\\x_{n}=d'_{n}\qquad &;&\ i=n.\\\end{array}}\,$

Последнее уравнение содержит только одно неизвестное. Решение его, в свою очередь, сводит следующее последнее уравнение к одному неизвестному, так что эту обратную замену можно использовать для нахождения всех неизвестных:

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

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

Варианты

В некоторых ситуациях, особенно в тех, которые связаны с периодическими граничными условиями , может потребоваться решение слегка возмущенной формы трехдиагональной системы:

${\begin{alignedat}{4}&a_{1}x_{n}&&+b_{1}x_{1}&&+c_{1}x_{2}&&=d_{1}\\&a_{i}x_{i-1}&&+b_{i}x_{i}&&+c_{i}x_{i+1}&&=d_{i}\,,\quad i=2,\ldots ,n-1\\&a_{n}x_{n-1}&&+b_{n}x_{n}&&+c_{n}x_{1}&&=d_{n}\,.\end{alignedat}}$

В этом случае мы можем использовать формулу Шермана-Моррисона , чтобы избежать дополнительных операций исключения Гаусса и при этом использовать алгоритм Томаса. Метод требует решения модифицированной нециклической версии системы как для входных данных, так и для разреженного корректирующего вектора, а затем объединения решений. Это можно сделать эффективно, если оба решения вычисляются одновременно, поскольку прямая часть алгоритма чистой трехдиагональной матрицы может использоваться совместно.

Если мы укажем: $$A={\begin{bmatrix}b_{1}&c_{1}&&&a_{1}\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\c_{n}&&&a_{n}&b_{n}\end{bmatrix}},x={\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\\\vdots \\x_{n}\end{bmatrix}},d={\begin{bmatrix}d_{1}\\d_{2}\\d_{3}\\\vdots \\d_{n}\end{bmatrix}}$$

Тогда решаемая система: $$Ax=d$$

В этом случае коэффициенты $a_{1}$ и $c_{n}$ вообще говоря, ненулевые, поэтому их наличие не позволяет напрямую применить алгоритм Томаса. Поэтому мы можем рассмотреть $B\in \mathbb {R} ^{n\times n}$ и $u,v\in \mathbb {R} ^{n}$ следующим образом: $$B={\begin{bmatrix}b_{1}-\gamma &c_{1}&&&0\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\0&&&a_{n}&b_{n}-{\frac {c_{n}a_{1}}{\gamma }}\end{bmatrix}},u={\begin{bmatrix}\gamma \\0\\0\\\vdots \\c_{n}\end{bmatrix}},v={\begin{bmatrix}1\\0\\0\\\vdots \\a_{1}/\gamma \end{bmatrix}}.$$ Где $\gamma \in \mathbb {R}$ — параметр, который необходимо выбрать. Матрица $A$ можно реконструировать как $A=B+uv^{\mathsf {T}}$. Тогда решение получается следующим образом: сначала решаем две трехдиагональные системы уравнений, применяя алгоритм Томаса: $$By=d\qquad \qquad Bq=u$$

Затем восстанавливаем решение $x$ используя формулу Шермана-Моррисона : $${\begin{aligned}x&=A^{-1}d=(B+uv^{T})^{-1}d=B^{-1}d-{\frac {B^{-1}uv^{T}B^{-1}}{1+v^{T}B^{-1}u}}d=y-{\frac {qv^{T}y}{1+v^{T}q}}\end{aligned}}$$

Реализация в виде функции C , которая использует рабочее пространство, чтобы избежать изменения входных данных для переменного тока, что позволяет их повторно использовать:

void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double u[restrict X]) {
    /*
     solves Ax = v, where A is a cyclic tridiagonal matrix consisting of vectors a, b, c
     X = number of equations
     x[] = initially contains the input v, and returns x. indexed from [0, ..., X - 1]
     a[] = subdiagonal, regularly indexed from [1, ..., X - 1], a[0] is lower left corner
     b[] = main diagonal, indexed from [0, ..., X - 1]
     c[] = superdiagonal, regularly indexed from [0, ..., X - 2], c[X - 1] is upper right
     cmod[], u[] = scratch vectors each of length X
     */

    /* lower left and upper right corners of the cyclic tridiagonal system respectively */
    const double alpha = a[0];
    const double beta = c[X - 1];

    /* arbitrary, but chosen such that division by zero is avoided */
    const double gamma = -b[0];

    cmod[0] = alpha / (b[0] - gamma);
    u[0] = gamma / (b[0] - gamma);
    x[0] /= (b[0] - gamma);

    /* loop from 1 to X - 2 inclusive */
    for (int ix = 1; ix + 1 < X; ix++) {
        const double m = 1.0 / (b[ix] - a[ix] * cmod[ix - 1]);
        cmod[ix] = c[ix] * m;
        u[ix] = (0.0f  - a[ix] * u[ix - 1]) * m;
        x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
    }

    /* handle X - 1 */
    const double m = 1.0 / (b[X - 1] - alpha * beta / gamma - beta * cmod[X - 2]);
    u[X - 1] = (alpha    - a[X - 1] * u[X - 2]) * m;
    x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;

    /* loop from X - 2 to 0 inclusive */
    for (int ix = X - 2; ix >= 0; ix--) {
        u[ix] -= cmod[ix] * u[ix + 1];
        x[ix] -= cmod[ix] * x[ix + 1];
    }

    const double fact = (x[0] + x[X - 1] * beta / gamma) / (1.0 + u[0] + u[X - 1] * beta / gamma);

    /* loop from 0 to X - 1 inclusive */
    for (int ix = 0; ix < X; ix++)
        x[ix] -= fact * u[ix];
}

Существует и другой способ решения рассмотренной выше слегка возмущенной формы трехдиагональной системы. Рассмотрим две вспомогательные линейные системы размерности $(n-1)\times (n-1)$: $${\begin{aligned}\qquad \ \ \ \ \ b_{2}u_{2}+c_{2}u_{3}&=d_{2}\\a_{3}u_{2}+b_{3}u_{3}+c_{3}u_{4}&=d_{3}\\a_{i}u_{i-1}+b_{i}u_{i}+c_{i}u_{i+1}&=d_{i}\\\dots \\a_{n}u_{n-1}+b_{n}u_{n}\qquad &=d_{n}\,.\end{aligned}}\quad i=4,\ldots ,n-1\qquad \qquad {\begin{aligned}\qquad \ \ \ \ \ b_{2}v_{2}+c_{2}v_{3}&=-a_{2}\\a_{3}v_{2}+b_{3}v_{3}+c_{3}v_{4}&=0\\a_{i}u_{i-1}+b_{i}u_{i}+c_{i}u_{i+1}&=0\\\dots \\a_{n}v_{n-1}+b_{n}v_{n}\qquad &=-c_{n}\,.\end{aligned}}\quad i=4,\ldots ,n-1$$

Для удобства дополнительно определим $u_{1}=0$ и $v_{1}=1$. Теперь мы можем найти решения $\{u_{2},u_{3}\dots ,u_{n}\}$ и $\{v_{2},v_{3}\dots ,v_{n}\}$ применение алгоритма Томаса к двум вспомогательным трехдиагональным системам.

Решение $\{x_{1},x_{2}\dots ,x_{n}\}$ тогда можно представить в виде: $$x_{i}=u_{i}+x_{1}v_{i}\qquad i=1,2,\dots ,n$$

Действительно, умножив каждое уравнение второй вспомогательной системы на $x_{1}$, сложив с соответствующим уравнением первой вспомогательной системы и воспользовавшись представлением $x_{i}=u_{i}+x_{1}v_{i}$, мы сразу видим, что число уравнений $2$ через $n$ исходной системы удовлетворены; осталось только удовлетворить номер уравнения $1$. Для этого рассмотрим формулу для $i=2$ и $i=n$ и заменить $x_{2}=u_{2}+x_{1}v_{2}$и $x_{n}=u_{n}+x_{1}v_{n}$ в первое уравнение исходной системы. Это дает одно скалярное уравнение для $x_{1}$: $$b_{1}x_{1}+c_{1}(u_{2}+x_{1}v_{2})+a_{1}(u_{n}+x_{1}v_{n})=d_{1}$$

Таким образом, мы находим: $$x_{1}={\frac {d_{1}-a_{1}u_{n}-c_{1}u_{2}}{b_{1}+a_{1}v_{n}+c_{1}v_{2}}}$$

Реализация в виде функции C , которая использует рабочее пространство, чтобы избежать изменения входных данных для переменного тока, что позволяет их повторно использовать:

void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double v[restrict X]) {
    /* first solve a system of length X - 1 for two right hand sides, ignoring ix == 0 */
    cmod[1] = c[1] / b[1];
    v[1] = -a[1] / b[1];
    x[1] = x[1] / b[1];

    /* loop from 2 to X - 1 inclusive */
    for (int ix = 2; ix < X - 1; ix++) {
        const double m = 1.0 / (b[ix] - a[ix] * cmod[ix - 1]);
        cmod[ix] = c[ix] * m;
        v[ix] = (0.0f  - a[ix] * v[ix - 1]) * m;
        x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
    }

    /* handle X - 1 */
    const double m = 1.0 / (b[X - 1] - a[X - 1] * cmod[X - 2]);
    cmod[X - 1] = c[X - 1] * m;
    v[X - 1] = (-c[0]    - a[X - 1] * v[X - 2]) * m;
    x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;

    /* loop from X - 2 to 1 inclusive */
    for (int ix = X - 2; ix >= 1; ix--) {
        v[ix] -= cmod[ix] * v[ix + 1];
        x[ix] -= cmod[ix] * x[ix + 1];
    }

    x[0] = (x[0] - a[0] * x[X - 1] - c[0] * x[1]) / (b[0] + a[0] * v[X - 1] + c[0] * v[1]);

    /* loop from 1 to X - 1 inclusive */
    for (int ix = 1; ix < X; ix++)
        x[ix] += x[0] * v[ix];
}

В обоих случаях решаемые вспомогательные системы действительно трехдиагональны, поэтому общая вычислительная сложность решения системы $Ax=d$ остается линейным относительно размерности системы $n$, то есть $O(n)$ арифметические операции.

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

учебнике «Числовая математика» В Альфио Квартерони , Сакко и Салери приведена модифицированная версия алгоритма, которая позволяет избежать некоторых делений (вместо этого используется умножение), что полезно для некоторых компьютерных архитектур.

Параллельные трехдиагональные решатели были опубликованы для многих векторных и параллельных архитектур, включая графические процессоры.

Подробное описание параллельных трехдиагональных и блочных трехдиагональных решателей см.

Ссылки

• Конте, SD; де Бур, К. (1972). Элементарный численный анализ . МакГроу-Хилл, Нью-Йорк. ISBN  0070124469 .
• Эта статья включает текст из статьи Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) на CFD-Wiki , которая находится под лицензией GFDL .

• Пресс, WH; Теукольский, С.А.; Феттерлинг, WT; Фланнери, BP (2007). «Раздел 2.4» . Численные рецепты: искусство научных вычислений (3-е изд.). Нью-Йорк: Издательство Кембриджского университета. ISBN  978-0-521-88068-8 .