Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 48
Текст из файла (страница 48)
Для любых двух решений х„и х„, полученных по этой схеме, имеет место соотношение 11х„+, — х„, 11 < 11х„— х„11, что сильнее требуемого свойства В-устойчивости. Если, кроме того, имеет место соотношение (52), то можно получать В-устойчивость на более широком классе систем — с положительной односторонней константой Липшица 1 = 0(1). В этом случае, очевидно, мы придем к соотношению (58) 11х„+, — х„~,)( < (1+ Ст) 11х„— х„11.
Интерполяционная схема с гауссовыми узлами. Выше был указан общий способ построения полностью неявных разностных схем типа Рунге — Кутты. Отмечалось, что выбор узлов интерполяционной формулы ~У является резервом, разумно распоряжаясь которым можно получить полезные свойства схемы. Подтвердим это положение, выбирая в качестве узлов ЦУ так называемые гауссовы узлы. Как известно, при любом выборе узлов ~~ (т'= 1, 2, ..., з) квадратурная формула точна в классе полиномов степени з — 1. Имея в своем распоряжении з свободных параметров $т, можно выбрать их так, что квадратурная формула станет точной в классе полиномов степени (з — 1) + з= 2з — 1, Такие узлы называются гауссовыми (для них существуют таблицы).
Если использовать схему интерполяционного типа именно с этими узлами, В-устойчивость такой схемы доказывается просто н не требуется проверять свойства сложно определяемой матрицы Д. Другими словами, В-устойчивость оказывается закономерным свойством схемы интерполяционного типа с гауссовыми узлами. Предварительно докажем одно полезное свойство схем интерполяционного типа, справедливое при любых, в сущности, узлах $У. Напомним, что в такой схеме фигурируют значения уз и ЧИСЛЕННЫЕ МЕТОДЫ МАТЕМАТИЧЕСКОЙ ФИЗИКИ 1Ч. и з40 Ут = У(УТ) (У = 1, 2, ..., «), определяемые решением системы уравнений (38), н интерполяционный полипом х.'(4), вычисляемый по узлам гу и значениям /т. Расширим эту информацию, добавив узел го= г„, и зададим значения Ув=х„, Ут (2 = 1, 2, ..., «).
по сеточной функции [Р, Ут) (У = О, 1, ..., «) построим интерполяционный полинам р(г) степени «: р(р)) = Уз (/= О, 1, .„«). Имеет место следующее утверждение. Утверждение. Полипом р(1) в точках 11 О=1,2,..., «) удовлетворяет исходному дифференциальному уравнению: (59) Докаэатиельстес) Имеем очевидное тождество р(г ) — р(т ) — $ р(г) ип ) Но р(Г)) = У', р(ГО) ° х„, р(1) — полипом степени « — 1; интеграл точно вычисляется по квадратурной формуле с коэффициентами та; . Итак, У' = х„+ т Ч' а„.
р(Р), 1= 1, 2, ..., «. г =! Видно, что величины р(Р) и гз =г(уз) = г"[р(гт)[ удовлетворяют одной и той же системе уравнений (38) (линейиой относительно этих величин). Предполагая невырожденность матрицы А, получаем совпадение: р(11) = У! = Л р(1!) 1, у = 1, 2, ..., «. Кроме того, очевидно, что р(1) вв.х(г), так как р н .2' — полиномы . степени « — 1, совпадающие в «точках ~Р (у =1, 2, ..., «). Таким образом, схему интерполяционного типа можно получить так называемым методом коллокации. Решение на интервале [г„, 1„+ т[ ищется (приближенно) в виде полинома степени «, для определения которого используются условия коллокации (т.е. выполнения уравнения в нескольких точках тт, число которых, естественно, связывается с числом свободных параметров полинома): р(«Ф) — Лр(е')[ 1=1 2, ...) «.
К этому добавляется, конечно, условие р(гв) = х„. $171 241 жесткие системы одг р(1.) = х. Р(1. + егт) У', р(1 + $ст) = г(У'), Такие же соотношения, естественно, выполняются и для второго решения. Рассмотрим величину г(1) = ПР(1) — Р(1)11'= (Р(1) — Р(1) Р(1) — Р(1)) Нас интересует сравнение величин г(1„) и г(1„+ т). В-устойчивость будет установлена, если будет показано (ограничимся классом дис- сипативных систем), что г(1„+ т) ц г(1„). Очевидно, с +т ~ — '„", 11 с г(1„ + т) — г(1„) = и надо оценить интеграл от г: г(1) = 2(р(1) — р(1), р(1) — р(1)). Заметим, что р и р — полиномы степени з; следовательно, г(1) — полипом степени 2г, а г(г) — полипом степени 24 — 1.
Для таких полнномов гауссова квадратура точна и интеграл ~ г сй можно вычислить по квадратурвой формуле: с +т г с(С = ~Х'„бС г(1„+ с'т) ~ ос' (Р— р, р — р) ~ с +1с,. с„ 1 с В силу (59) имеем (у = 1, 2, ..., г) г(1~+ ~1") 2(Р Р' Р Р)! с +1сс (у(Ус) — у(У1), Уг — Ус) и О. Здесь, конечно, использовалась диссипативность системы уравнений х = Дх). Так как коэффициенты гауссовых квадратур положитель- ны, пазучаем требуемый результат — атграктивность разностной схемы, Пусть теперь хв н х„— два разных решения, полученных по схеме (36), (38) с гауссовыми узлами 1'. Используем полиномы р(1) и .х(1), определенные выше для решения х„(аналогнчные полнномы для х„обозначим р(с) и х.(1)), Как оказалось, для этих полиномов точки сс явлюотся точками коллокации„т.е.
выполняются соотношения (у = 1, 2, ..., з) 242 пгинлижвнныв матовы вычислительной ьизикн 1ч. и !! 18. Жесткие линейные краевые задачи Рассмотрим численное решение специального класса линейных краевых задач, неявно содержащих большой параметр. Формально задача ставится стандартно: на интервале (О, Т! решается линейная задача — = Ах+ а(1), х Е 1т'", (1) с краевыми условиями (при г = 0 и г = Т соответственно) (1,, х(0)) = йп ю'= 1, 2, ..., й < л, (1н х(т)) = Ьп 1= й+ 1, й+ 2, ..., л, (2) 1; Е Я", !!1,!! = 1, 1= 1, 2, ..., н.
Матрицу А в теоретических рассуждениях будем считать постоянной, хотя все последующее относится н к случаю, когда А зависит от т', но изменяется со'временем в некотором смысле (ниже будет уточнено, в каком именно) «медленно», Жесткость системы (1) означает, что собственные значения матрицы А можно разделить на три насти: 1) левый жесткий спектр (собственные значения Л, ) характеризуется соотношениями ! 1ш Л, ! < 2„1 = 1, 2, ..., Г; КеЛ, < — Ь, 2) правый жесткий спектр (собственные значения Л,+.) характеризуется неравенствами «еЛ+ и +.1, !1ш Л,+ ! < Е, 1 = 1, 2, ..., 1+; 3) мягкий спектр образуют собственные значения Х,, удовлетво- Т ряющие условиям !Х.! <1, у'=1,2, ...,1, ) +)++У=и.
Число 2.П» 1 является характеристикой жесткости системы. Кроме того, мы считаем, что Т1 = 0(1), ТЕ ли 1, Качественная структура решения в этом случае хорошо известна, ее легко угадать нз общих соображений: решение содержит два пограничных слоя — левый и правый. Это следует хотя бы из вида общего решения однородной задачи х(~) =х~' С, ел 'Ф, + 2' С+е'~ ' Ф+ + 2,' с.е"' <р,, ! / где Ф,, Ф+, р — собственные векторы А, соответствующие выделенным выше частям ее спектра. При обобщении анализа на случай 243 жесткие линлйныл кгхлвыс зхдхчи З 1з1 переменной матрицы А(~) мы предполагаем, что спектр не претерпевает существенных изменений прн разных 1, т.е. число точек в каждой части спектра остается постоянным.
Область приложения численных алгоритмов, о которых пойдет речь, конечно, намного шире. В частности, точки спектра могут, так сказать, «непрерывно» заполнять интервал [-4„ 4,1, и разделение спектра на жесткую и мягкую части становится достаточно условным. Проблемы, которые рассматриваются ниже, отличаются от изложенных в 5 17. Предположим, что величина 11А11Т большая, но не очень. Использование в расчетах сетки с шагом ~1А!1 х«к1 (на практике «много меньше единицы» означает величину О.1 или 0.01, например) приемлемо с точки зрения объема необходимой памяти и количества арифметических действий. Будем считать, что именно такой шаг используется в простой разностной схеме (второго пор1щка) х+,— х х,„+х+, х Ат»нз 2 + Лт»нз (3) т ° 0,1,...,М вЂ” 1, Мх Т, Суть проблемы в том, как решать систему (2), (3). Жесткие краевые задачи возникают, например, в расчетах процессов прохождения излучения через слой большой оптической толщины.
При этом разные компоненты вектора х относятся к частицам разной энергии. Матрица А содержит члены, описывающие как поглощение излучения, так и переход частиц нз одной энергетической группы в другую («замедление»). В некоторых случаях может присутствовать и «рожденне» частиц под воздействием их поглощения («цепные реакции», характерные для процессов в ядерных реакторах).
Итак, величину 11А11Т (прнмерно совпадающую с 1Т) мы считаем не слишком большой, а вот величину еьг считаем очень большой. Характерной особенностью рассматриваемых задач является тот факт, что их решения суть ограниченные функции. Более аккуратно это можно сформулировать в виде требования ((х(~)з к С Да11 + ЙЦ), (4) где Ца11 — обычная норма правой части л(~), ~~Ь|) — норма правой части в краевых условиях. Постоянная С = О(1).