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

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

${\ 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 О (п)}$${\ Displaystyle О (п ^ {3})}$${\ displaystyle a_ {i}}$

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

Метод

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

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

а также

${\ displaystyle d '_ {i} = {\ begin {case} {\ 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 {case}}}$

Затем решение получается обратной подстановкой:

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

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

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

Для дела ${\ Displaystyle я = 2, 3, \ точки, п,}$

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

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

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

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

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

${\ displaystyle x_ {i} = {\ cfrac {d_ {i} -c_ {i} x_ {i + 1}} {b_ {i}}} \ quad {\ text {for}} i = n-1, п-2, \ точки, 1.}$

Реализация в подпрограмме VBA без сохранения векторов коэффициентов:

Sub TriDiagonal_Matrix_Algorithm(N%, A#(), B#(), C#(), D#(), X#())
    Dim i%, W#
    For i = 2 To N
        W = A(i) / B(i - 1)
        B(i) = B(i) - W * C(i - 1)
        D(i) = D(i) - W * D(i - 1)
    Next i
    X(N) = D(N) / B(N)
    For i = N - 1 To 1 Step -1
        X(i) = (D(i) - C(i) * X(i + 1)) / B(i)
    Next i
End Sub

Вывод

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

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

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

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

${\ displaystyle ({\ mbox {уравнение 2}}) \ cdot b_ {1} - ({\ mbox {уравнение 1}}) \ cdot 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 (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 ^ {th}}$${\ displaystyle n ^ {th}}$${\ displaystyle x_ {n}}$${\ Displaystyle (п-1) ^ {th}}$

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

${\ Displaystyle {\ тильда {а}} _ {я} = 0 \,}$

${\ Displaystyle {\ тильда {b}} _ {1} = b_ {1} \,}$

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

${\ Displaystyle {\ тильда {c}} _ {1} = c_ {1} \,}$

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

${\ Displaystyle {\ тильда {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 '_ {я} = 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} x_ {i} + c '_ {i} x_ {i + 1} = d' _ {i} \ qquad &; & \ i = 1, \ ldots, n -1 \\ x_ {n} = d '_ {n} \ qquad &; & \ i = n. \\\ end {array}} \,}$

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

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

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

Варианты

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

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

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

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

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

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

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

Рекомендации

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

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