Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 26
Текст из файла (страница 26)
ЗАДАЧА КОШИ Г46 зации экстраполяцнонных методов, однако, продолжается, н потому они должны учитываться при любом сравнении, 6.8. Подпрограмма йКЕ45 ККГ45 — это подпрограмма решения задач Коши для обык- новенных дифференциальных уравнений, основанная на формулах Рунге — Кутга, которые предложил в 1970 г.
Фельберг; она была составлена в 1974 г. Шампэнем и Уотсом. Подпрограмма' требует шести вычислений функции за шаг. Четыре из этих зна чений берутся с одним набором коэффициентов, приводя к методу четвертого порядка,, н все шесть значений комбинируются с дру-' гими коэффициентами, что дает метод пятого порядка.
Сравне- ние двух вычисленных величин позволяет получить оценку ошиб- ки, которая используется для контроля длины шага. Мы не станем подробно выводить полные формулы, а рас-: смотрим вместо этого более простой метод второго порядка. Этот метод использует лишь два вычисления функции за шаг. Первое из 'ннх .Ь,=йи)(у„, ти).
Затем предпринимается дробный шаг, основанный на й,. Введем два пока неопределенных коэффициента я и р. Пусть л,=п„((у„+фйь 1„+ал„). Наконец, для полного шага берется комбинация двух значений функции с еще двумя неопределенными коэффициентами: уии-1=уи+ г|Л т Гид Чтобы определить коэффициенты, разложим й, и й, в ряды Тейлора (с двумя переменными) относительно точки (1и, у„), что дает /г, =- л„)„, й.=". (7.+(1Мя..+ М ..+ ) У...=У.+(Ъ+У)й.г.+МИ„,.И+уи йи7а.+"..
Разложение у„„можно теперь сравнить с разложением точного локального решения и„(1), определенного в 9 6.3: и„(1„„) = и„(1„) + й„и,', (1„) + + ии ((„) +... А2 уи+ йи(и+ 2 (7у, и~и+ (~ и) + ' ' ' Приравнивая коэффициенты прп одинаковых степенях й„полу чаем три уравнения относительно четырех неизвестных коэффи 6.6. ПОДПРОГРАММА РКРг6 147 циентов: 1 7 а=в 2 Выбрав один из этих коэффициентов, например а, за параметр, приходим к однопараметрическому семейству методов Рунге — Кутта: й! М (ул 1в) йз = Ь„7" (у„+ ай„1„+ ай„), у„+,—— у„+(1 — — ) й,+ — й,.
1 Два очевидных выбора для а — это — и 1. Они определяют ме- тоды, тесно связанные с квадратурными формуламн соответст- венно прямоугольников и грапеций. Исследование первых чле- нов, отброшенных в вышеуказанных разложениях, показывает, что ни при каком выборе а не удается исключить эти члены для всех функций Г(у, 1); таким образом, метод имеет второй порядок при любом а. Формулы, используемые в подпрограмме ККР45, предпола- гают шесть вычислений функции за шаг: й„й„..., й,.
Они опре- деляются соотношениями ! -! йг=й„|(У„+ ч~'„Р;,йг, 1„+а,й„), 1=1, ..., б. 7=1 Новое значение УРР! вычисляется затем как взвешенная комби- нация шести величин /г;: 6 УР+ ! = УР + Х 7А. г=! Всего здесь 27 коэффициентов: б коэффициентов а;, 15 — 1)гг и 6 †Заметьте, что коэффициенты (3ы образуют нижний треугольный массив, так что каждое й, получается по предыдущим значениям й,.
Чтобы определить коэффициенты, все й! раскладываются в ряды Тейлора, эти разложения подставляются в формулу для у„+! и результат сравнивается с рядом Тейлора для и„(1„,). Можно найти такие значения коэффициентов, что разность между двумя разложениями начинается с члена порядка Ь„', а не Ь„'. Получающийся метод имеет, следовательно, пятый порядок. Хотя это и может показаться удивительным, но в то время как для р=1, 2, 3 и 4 можно построить метод порядка р за р вычисле- а ЗАДАчл Коши 148 таблица 6.! ний функции, для в=5 или 6 р вычислений функции дают лишь метод порядка р — 1. (Это отчасти объясняет популярность классического метода четвертого порядка, упомянутого в предыдущем параграфе; чтобы получить еще один дополнительный порядок точности, требуются два дополнительных вычисления функции.) Возможны различные способы выбора коэффициентов, приводящие к методам пятого порядка, Значения, предложенные Фельбергом, указаны в табл.
6.!. «Лишнее> вычисление функции, однако, также не пропадает. Фельберг нашел еще один набор коэффициентов у,', такой, что четыре из шести коэффициентов л; можно скомбинировать и получить второе значение у,*„,, точное до четвертого порядка. Таким образом, Коэффициенты у,' также приведены в таблице. Программа в действительности не вычисляет у,'„„а находит вместо этого оценку ошибки ~(у; — у,')йь используемую для контроля величины ;=1 шага. Иллюстрирующая основная программа в конце этого параграфа с помощью подпрограммы гсКГ45 определяет движение двух тел под действием взаимного гравитационного притяжения.
Пусть х(г) и у(г) обозначают положение одного тела в координатной системе, начало которой зафиксировано в другом теле. Из ньютоновых законов движения вытекают следующие дифференци- б.б. Подпгогплммл Ркгзб 149 альные уравнения: — еззх (О х" (1) = — азу (1) у (()=- где Р(1) =[х(1)'+у(1)'~" "-, а бх — константа, зависящая от гравитационной постоянной, масс обоих тел н выбранных единиц измерения. Если взять начальные условия х(0) =-1 — е, х'(0) =-О, у(0) =О, у (0) =а( — ~ уз х уз у у,=х, у,=у, Уравнения и начальные условия приобретают теперь вид (Уз+ У1)"' у,(О) =-1 — е, у, (0) =- О, уз =- уз Уз =Уз у Уз з== Ч у,(о) =-о, Уз Е!9 е~з/з Уз-= — —, у,(0) =а (— (! — е) Нормируя 1, можно исключить ц„но мы не делали этого, желая проиллюстрировать использование фортранной конструкции СОММОХ для передачи параметров типа а ат основной программы к подпрограмме, определяющей уравнения.
Параметр (ГБАО является важной контрольной переменной. При первом обращении к Й КГ45 ему следует присвоить значение 1. Обычно 11 КГ45 изменяет его значение на 2, и при последующих обращениях нужно сохранить это значение. Значения, не равные 2, полученные на выходе ГхКГ45, сигнализируют о наличии различных нерегулярностей или ошибок, подробно описываемых в комментариях. !ГБАО=-4 и 1Г(.АЙ=-7 — зто предупреждения, что программе будет очень трудно получить требуе- где е — некоторый параметр, такой, что 0 е(1, то решение оказывается периодическим с периодом 2п7а.
Орбита будет эллипсом с эксцентриситетом е, один фокус которого находится в начале координат. Чтобы записать это как систему из четырех уравнений первого порядка, положим Е. ЗАДАЧА КОШИ Таблица 6.2 Результаты, полученные иллюстрирующей программой. мую точность. Пользователь может продолжать процесс, либо увеличить границы погрешностей, либо перейти к подпрограмме многошагового метода. 1Н АО=З указывает, что затребована слишком высокая относительная точность, а 1Г).АО=-5 или б означают, что прежде, чем продолжать, нужно изменить допуски на ошибку. 1Г1.АС«=8 — указание 'на неправильность вызова тсКР45.
Мы настоятельно советуем пользователю включить в свою основную программу проверку параметра !г'1.АО '). В данном иллюстративном счете было взято е=0.25 и а=и(4; были напечатаны координаты тела для 0(1з-12 с шагом 0.5. Эти результаты приведены в табл. 5.2. Отметьте, что орбита периоднчна и период равен 8. ') Отметим уже здесь, что часто встречающиеся в комментариях к подпрограмме ЙКР46 термины «чнсло вычислений производных» нлн «чнсло значений функции» следует понимать как число обращений н подпрограмме !одновременного) вычисления правых частей системы.— Прим.
иерее. 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 0.7%000000 0.619768032 0.294417538 — 0.105!76382 — 0.490299793 — 0.8!3942832 — 1.0540315!7 -1.20073%42 -1.250000001 — 1.200735042 — 1.054031517 -0.813942832 -'0.490299793 — 0.105!76383 0.294417537 0.619768031 0.749999996 0.619768024 0.294417526 †' 0.105176395 -0.490299806 — 0.813942843 - 1.054031524 -1.200735047 -1.250000002 О.ОООРРОООО 0.477791373 0.812178519 0.958038092 0.939874996 0.799590802 0.575706078 0.300160708 -0.000000001 — 0.300160709 -0.575706079 -0.799590803 — 0.939874996 -0.958038092 -0.812178518 -0.477791370 ОЛЮ0000006 0.411791319 0.812178522 0.958038091 0.939874991 0 199590194 0.575706068 0.300! 60697 — 0.000000011 152 Ф ЗАДАЧА КОИИ с с с с с С С с с с с с с С С С С С С С С с С с с С С С с С С С |, с с С с с С с с с С с с с с 11 3! 41 71 81 ГОЙМАТ (Г5.1, 2Г16.9] ГОРМЛТ (17Н ТО!.ЕКАМСЕБ КЕБЕТ, 2Е |2.3) ГОКМАТ (1! Н МАМУ 5ТЕРБ) ГОЙМЛТ (12Н М(|СН ОПТРЦТ) ГОЙМАТ (14Н |МРКОРЕК СА1.Ц ЕМО БЦВКО()т!МЕ РКГ46 (Г, МЕЯМ, У, Т, ТОУТ, РЕЕЕКК, ! ЛВБГЙК,!Г1.АО, %ОРК, !%0ЙК) МЕТОД РУНГŠ— КУПА — ФЕЛЬБЕРГА ЧЕТВЕРТОГО-ПЯТОГО ПОРЯДКА СОСТАВИТЕЛИ ПРОГРАММЫ вЂ” Н.
А %АТТБ, Е. Г. БНЛМР1МЕ БЛМО(Л |ЛВОКАТОР1ЕБ АЕВ(|ЯЦЕКЯОЕ, МЕ% МЕХ1СО К КГ45 ПРЕДНАЗНЛЧГНЛ ГЛЛВНЫМ ОБРЛЗОМ ДЛЯ РЕШЕНИЯ НЕЖЕСТКИХ И СЛАБО ЖЕСТКИХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ, КОГДА ВЫЧИСЛЕНИЕ ПРОИЗВОДНЫХ НЕ СЛИШКОМ ДОРОГОСТОЯЩЕЕ йКГ45, ВООБЩГ ГОВОРЯ, НЕ СЛЕДУЕТ ИСПОЛЬЗОВАТЬ, ЕСЛИ ПОЛЬЗОВАТЕЛЮ ТРЕБУЕТСЯ ВЫСОКАЯ ТОЧНОСТЬ РЕЗЮМЕ ПОДПРОГРАММА РКГ45 ИНТЕГРИРУЕТ СИСТЕМУ ИЗ МЕЯМ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ Г! ЕРВОГО ПОРЯДКЛ СЛЕДУЮЩЕГО ВИДА. Оу(!):ОТ =Г(т, у(1). у(2), ..., у(меями, ГДГ У(1) ЗАДАНЫ В Т. ОБЫЧНО ПОДПРОГРЛММУ ПРИМЕНЯЮТ ДЛЯ ИНТГГРИРО- влния от т до тоцт, одн Ако ее л(ожно использоВАть И КАК ОДНОШЛГОВ ЫЙ ИНТЕГРАТОР, ЧТОБЫ ПРОДОЛЖИТЬ РЕШЕНИЕ НЛ ОДИН ШЛГ В НАГРАВЛГНИИ ТО(|Т.
НЛ ВЪ| ХОДЕ ПАРЛМЕТГ'АМ, ФИГУР!ЛРУЮЩИМ В СПИСКЕ ВЫЗОВА, ПРИСВЛИВЛ|ОТСЯ ЗНАЧЕНИЯ, НЕОБХОДИМЫЕ ДЛЯ ПРОДОЛЖЕПГ|Я ИНТЕГРИРОВАНИЯ. ПОЛЬЗОВАТЕЛК) ЙУЖНО ЛИШЬ ЕЩЕ РЛЗ ОБРАТИТЬСЯ К ЙКГ46 (И, ВОЗМОЖ1!О, ОПРЕДЕЛИТЬ НОВОЕ ЗНАЧЕНИЕ ДЛЯ Т01|Т). В ДЕЙСТВИТЕЛЬНОСТИ ЙКГ45 — ЗТО ПРОГРЛММЛ ИНТЕР. ФЕЙСА, КОТОРЛЯ ВЫЗЫВАЕТ ПОДПРОГРАММУ ККГБ, ОСУШЕСТВЛЯЮЦ[УЮ ПРОЦЕСС РЕЦ|ЕНИЯ. ЙКГ5 В СВОЮ ОЧЕРЕДЬ ВЬ!ЗЪ|ВЛЕТ ПОДПРОГРЛММУ ГЕН1., КОТОРАЯ ВЫЧИСЛЯЕТ ПРИБЛИЖЕННОЕ РГШЕНИЕ НЛ ОДИН ШЛГ. ККГ45 ИЯПОЛЬЗУЕТ МЕТОД РУНГŠ— КУТТЛ вЂ” ФЕЛЬБЕРГА, ОПИСАННЫЙ В СЛГДУ|01ДЕЙ ПУБЛИКЛЦИИ: Е. ГЕНЕВЕКО, 1..0%-0йРЕй С1.Л551СЛ1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
