1626435587-51311eae4652e8ad616b5bdef025cbb3 (844239), страница 18
Текст из файла (страница 18)
е. вместо величины ℎ ( ) вычисляя площадь трапеции [ ( ) + ( + ℎ )] ℎ /2. Тот же порядок точности (Σ = (ℎ2 ))можно получить, аппроксимируя интеграл площадью прямоугольников с высотой, равной значению функции посередине между соседнимиузлами сетки, ℎ ( + 21 ℎ ). Эти же идеи лежат в основеиметодов Эйлера, которые мы сейчас рассмотрим.состоит в построении численного решения задачи Коши (70) по схеме[︂(︁)︁]︂ℎ ( , ) + +1 , + ℎ ( , ) .(77)+1 = +2го модифицированногоИсправленный метод Эйлераисправленно-бГеометрический смысл разностной схемы (77) представлен на рис. 7 ( ).Из точки ( , ) выполняется шаг ℎ методом Эйлера (73) по касательной к интегральной линии уравнения (70), в результате чего мы приходим в точку ( +ℎ, +ℎ ( , )) ≡ (+1 , ˜+1 ), см. рис.
7 ( ). В полученной точке (+1 , ˜+1 ) вычисляется правая часть уравнения (70)и находится среднее арифметическое между ( , ) и (+1 , ˜+1 )аналогично тому, как мы вычисляли полусумму значений функциина правом и левом концах отрезка при вычислении интеграла методомтрапеций. После этого делается шаг ℎ из точки ( , ) вдоль прямой,имеющей наклон 21 [ ( , ) + (+1 , ˜+1 )], что приводит нас к следующему ( + 1)-му узлу сетки, в точку (+1 , +1 ).
На рис. 7 ( )видно, что ошибка, допускаемая на шаге численного интегрирования91ббисправленным методом Эйлера (77), значительно меньше аналогичнойошибки в схеме ломаных (73).Определим порядок точности исправленного метода Эйлера (77).Для этого вычислим невязку , подставляя точное решение задачи(70) в численную схему исправленного метода Эйлера (77), предварительно переписав её в виде, соответствующем дифференциальномууравнению (70):[︂(︁)︁]︂1+1 − = ( , ) + +1 , + ℎ ( , ) .ℎ2Для невязки имеем:[︂(︁)︁]︂+1 − 1 =− ( , ) + +1 , + ℎ ( , ) ,ℎ2где, как и раньше, используется обозначение ≡ ( ).
Подставляяв полученное равенство разложение () в ряд Тейлора в точке , атакже раскладывая (+1 , + ℎ ( , )), имеем:2 =]︀ + ′ ℎ + ′′ ℎ2 + (ℎ3 ) − 1 [︀− + + ℎ + ℎ + (ℎ2 ) .ℎ2В правой части последнего равенства значения функции и её частныепроизводные вычисляются в точке ( , ). Используя разложение (71)решения () в ряд Тейлора, получаем окончательно = (ℎ2 ), т. е.исправленный метод Эйлера имеет второй порядок точности.Перейдем теперь к рассмотрению, являющегося аналогом метода центральных прямоугольниковдля вычисления интегралов (модифицированный метод Эйлера в точности совпадает с методом центральных прямоугольников для задачи(75), см.
п. 4.3 на с. 60):(︂)︂ℎℎ(78)+1 = + ℎ · + , + ( , ) .22модифицированного метода Эй-лераГеометрический смысл выражения (78) состоит в следующем. Вначалемы делаем шаг ℎ/2 из текущей точки ( , ) методом Эйлера (73),получая точку (+1/2 , +1/2 ) = ( + 21 ℎ, + 21 ℎ ( , )). В этойпромежуточной точке, лежащей посередине между узлами сетки, мывычисляем функцию (+1/2 , +1/2 ) и используем полученное значение вместо ( , ) для выполнения шага в методе Эйлера (73).92Определим порядок точности модифицированного метода Эйлера(78). Для невязки имеем:(︂)︂+1 − ℎℎ =− + , + ( , ) ,ℎ222]︁ + ′ ℎ + ′′ ℎ2 + (ℎ3 ) − [︁ℎℎ− + + +(ℎ2 ) = (ℎ2 ).ℎ22Следовательно, модифицированный метод Эйлера также имеет второйпорядок точности.
=6.3. Методы Рунге — Кутты 2-го порядкаРассмотренные выше исправленный и модифицированный методыЭйлера являются представителями. Покажем, как оба этих метода могут быть получены без помощи аналогий с задачей о вычислении интегралов, и укажемнаиболее точную схему в семействе методов Рунге — Кутты второго порядка. Будем искать численное решение задачи Коши (70) в виде[︁(︀)︀]︁+1 = + ℎ · 1 (, ) + 2 + ℎ1 , + ℎ2 (, ) ,(79)ты второго порядкасемейства методов Рунге — Кут-где 1 , 2 , 1 , 2 — некоторые коэффициенты, ℎ ≡ +1 − .
Подставляяв (79) вместо точное решение (), получаем невязку: =+1 − − 1 ( , ) − 2 ( + ℎ1 , + ℎ2 ).ℎЗдесь и далее мы обозначаем ( , ) за для краткости. Раскладывая в ряд функцию двух переменных ( + ℎ1 , + ℎ2 ) в точке( , ), имеем для невязки:(︂)︂1ℎ2ℎ34 = + ′ ℎ + ′′+ ′′′+(ℎ)− −ℎ26[︂]︂− 1 − 2 + ℎ1 + ℎ2 −[︂]︂ℎ2 21 + 2 1 2 + 2 22 + (ℎ3 ).(80)− 22Подставляя разложение (72) точного решения () с производными(71) и приравнивая нулю коэффициенты при ℎ0 , ℎ1 и ℎ1 , по93лучаем систему из⎧⎨ 1 − 1 − 21− 2 · 1⎩ 212 − 2 · 2трёх уравнений:===000⇒2 ≡ , 1 = 1 − ,1.1 = 2 = 2(81)При выполнении условий (81) невязка (80) занулится в нулевом и первом порядках по шагу сетки ℎ — схема (79) будет иметь второй порядокточности при любых .
С учетом (81) можно переписать схему (79) ввиде:[︂(︂)︂]︂ℎℎ, + ( , ) . (82)+1 = + ℎ (1 − ) (, ) + +22однопараметрическое семейство методов Рунге — Кутты второго порядка. В частности, при = 21 выражеТаким образом, мы получилиние (82) даёт исправленный метод Эйлера (77), при = 1 — модифицированный метод Эйлера (78).Зададимся теперь вопросом, какое значение параметра обеспечивает наиболее высокую точность расчётов в методах Рунге — Куттывторого порядка (82)? Для ответа на этот вопрос выпишем невязку (80)при выполнении условий (81):[︂]︂ℎ222 = + 2 + + + −6[︂]︂121ℎ22 + 2 + 2 + (ℎ3 ).−2 4244Группируя в полученном выражении частные производные и выносяобщий множитель, приходим к ответу:(︂)︂ [︂]︂ℎ232 =1− + 2 + +64[︂]︂ℎ22 + + (ℎ3 ).(83)+6Из полученного выражения (83) видно, что при любых невязка необращается тождественно в ноль во втором порядке по шагу сетки ℎ,однако при = 3/4 зануляется первая группа слагаемых, пропорциональная ( + 2 + 2 ).
Таким образом, все схемы семейства(82) имеют второй порядок точности, однако из (83) можно ожидать,что для минимизации численного коэффициента при ℎ2 в остаточномчлене в общем случае следует использовать = 3/4.94В частных случаях = 21 и = 1 (исправленный и модифицированный метод Эйлера соответственно) коэффициент в круглой скобкеперед первым слагаемым в (83) равен − 21 и 14 соответственно, что позволяет ожидать, что модифицированный метод Эйлера в общем случаебудет давать несколько меньшую погрешность, чем исправленный метод Эйлера (для квадратурных формул отношение погрешностей метода трапеций к погрешности метода центральных прямоугольниковбыло равно 2 + (ℎ)).Заметим, что в случае решения задачи Коши (75) с правой частью,зависящей только от (не зависящей от ), выражение для невязки(83) зануляется во втором порядке точности при выборе = 43 , чтопозволяет вычислять интегралы по схеме (82) с погрешностью (ℎ3 ):∫︁() =0−1∑︁[︂ℎ=0(︂)︂]︂(︁)︁321( ) + + ℎ+ max ℎ3 ,443(84)где ℎ ≡ +1 − .
Заметим, что формула (84) занимает промежуточноеместо между обобщённой формулой трапеций (Σ = (ℎ2 )) и формулой Гаусса при = 2 (Σ = (ℎ4 )).6.4. Метод Рунге — Кутты 4-го порядкаИспользованный выше подход — суммирование значений функции (, ), вычисленных в промежуточных узлах сетки и умноженных направильным образом подобранные коэффициенты, — позволяет повысить порядок точности интегрирования ОДУ. Рассмотрим без доказательства способ построения численного решения (70) с четвёртым порядком точности.требует вычисленияправой части уравнения (70) в четырёх точках на каждом шаге:[︂]︂ℎ1 + 22 + 23 + 4 ,+1 = +6Метод Рунге — Кутты четвёртого порядка1234= ( , ),(︂)︂ℎℎ= + , + 1 ,22(︂)︂ℎℎ= + , + 2 ,22= ( + ℎ, + ℎ3 ).95(85)Можно заметить, что в частном случае (75), т.
е. при (, ) = (),имеем 2 = 3 , и схема Рунге — Кутты 4-го порядка (85) переходит вметод Симпсона для вычисления интегралов.Схема (85) Рунге — Кутты 4-го порядка реализована в большом количестве программ и пакетов, выполняющих численное интегрирование ОДУ и их систем. Методы Рунге — Кутты более высокого порядкаточности практически не употребляются. Пятичленные формулы имеют четвёртый порядок точности; шестичленные имеют шестой порядокточности, но слишком громоздки. Кроме того, высокий порядок точности схем реализуется лишь при наличии у правой части ОДУ (, )непрерывных производных соответствующего порядка.6.5. Использование адаптивного шагаВ различных физических задачах возникают дифференциальныеуравнения, правая часть (, ) которых может существенно (на порядок величины и более) изменяться на промежутке интегрирования[, ].
Например, при вычислении траектории движения кометы гравитационные силы будут возрастать обратно пропорционально квадрату расстояния при приближении к Солнцу и планетам и убывать привыходе на периферию Солнечной системы. Для повышения скоростирасчётов естественно использовать переменный шаг сетки, увеличиваяего там, где решение меняется медленно, и уменьшая в тех областях,где правая часть (70) принимает относительно большие значения.96Одной из наиболее простых и притом достаточно универсальныхвозможностей для автоматической коррекции шага численного интегрирования является пересчёт на сетке с удвоенным шагом.