Федоренко Р.П. Введение в вычислительную физику (Федоренко Р.П. Введение в вычислительную физику.djvu), страница 15

Исследование сходимости методов Рунге-Кутты и=О л = 1, 2, ..., Ф. х — х*, о о х„— х„ —,г(х„, г„,), [Еъ(хт) ] Очевидно, теперь запись Е'(х') = 0 эквивалентна разностной зада- че. Вычисление погрешности аппроксимации состоит в том, что в уравнения С'(х') = 0 подставляегся ограничение на сетку точного решения (Я'„), Я'„= 2'(~„). Получим невязку .и'(О) — хв — — О, ю л-! " ' — У(й'„н г„$) = О(т), л=О, л=1,2,...,Ж.

Таким образом, можно сформулировать следующее утверждение. Ут ве ржд ение. Пусть решение дифференциального уравнения Я'(г) имеет две ограниченные производные (для этого достаточно, чтобы Дх, г) имела ограниченные первые производные). Тогда явная схема Эйлера аппроксимирует дифференциальную задачу и имеет первый порядок аппроксимации. Для неявной схемы (5.3) имеем, очевидно, тот же самый результат.

Несколько сложнее вопрос с третьей из простейшей схем (5.4). Она может быть оформлена в абстрактной форме: хв — хв, х, — х'„ л=1, " ' — У(х„„г„,), л= 2, 3, ..., М. (~ (х ))л Но величина х', не входит в постановку задачи, ее нужно как-то определить. В зависимости от того, как это будет сделано, порядок аппроксимации будет разный.

Применим описанную выше абстрактную схему исследования к конкретному вопросу — к обоснованию сходимости метода конечных разностей для системы обыкновеннмх дифференциальных уравнений (5.1). Мы должны доказать аппроксимацию и устойчивость, а не принимать их как предположения. Начнем с аппроксимации для простейшей схемы — явного метода Эйлера (5.2).

Прежде всего, нужно четко оформить разностные уравнения в абстрактной форме 2Р(х') = О. Под х' мы будем понимать совокупность (х„), л = О, 1, ..., М = Т/т. Оператор Е.' отображает сеточную функцию х' в такую же сеточную функцию, и надо определить все ее компоненты. Положим 71 исслвдовАнка сходимости мвтодов и нгв-кгггы Почти очевидны следующие факты для погрешности аппроксимации: "о=0 Оценим г, для нескольких способов определения х~.' а) если х", — какая угодно величина, то г, = 2'(т) — х", = О(1), схема имеет нулевой порядок аппроксимации и для расчетов непригодна; б) если х', = х~, то г, = .и'(т) — .в'(0) = О(т) и схема имеет первый порядок аппроксимации, в) если х', = хо+ тУ(хо, 1ь), то г, = Я'(г) — хо — тУ(хо, 1о) = = О(тз) и схема имеет второй порядок аппроксимации.

Разумеется, все это верно лишь в предположении, что решение дифференциальной задачи 2'(г) имеет трн ограниченных производных (для чего достаточно, чтобы Дх, ~) имела две ограниченных производных в интересуюшем нас диапазоне изменения х и 1). Подчеркнем, что обеспечение второго порядка аппроксимации третьей схемы потребовало достаточно аккуратного (хотя н не очень сложного в данном случае) определения дополнительных начальных данных, тг Это типичное обстоятельство для схем, формальный порядок которых превышает формальный порядок дифференциального уравнения. Такие схемы используются именно с целью получить более высокий порядок аппроксимации, и недостаточно внимательное решение вопроса о дополнительных начальных данных может лишить схему желаемого порядка аппроксимации.

Устойчивость разиостных схем. Рассмотрим схемы следующей формы: "=Р(к„), п=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.