Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 15
Текст из файла (страница 15)
Устойчивость разиостных схем. Рассмотрим схемы следующей формы: "=Р(к„), п=0,1,...,У вЂ” 1, М=7Ут. Как мы увидим в дальнейшем, такая форма охватывает важные классы схем Рунге — Кутты. Функция Цх), конечно, связана с правой частью уравнения у(х, г) (эта связь будет конкретизирована ниже). Все дальнейшее не претерпит никаких изменений, если вместо Р(х) будут использоваться функции Р(х, г, т), но ради простоты записи мы выбросим несущественные аргументы. Во все выкладки их прн желании можно вписать механически, ничего не меняя. Проверка устойчивости связана со сравнением решений двух «почти совпадающих» систем уравнений: " = Р(х„) + з'„, х = х,'„ уо= уо. основы вычислитвльноа мАтвмлтики 1ч. ! (Таким образом мы докажем устойчивость разностных уравнений по начальным данным и правым частям.) Доказательство.
Перепишем уравнения в виде х«+1 хи+ т~(х«) + твв~ у.,ь = у„+ ту(у„) + Вычитая второе уравнение из первого, получаем (г) Цх„+, — у„„,Ц ж Цх„— у„Ц+ «ЦР(х«) — Р(у„)Ц+ 2тв. Используя условие Липшипд, преобразуем оценку Цх„, — у„+,Ц а (1 + Ст) Цх, — у„Ц + 2«в. Применим зту основную оценку последовательно: Цх, — у~ Ц % (1 + Ст) Цхо — увЦ + 2«а, Цхз — узЦ и (1+ Ст) Цх1 — уД + 2тв и и (1+ Ст)зЦх„— увЦ + 2««[1 + (1 + Ст) [, Цх — узЦ и (1 + Ст) Цх — у Ц + 2тв < и (! + Ст)зЦхо уоЦ + 2те[1+ (1+ Ст) + (1+ Ст)з]. Легко угадывается и доказывается по индукции общая формула: Цх„— у„Ц а (1 + Ст)" Цхо — увЦ + + 2тв[1+ (1 + Ст) + ... + (1 + Ст)" Относительно в'„и а'„' предположим, что онн ограничены обшей констзнтой: Цв'„Ц и в, Цв'„'Ц и в.
Еще раз подчеркнем, что на самом деле ниже речь пойдет не о двух системах, а о семействах систем с параметром т. Чем меньше т, тем больше уравнений в системах, а нас будут интересовать оценки, равномерные по т. Ради простоты мы не пишем х'„, Р(х'„, т) и т.д., но «скрытый» аргумент т следует всегда иметь в виду. Теор,ема.
Пусть функция Г(х) удовлетворяет условию Липшица с постоянной С: ЦР(х) — Р(у)Ц а СЦх — уЦ. (Заметим, что это равномерная по т оценка, С от т не должна зависеть.) Пусть шаг т мал: Ст~1. Тогда система разностных уравнений устойчива и имеет место оценка Цх„— у Ц и ест Цх,',— у'Ц+2«ее~/С, т'и жТ/т. (1) исследовАние сходимости 'методов егнге-ютгы 73 1 71 илн (после суммирования прогрессии) 11х„— у„11 ж(1+ Ст)" +2те .' < ц(1+ Ст)" ~~»,— ув1~+ — ' . (3) Заметим, что л и Тут; следовательно, (1 + Ст)" и (1 + Ст) ™ ~ ест при Се ~ 1.
Используя эту оценку в (3), получаем (1). Кстати, условие Ст и 1 не следует считать очень жестким, так как, например, 1.5 ж ев 4, а 1.1 ж ев вм Теорема об устойчивости доказана. Применим ее к некоторым широко используемым на практике схемам. Метод Эйлера с пересчетом.
Переход от известного х„к новому х„+~ делается в два этапа: а) находится значение х„~,~ — — х„+ 0.5т Дх„, 7„); б) вычисляется хв м = х„+ т Дх„+пз, 7„+ 0.5т). Эту схему можно записать в обшей форме: х„, = х„+ тР(х„, г„, т), где г(х, г, т) ея У(х+ 0.5т У(х), 1+ 0.5т). Легко проверить, что если У удовлетворяет условию Липшица (по х) с константой С, то и г удовлетворяет этому условию с несущественно большей константой С(1 + Ст).
Таким образом, схема метода Эйлера с пересчетом устойчива. Проверим порядок аппроксимации, т.е. оценим выражение Т ( и+! л) ( и' а' Используем ряд Тейлора: + ) ~+ ф +т ь' .1О(з) В силу уравнения .Ф'„= У(л'„, 7„) имеем ~"в = Ухй'в + А~ = Ук(Я'в 7в) Лй'ю 7в) + А(й'в~ 7л) Итак, С другой стороны, =У(й; 7.)+ЗУ,(й'в 1„) У(й; 7,)+$У,(Л'„7„)+ О(т'). ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ 1ч.7 74 Обьединяя оба результата, получаем ТЕ+7 = — (Я'„+, — Я'„) — Р(е „, ГВ, т) = О(тз). Таким образом, схема Эйлера с пересчетом имеет второй порядок аппроксимации и в силу теоремы б б — второй порядок точности (в предположении, что 1(х, 7) имеет две ограниченных производных и, следовательно, в."(г) — трн).
Методы Рунге-Кутты. Метод Эйлера с пересчетом является простейшим вариантом одной из наиболее распространенных в современной вычислительной практике схем численного интегрирования обыкновенных дифференциальных уравнений, объединяющей семейство методов с общим названием «методы Рунге — Кутты». Основу этих методов составляет ряд Тейлора. Связь и.'„+, —— = Я'(г„+ т) с Х„имеет форму й'.„+, =й'„+ .й'„+ — й'„+...
+ ~ й77+ О(те ). (4) Приближенное решение находится из того же выражения, но с вы- брошенным остаточным членом: Чтобы иметь вычислительную схему, нужно дать выражение для производных. В принципе здесь нет серьезных проблем: й „=1(х„, г„), й'„= 1„(х„, г„) 1(х„, г„) +1,(х„, г„) и т.д. Однако аналитические выражения начинают катастрофически усложняться в результате дифференцирования, и этот путь оказывается крайне неудобным. В методах Рунге-Кутты строится способ вычисления отрезка ряда Тейлора, требующий лишь вычисления 1(х, г) в разных точках. Вот одна из распространенных схем. ПеРехоД от (7„, х„) к (Г„+и хВ ы) начинаетса с вычислении вспомогательных величин: х, =1(х„, 1„) т, Аз=1(х„+0.5йи т„+0.5т) т, ~3 ' 1(хл + 0,5Ам ~„+ 0.5Т) т, й = 1(х + хз, г„+ т) т. Затем делается собственно шаг иитшрирования хВ м х„+ (11б)(В7+2ез+ 2йз+ В4) исслвдовхнив сходимости методов гтнгв-ялты Таким образом, один шаг требует четырехкратного вычисления правой части.
Мы не станем точно вычислять погрешность аппроксимации, это требует громоздких выкладок. Поясним, однако, основную идею. Если провести разложения х„йм Ам й по малым параметрам с достаточным числом членов и вычислить после этого выражение к„+ (1/б)(И, + 2кз+ 2/сз+ ««), то оно совпадет с отрезком ряда Тейлора (4) вплоть до членов порядка т4. Расхождение начнется только в членах О(тз). Метод Рунге — Кутты можно записать в стандартной форме: х„, — х„ " = Г(х„, /„, т), где Г(х, /, т) — суперпозиция («многоэтажнаяь) функций /(х, /).
Легко проверить, что константы Липшица для / и Г почти (с точностью до О(т)) совпадают. Что касается порядка аппроксимации, то, как уже было объяснено, он в данном случае четвертый. Существуют схемы Рунге— Кутты разных порядков, причем порядок аппроксимации равен числу вычислений правой части на один шаг процесса. Сказанного выше достаточно, чтобы сформулировать следующее утверждение. Утверждение 1. Если функция /(х, /) имеет ограниченные производные /с-го порядка (следовательно, решение .Ф'(/) производные (х + 1)-го порядка), то метод Рунге — Кутты х-го порядка имеет х-й порядок аппроксимации и к-й порядок сходимости. Однако стоит еще раз напомнить, что, кроме погрешности аппроксимации, есть еше погрешность машинного представления чисел, т.е.
фактически после подстановки в разностные уравнения машинных чисел Я«„" получим — (й'"„+, — .«'„") — Р( «'„", /„, т) = О(т~) + — !.Ф' ~ с, и при т < ()2'~ с)щ«+ц результаты численного интегрирования с уменьшением шага начнут ухудшаться, а не улучшаться. Обратим внимание на то, что в теореме об устойчивости в оценке фигурирует крайне неприятный множитель ест. Конечно, оценки очень грубы, но простые примеры показывают, что в общем случае, если, кроме условия Лившица для /, других предположений не делать, эта оценка не улучшаема. Но в важных частных случаях онз может быть намного улучшена, и зто существенно, так как часто приходится проводить численное интегрирование в ситуации, (ч.
! ОСНОВЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ ТВ когда тгесг (к — порядок аппроксимации метода) есть величина много ббльшая требуемой точности, а попытки достичь этой точности за счет уменьшения т приводят к непосильному для ЭВМ обьему вычислений. Простые примеры более точных оценок мы сейчас получим. Рассмотрим два случая. 1. Пусть определяемая численным интегрированием система такова, что матрица А(х) = — (Е„(х) + Е'(х)) в нужном нам диапа! зоне изменения х строго отрицательна, т.е. (А(х)Ц, с) ч — а(с, Ц), »( $, х (а > О). Траекторию, в окрестности которой выполняется это условие, будем называть «устойчивой».
У т в е р жд е ни е 2. Прн интегрировании устойчивой траектории методом Рунге — Кутты л-го порядка аппроксимации погрешность приближенного решения есть О(т") при всех ! > О. Утверждение будет доказано, если мы получим оценку устойчивости разностной схемы, не содержащую экспоненциального множителя типа ест. Доказан!ел»стива. Оценим сначала норму ((е+тг„)ф, (е-ьтг„)1) !!Е+ т,д !!г $ (1, Е) + гт (А1, 1) + тг(Р'„1, Л„1) =зпр ' ' " " а1 — 2та+ О(тг). При достаточно малых т можно пренебречь величиной О(тг) и пользоваться оценкой ~~Е+ тЕ„(х))~ < 1 — С,т, С! >О.
Теперь обратимся к оценке устойчивости. Как и раньше, из уравнений (2) имеем ~)х„+! — у»+!~~ И ~)(х„— у„) + т(Е(х„) — )т(у„))~~ + 2те. Оценим норму в правой части более аккуратно: Е(х) — Е(у) = Г(у+ з(х — у)) ~, «' — г(у+ з(х — у)) !(з= ~ Е„(у+ з(х — у))(х — у) в(з. исслвдовлннв скодимостн мвтодов и нгв-ютгы 1 Так как х — у= )(х — у) ~!з, то о 11(х — у) + т(Е(х) — Р(у))11 = 1 = !! ~ (Е+ тР„(у+ з(х — у)))(х — у) с(з!! ц а 1 ц ~ 11Е+ тР„!! с(з !1х — у1! ц (1 — Сг с)11х — у1!.
о В итоге мы имеем оценку 11х„~, — у„+~!! ц (1 — С~т)11х„— у„!! + 2тв. Отсюда, как н раньше, 1 — ц — С,т)" 11х„— У„!1 ц (1 — С,т)" !!хо Уо11 + 2те 1 Таким образом, 11» У 11 ц !!хо Уо11 + 2сЖ и эта оценка не ухудшается при всех п>0, Имея такую оценку устойчивости, получаем утверждение о порядке скодимости, совпадающем с порядком аппроксимации при интегрировании на сколь угодно большом интервале времени (конечно, только при интегрировании уетойчивой задачи).
2. Более интересный н тонкий результат можно получить для «не неустойчивых» систем, т.е. для систем, у которых матрица А(х) (симметричная часть Е,) только неположительна, т.е. (АС, г) ц О, 'т' $, х. В этом случае оценка 11Е+ тР„(х) 1! находится так же, как это делалось выше; 11 Е+ тр„(х) 11 ц ! + С тз. Повторяя далее оценки для решений двух разностнык уравнений, имеем !!х„тз — уам11 ~ (1+ Стт )11х„— у„11+ 2тв, откуда уже известным способом получаем и+с т')" — ! 11х„— У„11 ц (! + Стт ) !!хо Уо11+2те Стт~ основы вычислительноя млтвмдтикн 'та Пусть теперь используется метод х-го порядка аппроксимации, т.е. е = О(т~) (й в 2). Тогда, как нетрудно понять, можно положить н и 1/тз (что соответствует времени г„= нт ж С/т); при этом (1 + Сзт~)" ж ессь Таким образом, для и и Сз/т~, Сз = О(1), имеем их — у и и ес,сз ~~х — у ) в + О(тг 1) ес,с, Итак, прн интегрировании «не неустойчивой» системы методом й-го (х в 2) порядка аппроксимации с шагом т на интервале времени О < 7 и С /т численное решение имеет точность О(т~ ').
Разумеется, на конечном интервале времени О а г а т точность метода есть О(т"). В дальнейшем она понижается на один порядок. С такими задачами вычислители встречаются при расчете процессов, носящих характер «вращений», колебаний и т.п. Их приходится рассчитывать на длительных интервалах времени, так как физическое время, интересное для приложений, обычно бывает очень большим для системы в том смысле, что оно содержит большое число колебаний или оборотов. Поэтому приведенная выше оценка точности численного интегрирования очень важна. Мы рассмотрели простые варианты теорем, дающих оценки погрешности численного интегрирования существенно более точные, чем стандартные (опирающиеся только на условие Липшица для правой части).