Алгоритм трехдиагональной матрицы
Улучшенное сокращение для определенных матриц

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

а я х я − 1 + б я х я + с я х я + 1 = г я , {\displaystyle a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},}

где и . а 1 = 0 {\displaystyle a_{1}=0} с н = 0 {\displaystyle c_{n}=0}

[ б 1 с 1 0 а 2 б 2 с 2 а 3 б 3 ⋱ ⋱ ⋱ с н − 1 0 а н б н ] [ х 1 х 2 х 3 ⋮ х н ] = [ г 1 г 2 г 3 ⋮ г н ] . {\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}}.}

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

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

Метод

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

с я ′ = { с я б я , я = 1 , с я б я − а я с я − 1 ′ , я = 2 , 3 , … , н − 1 {\displaystyle 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}}}

и

г я ′ = { г я б я , я = 1 , г я − а я г я − 1 ′ б я − а я с я − 1 ′ , я = 2 , 3 , … , н . {\displaystyle 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 ′ , {\displaystyle x_{n}=d'_{n},}
x i = d i ′ − c i ′ x i + 1 , i = n − 1 , n − 2 , … , 1. {\displaystyle x_{i}=d'_{i}-c'_{i}x_{i+1},\quad i=n-1,n-2,\ldots ,1.}

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

Для делать i = 2 , 3 , … , n , {\displaystyle i=2,3,\dots ,n,}

w = a i b i − 1 , {\displaystyle w={\cfrac {a_{i}}{b_{i-1}}},}
b i := b i − w c i − 1 , {\displaystyle b_{i}:=b_{i}-wc_{i-1},}
d i := d i − w d i − 1 , {\displaystyle d_{i}:=d_{i}-wd_{i-1},}

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

x n = d n b n , {\displaystyle x_{n}={\cfrac {d_{n}}{b_{n}}},}
x i = d i − c i x i + 1 b i for  i = n − 1 , n − 2 , … , 1. {\displaystyle x_{i}={\cfrac {d_{i}-c_{i}x_{i+1}}{b_{i}}}\quad {\text{for }}i=n-1,n-2,\dots ,1.}

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

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 ]) { /*  решает Ax = d, где A — трехдиагональная матрица, состоящая из векторов a, b, c  X = количество уравнений  x[] = изначально содержит входные данные v и возвращает x. индексируется из [0, ..., X - 1]  a[] = поддиагональ, индексируется из [1, ..., X - 1]  b[] = главная диагональ, индексируется из [0, ..., X - 1]  c[] = наддиагональ, индексируется из [0, ..., X - 2]  scratch[] = пространство для царапин длиной X, предоставленное вызывающей стороной, позволяющее сделать a, b, c константами  в этом примере не выполняется: ручное дорогостоящее исключение общей подвыражения  */ scratch [ 0 ] = c [ 0 ] / b [ 0 ]; х [ 0 ] = х [ 0 ] / b [ 0 ];                                  /* цикл от 1 до X - 1 включительно */ 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 ]); }                                             /* цикл от X - 2 до 0 включительно */ for ( int ix = X - 2 ; ix >= 0 ; ix -- ) x [ ix ] -= scratch [ ix ] * x [ ix + 1 ]; }

Вывод

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

Предположим, что неизвестные равны , а уравнения, которые необходимо решить, следующие: x 1 , … , x n {\displaystyle x_{1},\ldots ,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 , i = 2 , … , n − 1 a n x n − 1 + b n x n = d n . {\displaystyle {\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 {\displaystyle i=2}

( equation 2 ) ⋅ b 1 − ( equation 1 ) ⋅ a 2 {\displaystyle ({\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 . {\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}.}

Обратите внимание, что было исключено из второго уравнения. Использование аналогичной тактики с измененным вторым уравнением в третьем уравнении дает: x 1 {\displaystyle 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 . {\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}.\,}

Это время было устранено. Если эту процедуру повторять до строки; (измененное) уравнение будет включать только одно неизвестное, . Это может быть решено для и затем использовано для решения уравнения, и так далее, пока все неизвестные не будут решены для. x 2 {\displaystyle x_{2}} n t h {\displaystyle n^{th}} n t h {\displaystyle n^{th}} x n {\displaystyle x_{n}} ( n − 1 ) t h {\displaystyle (n-1)^{th}}

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

a ~ i = 0 {\displaystyle {\tilde {a}}_{i}=0\,}
b ~ 1 = b 1 {\displaystyle {\tilde {b}}_{1}=b_{1}\,}
b ~ i = b i b ~ i − 1 − c ~ i − 1 a i {\displaystyle {\tilde {b}}_{i}=b_{i}{\tilde {b}}_{i-1}-{\tilde {c}}_{i-1}a_{i}\,}
c ~ 1 = c 1 {\displaystyle {\tilde {c}}_{1}=c_{1}\,}
c ~ i = c i b ~ i − 1 {\displaystyle {\tilde {c}}_{i}=c_{i}{\tilde {b}}_{i-1}\,}
d ~ 1 = d 1 {\displaystyle {\tilde {d}}_{1}=d_{1}\,}
d ~ i = d i b ~ i − 1 − d ~ i − 1 a i . {\displaystyle {\tilde {d}}_{i}=d_{i}{\tilde {b}}_{i-1}-{\tilde {d}}_{i-1}a_{i}.\,}

Чтобы еще больше ускорить процесс решения, можно разделить (если нет деления на нулевой риск), новые модифицированные коэффициенты, каждый из которых обозначен штрихом, будут иметь вид: b ~ i {\displaystyle {\tilde {b}}_{i}}

a i ′ = 0 {\displaystyle a'_{i}=0\,}
b i ′ = 1 {\displaystyle b'_{i}=1\,}
c 1 ′ = c 1 b 1 {\displaystyle c'_{1}={\frac {c_{1}}{b_{1}}}\,}
c i ′ = c i b i − c i − 1 ′ a i {\displaystyle c'_{i}={\frac {c_{i}}{b_{i}-c'_{i-1}a_{i}}}\,}
d 1 ′ = d 1 b 1 {\displaystyle d'_{1}={\frac {d_{1}}{b_{1}}}\,}
d i ′ = d i − d i − 1 ′ a i b i − c i − 1 ′ a i . {\displaystyle d'_{i}={\frac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}}.\,}

Это дает следующую систему с теми же неизвестными и коэффициентами, определенными через исходные, указанные выше:

x i + c i ′ x i + 1 = d i ′ ;   i = 1 , … , n − 1 x n = d n ′ ;   i = n . {\displaystyle {\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 ′ {\displaystyle x_{n}=d'_{n}\,}
x i = d i ′ − c i ′ x i + 1 ;   i = n − 1 , n − 2 , … , 1. {\displaystyle x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1.}

Варианты

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

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 , i = 2 , … , n − 1 a n x n − 1 + b n x n + c n x 1 = d n . {\displaystyle {\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 = [ b 1 c 1 a 1 a 2 b 2 c 2 a 3 b 3 ⋱ ⋱ ⋱ c n − 1 c n a n b n ] , x = [ x 1 x 2 x 3 ⋮ x n ] , d = [ d 1 d 2 d 3 ⋮ d n ] {\displaystyle 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}}}

Тогда система, которую нужно решить, имеет вид: A x = d {\displaystyle Ax=d}

В этом случае коэффициенты и , вообще говоря, не равны нулю, поэтому их наличие не позволяет напрямую применить алгоритм Томаса. Поэтому мы можем рассматривать и следующим образом: Где — параметр, который необходимо выбрать. Матрица может быть восстановлена ​​как . Тогда решение получается следующим образом: [4] сначала решаем две трехдиагональные системы уравнений, применяя алгоритм Томаса: a 1 {\displaystyle a_{1}} c n {\displaystyle c_{n}} B ∈ R n × n {\displaystyle B\in \mathbb {R} ^{n\times n}} u , v ∈ R n {\displaystyle u,v\in \mathbb {R} ^{n}} B = [ b 1 − γ c 1 0 a 2 b 2 c 2 a 3 b 3 ⋱ ⋱ ⋱ c n − 1 0 a n b n − c n a 1 γ ] , u = [ γ 0 0 ⋮ c n ] , v = [ 1 0 0 ⋮ a 1 / γ ] . {\displaystyle 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}}.} γ ∈ R {\displaystyle \gamma \in \mathbb {R} } A {\displaystyle A} A = B + u v T {\displaystyle A=B+uv^{\mathsf {T}}} B y = d B q = u {\displaystyle By=d\qquad \qquad Bq=u}

Затем реконструируем решение, используя формулу Шермана-Моррисона : x {\displaystyle x} x = A − 1 d = ( B + u v T ) − 1 d = B − 1 d − B − 1 u v T B − 1 1 + v T B − 1 u d = y − q v T y 1 + v T q {\displaystyle {\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 , которая использует рабочее пространство, чтобы избежать изменения своих входных данных для ac, что позволяет использовать их повторно:

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 ]) { /*  решает Ax = v, где A - циклическая трехдиагональная матрица, состоящая из векторов a, b, c  X = количество уравнений  x[] = изначально содержит входные данные v и возвращает x. индексируется из [0, ..., X - 1]  a[] = поддиагональ, регулярно индексируется из [1, ..., X - 1], a[0] - нижний левый угол  b[] = главная диагональ, индексируется из [0, ..., X - 1]  c[] = наддиагональ, регулярно индексируется из [0, ..., X - 2], c[X - 1] - верхний правый  cmod[], u[] = векторные царапины, каждый из которых имеет длину X  */                           /* нижний левый и верхний правый углы циклической трехдиагональной системы соответственно */ const double alpha = a [ 0 ]; const double beta = c [ X - 1 ];             /* произвольно, но выбрано так, чтобы избежать деления на ноль */ const double gamma = - b [ 0 ];      cmod [ 0 ] = альфа / ( b [ 0 ] - гамма ); u [ 0 ] = гамма / ( b [ 0 ] - гамма ); x [ 0 ] /= ( b [ 0 ] - гамма );                   /* цикл от 1 до X - 2 включительно */ 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 ; }                                                      /* обработка X - 1 */ const double m = 1.0 / ( b [ X - 1 ] - альфа * бета / гамма - бета * cmod [ X - 2 ]); u [ X - 1 ] = ( альфа - a [ X - 1 ] * u [ X - 2 ]) * m ; x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;                                                      /* цикл от X - 2 до 0 включительно */ 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 ] * бета / гамма ) / ( 1,0 + u [ 0 ] + u [ X - 1 ] * бета / гамма );                         /* цикл от 0 до X - 1 включительно */ for ( int ix = 0 ; ix < X ; ix ++ ) x [ ix ] -= fact * u [ ix ]; }

Существует также другой способ решения слегка возмущенной формы трехдиагональной системы, рассмотренной выше. [5] Рассмотрим две вспомогательные линейные системы размерности : ( n − 1 ) × ( n − 1 ) {\displaystyle (n-1)\times (n-1)}           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 … a n u n − 1 + b n u n = d n . i = 4 , … , n − 1           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 … a n v n − 1 + b n v n = − c n . i = 4 , … , n − 1 {\displaystyle {\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 {\displaystyle u_{1}=0} v 1 = 1 {\displaystyle v_{1}=1} { u 2 , u 3 … , u n } {\displaystyle \{u_{2},u_{3}\dots ,u_{n}\}} { v 2 , v 3 … , v n } {\displaystyle \{v_{2},v_{3}\dots ,v_{n}\}}

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

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

Таким образом, мы находим: x 1 = d 1 − a 1 u n − c 1 u 2 b 1 + a 1 v n + c 1 v 2 {\displaystyle x_{1}={\frac {d_{1}-a_{1}u_{n}-c_{1}u_{2}}{b_{1}+a_{1}v_{n}+c_{1}v_{2}}}}

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

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 ]) { /* сначала решим систему длины X - 1 для двух правых частей, игнорируя ix == 0 */ cmod [ 1 ] = c [ 1 ] / b [ 1 ]; v [ 1 ] = - a [ 1 ] / b [ 1 ]; x [ 1 ] = x [ 1 ] / b [ 1 ];                                          /* цикл от 2 до X - 1 включительно */ 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 ; }                                                      /* обработка 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 ;                                                           /* цикл от X - 2 до 1 включительно */ 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 ]);                         /* цикл от 1 до X - 1 включительно */ for ( int ix = 1 ; ix < X ; ix ++ ) x [ ix ] += x [ 0 ] * v [ ix ]; }

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

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

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

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

Для подробного рассмотрения параллельных трехдиагональных и блочных трехдиагональных решателей см. [9]

Ссылки

В Wikibook Algorithm Implementation есть страница по теме: Алгоритм трехдиагональной матрицы
  1. ^ Прадип Нийоги (2006). Введение в вычислительную гидродинамику . Pearson Education India. стр. 76. ISBN 978-81-7758-764-7.
  2. ^ ab Biswa Nath Datta (2010). Численная линейная алгебра и приложения, второе издание . SIAM. стр. 162. ISBN 978-0-89871-765-5.
  3. ^ Николас Дж. Хайэм (2002). Точность и устойчивость численных алгоритмов: Второе издание . SIAM. стр. 175. ISBN 978-0-89871-802-7.
  4. ^ Батиста, Милан; Ибрагим Каравия, Абдель Рахман А. (2009). «Использование формулы Шермана–Моррисона–Вудбери для решения циклических блочно-трехдиагональных и циклических блочно-пятидиагональных линейных систем уравнений». Прикладная математика и вычисления . 210 (2): 558–563. doi :10.1016/j.amc.2009.01.003. ISSN  0096-3003.
  5. ^ Рябенький, Виктор С.; Цынков, Семен В. (2006-11-02), "Введение", Теоретическое введение в численный анализ , Chapman and Hall/CRC, стр. 1–19, doi :10.1201/9781420011166-1, ISBN 978-0-429-14339-7, получено 2022-05-25
  6. ^ Квартерони, Альфио ; Сакко, Риккардо; Салери, Фаусто (2007). «Раздел 3.8». Численная математика . Спрингер, Нью-Йорк. ISBN 978-3-540-34658-6.
  7. ^ Chang, L.-W.; Hwu, W,-M. (2014). "Руководство по реализации трехдиагональных решателей на графических процессорах". В V. Kidratenko (ред.). Численные вычисления на графических процессорах . Springer. ISBN 978-3-319-06548-9.{{cite conference}}: CS1 maint: multiple names: authors list (link)
  8. ^ Венетис, IE; Курис, А.; Собчик А.; Галлопулос, Э.; Самех, А. (2015). «Прямой трехдиагональный решатель, основанный на вращении Гивенса для архитектур графических процессоров». Параллельные вычисления . 49 : 101–116. doi :10.1016/j.parco.2015.03.008.
  9. ^ Галлопулос, Э.; Филипп, Б.; Самех, А.Х. (2016). "Глава 5". Параллелизм в матричных вычислениях . Springer. ISBN 978-94-017-7188-7.
  • Конте, SD; де Бур, C. (1972). Элементарный численный анализ . McGraw-Hill, Нью-Йорк. ISBN 0070124469.
  • В данной статье использован текст из статьи Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) на CFD-Wiki, которая находится под лицензией GFDL .
  • Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Раздел 2.4". Numerical Recipes: The Art of Scientific Computing (3-е изд.). Нью-Йорк: Cambridge University Press. ISBN 978-0-521-88068-8.
Retrieved from "https://en.wikipedia.org/w/index.php?title=Tridiagonal_matrix_algorithm&oldid=1217946723"