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

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

${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{cases}{egin{array}{lcl}{cfrac {c_{i}}{b_{i}}}&;&i=1{cfrac {c_{i}}{b_{i}-a_{i}c'_{i-1}}}&;&i=2,3,dots ,n-1end{array}}end{cases}},}$

и

${displaystyle d'_{i}={egin{cases}{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,dots ,n.end{array}}end{cases}},}$

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

${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}}}_{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{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}}}$

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

${displaystyle ({mbox{equation 2}})cdot b_{1}-({mbox{equation 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-1x_{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{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}}}$

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

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

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

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

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

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

1. Прадип Нийоги (2006). Введение в вычислительную гидродинамику. Pearson Education India. п. 76. ISBN  978-81-7758-764-7.
2. Бисва Натх Датта (2010). Численная линейная алгебра и приложения, второе издание. СИАМ. п. 162. ISBN  978-0-89871-765-5.
3. Николас Дж. Хайэм (2002). Точность и стабильность численных алгоритмов: второе издание. СИАМ. п. 175. ISBN  978-0-89871-802-7.
4. Quarteroni, Alfio; Сакко, Риккардо; Салери, Фаусто (2007). «Раздел 3.8». Вычислительная математика. Спрингер, Нью-Йорк. ISBN  978-3-540-34658-6.
5. Chang, L.-W .; Hwu, W, -M. (2014). «Руководство по реализации трехдиагональных решателей на графических процессорах». В В. Кидратенко (ред.). Численные вычисления на GPU. Springer. ISBN  978-3-319-06548-9.
6. Venetis, I.E .; Курис, А .; Собчик, А .; Gallopoulos, E .; Самех, А. (2015). «Прямой трехдиагональный решатель, основанный на вращении Гивенса для архитектур GPU». Параллельные вычисления. 49: 101–116.
7. Gallopoulos, E .; Philippe, B .; Самех, А.Х. (2016). «Глава 5». Параллелизм в матричных вычислениях. Springer. ISBN  978-94-017-7188-7.

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

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