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

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

${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 {egin {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}} {egin {bmatrix} {x_ {1}} {x_ {2}} {x_ {3}} vdots {x_ {n}} end {bmatrix}} = {egin {bmatrix} {d_ {1}} {d_ {2}} {d_ {3}} vdots {d_ {n}} end {bmatrix}}.}$

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

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

Метод

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

${displaystyle c '_ {i} = {egin {case} {egin {array} {lcl} {cfrac {c_ {i}} {b_ {i}}} &; & i = 1 {cfrac {c_ {i} } {b_ {i} -a_ {i} c '_ {i-1}}} &; & i = 2,3, точки, n-1 end {array}} end {case}},}$

и

${displaystyle d '_ {i} = {egin {case} {egin {array} {lcl} {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, точки, n. end {массив }} конец {случаи}},}$

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

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

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

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

За ${displaystyle i = 2,3, dots, n,}$ делать

${displaystyle {egin {array} {lcl} w = {cfrac {a_ {i}} {b_ {i-1}}} b_ {i}: = b_ {i} -wc_ {i-1} d_ { i}: = d_ {i} -wd_ {i-1} end {array}}}$

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

${displaystyle {egin {array} {lcl} x_ {n} = {cfrac {d_ {n}} {b_ {n}}} x_ {i} = {cfrac {d_ {i} -c_ {i} x_ {) i + 1}} {b_ {i}}}, quad {ext {for}} i = n-1, n-2, dots, 1end {array}}}$

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

Sub TriDiagonal_Matrix_Algorithm(N%, A #(), B #(), C #(), D #(), ИКС#())    Тусклый я%, W #    За я = 2 К N        W = А(я) / B(я - 1)        B(я) = B(я) - W * C(я - 1)        D(я) = D(я) - W * D(я - 1)    Следующий я    Икс(N) = D(N) / B(N)    За я = N - 1 К 1 Шаг -1        Икс(я) = (D(я) - C(я) * Икс(я + 1)) / B(я)    Следующий яКонец Sub

Вывод

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

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

${displaystyle {egin {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} ,. конец {выравнивание}}}$

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

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

${displaystyle {ilde {a}} _ {i} = 0,}$

${displaystyle {ilde {b}} _ {1} = b_ {1},}$

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

${displaystyle {ilde {c}} _ {1} = c_ {1},}$

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

${displaystyle {ilde {d}} _ {1} = d_ {1},}$

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

Чтобы еще больше ускорить процесс решения, ${displaystyle {ilde {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 {egin {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 {egin {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}}}$

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

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

Учебник Вычислительная математика Авторы Quarteroni, Sacco и Saleri перечисляют модифицированную версию алгоритма, которая избегает некоторых делений (вместо этого используется умножение), что полезно для некоторых компьютерных архитектур.

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

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

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

Викибук Реализация алгоритма есть страница по теме: Алгоритм трехдиагональной матрицы

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

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