Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Для систем уравнений с якобиевой матрицей был разработан специальный метод прогонки, требующий O(N) операций и O(N) ячеек памяти. Этот метод был разработан почти одновременно в нескольких местах учеными, занимавшимися в сущности одной и той же проблемой. Она была связана с закрытыми работами, поэтому публикации последовали спустя годы после того, как была придумана и весьма эффективно применена прогонка.
Алгоритм прогонки. Решение ищется в форме прогоночного соотношения:
хп— і ^nXп + Qn’ ^ I > 2, ..., N,
где Pn, Qn — неизвестные пока прогоночкые коэффициенты. Нетрудно видеть, что, определив Pn, Qn, мы в сущности приведем систему с трехдиагональной матрицей к системе с двухдиагональной матрицей.
ПРОГОНКА В РАЗНОСТНОЙ ЗАДАЧЕ ШТУРМА-ЛИУВИЛЛЯ
95
Алгоритм начинается с того, что левое краевое условие записывается в форме прогоночного соотношения:
X0 = (C0Zb0) X1 Cl0Zb0,
т.е. P1 = C0Zb0 Q1 = -Cl0Zb0 (отметим, что P1 < 1). Далее следует процесс прямой прогонки: последовательно (по рекуррентным формулам) вычисляются P2, Q2, затем P3, Q3, и т.д. вплоть до PN, Qn.
Выведем эти рекуррентные формулы. Пусть прогоночные коэффициенты Pn, Qn уже вычислены. Подставляя xn_1 = PnXn + Qn в n-е уравнение anxn_l — bnxn + Cnxn+1 = dn, получаем
“п(рп*п + - ЬпХп + СЛ+1 = dn-
Соотношение между хп и хп+1 представим в форме прогоночного соотношения, т.е. разрешим его относительно хп:
с a Q — d
п . , Я^/1 п
b -a P ’ b -a P *
п п » п п п
Оно примет стандартный вид хп = Pn+Jxn+1 -h Qn+1, если положить
с а О —d
П — пХП я
P =--------------— о
п + 1 b —a P '
b -a P ’ ^n-H b -a P '
п п п п п п
Это и есть рекуррентные формулы прямой прогонки. По ним вычисляем Pn, Qn вплоть до PN, Qn. (Для этой цели мы последний раз можем использовать стандартное трехчленное уравнение с номером N-I.)
Имея неиспользованное пока правое краевое условие (N-e уравнение) aNxN_i — bNxN = dN и прогоночное соотношение xN_l = = PnXn + Qn, можно найти величину xN (разрешение правого краевого условия).
Неизвестные хп последовательно определяются справа-налево по формулам (обратная прогонка)
Хп-1 = PnXn + Qn’ n = N,N- 1,...,1.
Исследование устойчивости прогонки. Проведем исследование вычислительной устойчивости прогонки, т.е. покажем, что погрешности вычислений, связанные с конечной разрядностью машинных чисел (погрешности округления), и погрешности, связанные с машинной реализацией арифметических операций, накапливаются в такой мере, что это не приводит к существенным погрешностям в результате.
96
ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ
14.1
Рассмотрим сначала процесс прямой прогонки коэффициента Pn. Введем величины Pn, вычисленные по формулам прогонки Pn = cj(bn — апРп) в идеальной арифметике, т.е. без погрешностей округления и погрешностей в выполнении операций. Реальные расчеты на ЭВМ дают Pn = Pn + Ьп, где &п — погрешность предшествующих получению Pn вычислений.
Предположим, что дп — малая погрешность (такова она, во всяком случае, на начальном этапе прогонки) и проанализируем процесс ее эволюции, записывая соотношение между дп и 6п+1 в виде
сп
Pn+i + 6n+l = b -а (Р +6 ) + гп-
п п п п'
Здесь Zn — суммарная погрешность, связанная с выполнением операций в правой части формулы. В нее включаются также погрешности машинного представления коэффициентов с„, Ьп, ап, погрешность выполнения машинных операций и погрешность округления, связанная с записью Рп+1 в памяти в виде Р“+1. Последняя погрешность очень мала и зависит от разрядности представления чисел в ЭВМ.
Таким образом, погрешность 6п+1 есть следствие двух погрешностей: наследственной погрешности Ьп, в которой суммируются погрешности предшествующих вычислений, и локальной погрешности zn, для которой нетрудно указать хорошую оценку через погрешности сп, Ъп и т.д.
Воспользуемся тем, что IdnI^Pn, и применим линейный анализ, пренебрегая величинами 0(62). Используя обычное исчисление дифференциалов, получаем
Из этого соотношения следует формула
К+1=Т(Рп + і)2 К + V « = 0. 1. -- N- 1.
Заметим, что число шагов, выполняемых по рекуррентной формуле, достаточно велико (порядка Т/х), и нас будет интересовать асимптотическое поведение погрешности (при N= Т/ T-* OO).
Сначала получим важные для анализа свойства прогоночных коэффициентов.
Л емма 1. При значениях Ьп > ап 4* C11 и P1 < 1 для всех п имеем 0 < Р„ ^ 1.
ПРОГОНКА В РАЗНОСТНОЙ ЗАДАЧЕ ШТУРМА-ЛИУВИЛЛЯ
97
В самом деле, если О Pn < 1, то
а) Рп + 1 ~ Ь-аР'^Ь-а. >
П П
б Л P , = --------Sj-----------=------------- ^ 1
) п + 1 Ьп-аРп ап + сп-апРп сп + ап(1-Рп)
Разумеется была использована положительность локальных коэффициентов ап, Ьп, сп.
Лемма 2. Пусть p(t) удовлетворяет условию Липшица. Тогда I a Jcn I ^ 1 + Ст, где постоянная С не зависит от т.
Действительно,
aJl - Р«п-х12)-х4« _ , , П/ ч сп p(tn + xh)+xqn
С учетом приведенных оценок имеем
|йв + 1|«Е(1+Ст)|&в| +є, |е„| =S є.
Используя это соотношение и стандартные рассуждения (см. § 5, 8), получаем
