Введение в вычислительную физику - Федоренко Р.П.
ISBN 5-7417-0002-0
Скачать (прямая ссылка):
В итоге мы получаем уравнения для ф и р:
ф = д sin2 ф — (Mp) cos2 ф, ф(0) = а,, (22)
P= Р(1/р+ ф) cos2 ф/sin2 ф + / sin ф, P(O) = ^1. (23)
Деление на sin ф не приводит к неприятностям, так как из (22) следует
(Mp+ ф) = (q + I/р) sin2 ф.
Уравнение для (3 лучше использовать в виде
P = / sin ф + р sin ф cos ф (q + 1 /р). (24)
Заметим, что нас интересуют задачи с большим параметром. Большой величиной считаем q, причем знак q(t) не определяем однозначно. При д> О однородное уравнение (17) имеет решения экспоненциального типа (одно быстро убывающее, другое быстро растущее вправо). При д< О уравнение (17) имеет решения колебательного характера, причем, если \д\Т^>\ (1000, например), на 10, Т] укладываются десятки колебаний. В некоторых задачах (см. § 15) g(t) может иметь разные знаки на разных частях интервала [О, Г].
Что можно сказать в этих условиях об уравнении (22) для ф и о возможности достаточно аккуратного численного интегрирования его на большом интервале времени? А интервал действительно большой, так как величина
T \д Sin2 ф + (l/р) cos2 ф]^ « 2Tg sin ф cos ф
в принципе может быть достаточно большой, если только функция Ф(і) не «застревает» надолго в окрестности таких значений, где sin ф COS ф « 0.
9 — 1833
258
ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ
[Ч. II
Пытаясь разобраться в характере уравнения (22), заметим, что ІФІ ^tf- (Мы пользуемся не очень строгими оценками; речь идет о
грубом анализе, в котором величиной (1/р) cos2 ip O(I) можно пренебречь.) Следовательно, все траектории (22) проходят в конусе с раствором, определяемым величиной |tf|, т.е. в целом траектории (22) не имеют экспоненциального роста, траекторий типа e^\t, среди решений (23) нет, хотя на малых частях интервала [О, T] такие решения могут и появиться.
В § 7 специально подчеркивалось, что процесс накопления погрешностей при численном интегрировании задачи Коши существенным образом зависит от устойчивости вычисляемой траектории. Устойчивость же определяется уравнением в вариациях, которое для (22) имеет вид (бір — малое возмущение траектории <р)
6ф = sin 2ip (tf — l/р) 6<р. (25)
Устойчивость траектории зависит от знака и модуля q sin 2ір. В принципе возможна ситуация, когда sin 2ip ^ 1 и траектория сильно
неустойчива: решение (25) ведет себя, как е^ч\'К Однако из уравнения (22) видно, что <р быстро уходит от такого значения и «неустойчивый» участок на траектории не может быть длительным. Ограничимся этими простыми соображениями, которые, видимо, можно превратить в достаточно строгий анализ.
Периодическая прогонка. Опишем полезный в приложениях алгоритм решения специальной системы уравнений высокого порядка, возникающей при решении краевой задачи для уравнения Штурма—JIиувилля (10.1) с периодическими краевыми условиями
X(O) = X(T), Jc(O) = Jc(T).
Разностное уравнение (10.2) можно считать определенным при всех, значениях и = 0, 1, ..., N, если реализовать условие периодичности, отождествив выходящие за пределы сетки значения с сеточными:
xN +1 = Х0> X~l ~ xN'
Итак, мы приходим к системе уравнений, аналогичной (10.3):
aOxN ^OxO cOxI = f 0’ anxn-1 ^nXll + CnXn + l = /„, aNxN-1 — ^NxN + CNX 0 = f N’
где п = 1, 2, ..., N — 1. Матрица системы отличается от знакомой нам, трехдиагональной наличием двух ненулевых элементов: на последнем месте первой строки И на первом месте последней. Для ЭКОНОМНОГО! (требующего O(N) операций) решения такой системы А. А. Абрамо-,
ЖЕСТКИЕ ЛИНЕЙНЫЕ КРАЕВЫЕ ЗАДАЧИ
259
вым построен алгоритм, обобщающий классическую прогонку. Он часто используется в практической вычислительной работе.
Решение ищем в форме «прогоночного соотношения»
*„-1 ~ PпХп + Qn + ^nxN-
Очевидно, первое уравнение системы можно записать в этой форме, и мы получаем явные выражения для стартовых значений прогоноч-ных коэффициентов:
P К = C(/^0> Qn — f </^0» = а(/^0*
Теперь стандартная процедура позволяет получить рекуррентную формулу. Пусть Pn, Qn, Rn известны. Исключая из уравнения с индексом п значение имеем уравнение
ап(РпХх + Qn + KXn) - Mn + СпХп +1 = /„.
связывающее хп, хп+1, xN+l. Этому уравнению можно придать стандартную прогоночную форму, разрешив его относительно хп. В результате мы получаем искомые соотношения:
А = Ьп - апРп>
„ _сп п _“А л _anQn-fn
n + l А ’ п+1 — А ’ Уп+1— А
Эту операцию можно продолжить вплоть до значения п = N — 1. Прогоночное соотношение
xN-I = PNXN Qn ^NxN
после подстановки в N-Q уравнение системы даст соотношение, связывающее Xn с X0. Придадим ему форму Xn = CtivX0 + рд, и будем искать решение в виде хп = апх0 + р„.
Новые прогоночные коэффициенты ап, р„ (п = N, N — 1, ..., 1) находим по рекуррентным формулам справа-налево, имея их стартовые значения aN, fiN. Для этого из прогоночного соотношения *„-i = PnXn + Qn + RnxN' считая, что значения ап и Pn известны, исключим хп, xN:
Хп і = Рп(апхо + Pn) + Qn + KiaNxо + Pjv)*
Приводя подобные члены, получаем рекуррентные соотношения «»-1 = Рп«п + KaN' Pn-I = Qn+ PnK + KPn-Последнее такое соотношение имеет вид
