Федоренко Р.П. Введение в вычислительную физику (1185915), страница 49
Текст из файла (страница 49)
Если использовать схему интерполяционного типа именно с этими узлами, В-устойчивость такой схемы доказывается просто н не требуется проверять свойства сложно определяемой матрицы Д. Другими словами, В-устойчивость оказывается закономерным свойством схемы интерполяционного типа с гауссовыми узлами. Предварительно докажем одно полезное свойство схем интерполяционного типа, справедливое при любых, в сущности, узлах $У. Напомним, что в такой схеме фигурируют значения уз и ЧИСЛЕННЫЕ МЕТОДЫ МАТЕМАТИЧЕСКОЙ ФИЗИКИ 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). Здесь необходимы пояснения. Оценка (4), если не оговаривать требований к С, тривиальна: она выполняется всегда (за редким исключением задач, как говорят, «на спектре», т.е. когда нарушены условия единственности и существования решения при любых правых частях).
Однако «универсальная» оценка (4) в общем случае имеет постоянную с ехр ЦА11т). мы же рассматриваем класс задач, в котором С = О(1) ~ ехр (11А11Т). Такие жесткие краевые задачи играют осо- нгивлнженные методы вычислительной»изики 1ч. и бенно большую роль в приложениях (им присвоено специальное наименование «вычисдительно корректные задачи»). Приведем самые общие сведения из теории, позволяющей выделять эти задачи.
Оказывается, что далеко не все формально возможные постановки задач для жесткой системы (1) приводят к вычислительно корректной задаче. Определяющую роль играет такой сравнительно простой и легко контролируемый фактор, как число краевых условий на левом и правом концах интервала [О, Т1, Это число должно находиться в определенном соотношении с числом точек в разных частях спектра.
Точнее, необходимым условием вычислительной корректности краевой задачи являются следующие неравенства: а) й ~ 1, т.е. число краевых условий на левом конце интервала должно быть не меньше числа сильно убывающих вправо решений; б) и — и и!+, т,е. число краевых условий на правом конце интервала должно быть не меньше числа сильно убывающих влево решений.
При нарушении этих условий краевая задача (1) оказывается вычислительио некорректной, т,е. в оценке (3) постоянная С имеет недопустимую с принятой здесь точки зрения величину 0(е~г). Докажем этот факт. (Доказательство технически простое; оно поучительно, так как вскрывает механизм вычислительной некорректности.) Построим конечное (О(1)) возмущение решения краевой задачи, приводящее к очень малому возмущению краевых условий.
Будем истолковывать эту конструкцию «в обратном направлении» вЂ” как демонстрацию того, что очень . малое возмущение краевых условий приводит к конечному возмущению решения. Итак, пусть к(г) — решение краевой задачи (1), (2) и пусть А < г . Рассмотрим множество всех сильно убывающих вправо решений однородной системы (1).