Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
Вычисление точек комплексного спектра. Спектральная задача существенно осложняется в случае несамосопряженного дифференциального оператора, когда собственные значения могут оказаться комплексными. Рассмотрим несложную задачу
-il + a^L + bx=Xx, (8)
dt2 dt
с краевыми условиями х(0) = х(1) = 0. Пусть коэффициенты а, b — какие-то несложные функции или даже постоянные комплексные числа (или х — вектор, а, b — матрицы).
Исследование спектра (вообще говоря, комплексного) несколько облегчается тем, что обычно задачей численного анализа является определение точек спектра, расположенных в некоторой окрестности нуля. Размеры этой окрестности, конечно, зависят от модуля а. В области IЯI » I a I часто можно воспользоваться подходящими асимптотическими методами, т.е. изучать спектр задачи (8) как слабое возмущение хорошо изученного спектра задачи при а — 0. Эти аналитические методы достаточно эффективны при больших IЯI, но теряют точность и иногда просто непригодны в области малых | Я |.
Однако именно точки спектра с малыми | Я | нередко представляют наибольшую прикладную ценность, и численные методы в этой области удачно дополняют аналитические исследования. Таким образом, речь идет о приближенном вычислении относительно небольшого числа точек спектра. И здесь в зависимости от ситуации может быть использован метод «пристрелки», когда интегрируется задача Коши с начальными данными х(0, Я) = 0, х(0, Я) — 1 и этим алгоритмом численного интегрирования определяется функция комплексного переменного Ф(Я) = х(1, Я).
В другой ситуации может оказаться целесообразным применение Метода прогонки, когда, например, решение ищется в форме x(t) = a(t) x(t), а для «прогоночного коэффициента» a(t) обычным способом получается уравнение (типа уравнения Риккати), содержащее параметр Я. Интегрируя эту задачу (конечно, численно), вычисляем a(t, Я) и определяем Ф(А) = а(1, Я). Если на правом конце за-
§ 15]
СПЕКТРАЛЬНАЯ ЗАДАЧА ШТУРМА-ЛИУВИЛЛЯ
187
дано более общее краевое условие (например, х(1) + Р*(1) =0), оно вместе с прогоночным соотношением х(1) = а(1) X(I) образует систему линейных уравнений относительно х(1) и х(1), а Ф(Х) определяется как детерминант этой системы. Так или иначе мы получаем функцию комплексного переменного, значения которой вычисляются процедурой интегрирования задачи Коши (в комплексных числах). Точки спектра исходной задачи суть нули этой функции.
При определенных условиях, которые хорошо изучены в теории обыкновенных дифференциальных уравнений, Ф(Х) — аналитическая функция (это следствие простой формы зависимости уравнения от X), не имеющая полюсов при ограниченных X. Поэтому для подсчета числа ее нулей в некоторой области G следует вычислить известный интеграл по контуру или, проще говоря, вращение векторного поля (Re Ф(Х), Im Ф(Х)} при обходе контура 3G. Конечно, это — громоздкая операция: надо «покрыть» контур сеткой точек (XJ и в каждой точке Xi вычислить Ф(А.-), т.е. проинтегрировать задачу Коши.
Несколько облегчает работу то, что сетка (XiJ не должна быть особенно густой. Точнее, дело обстоит так. Расчет начинается с достаточно широкой области G, имеющей (для простоты и определенности) форму прямоугольника. Вычислив вращение вдоль контура (оно будет равно 2 лп, где п — число точек спектра в G), делят его пополам (пополам делится та сторона прямоугольника, которая на данном этапе процесса локализации корней Ф(Х) длиннее).
Вычислив вращение вдоль контура одного из полученных меньших прямоугольников, определяют число точек спектра в двух частях исходной области, и т.д. Когда прямоугольник велик (и его контур достаточно длинен), шаг сетки на контуре может быть взят достаточно большим. Например, в некоторых расчетах, проводившихся по этой схеме, было принято считать «нормальной» ситуацию, в которой при переходе от Xi к Xi+1 значения arg Ф(Х) изменялись в пределах интервала [л/6, л/3]. Если изменение было меньшим, шаг увеличивался, если ббльшим, — происходил возврат в точку Xi, шаг А изменения X уменьшался (например, Д — А/2) и делался переход в точку Xi+1 = Xi + Д на контуре.
Таким образом, сетка (X1) не задавалась заранее, а «генерировалась» простым алгоритмом с адаптацией (с регулированием шага в зависимости от градиента функции arg Ф(А.)). Конечно, такая тактика сопряжена с некоторым риском: при вычислении Ф(Х^) и Ф(ХІ+1) приращение аргументов определено с точностью до 2ки (к — любое целое), причем число к вычислитель назначает сам.
Поясним сказанное подробнее. Вычислив комплексное число х + іу, обычно обращаются к подпрограмме, входящей в стандартную
188
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
библиотеку и имеющей на языке FORTRAN имя ATAN2(*, у). Результатом является главное значение arctg(y/x); к нему из каких-то дополнительных соображений нужно добавить 2 к л.
В рассматриваемом случае, используя предположение о том, что шаг А = Xj+1 — Xi «достаточно мал», число к выбирают таким, чтобы изменение arg Ф(А/+1) по сравнению с arg Ф(А^) было минимальным. Здесь, конечно, есть риск ошибиться на 2кп. Можно получить достаточно надежное подтверждение правильности этого решения (или обнаружить его ошибочность), «пройдя» участок ( Xj> Xi+l) с существенно меньшим шагом, HO это слишком «дорого» и не делается без достаточных к тому оснований.
