Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 21
Текст из файла (страница 21)
перев. 5.5. подпРОГРАммл яовхсв 0.0349 О.Ю77 -0.0327 0.3702 -ОЛ402 0.3702 -О.О327 О Ю77 0.0349 рациональные выражения для действительных весов (умноженные на 8) включены в саму подпрограмму. Поскольку формула точна для 1(х)=1, то веса удовлетворяют равенству в ~ч'., ш» =1. в=о Однако некоторые из весов отрицательны, и в ~~.", ! ш„! = 1.4512. «=о Этот фактор важен при анализе ошибок. Например, предположим, что функция 1(х)=-1 известна лишь с некоторой ошибкой, так что 1(х)=0.99 при х=2/8, 4/8, б/8, в то время как в прочих узлах 1(х)=!.01.
Тогда формула в качестве значения интеграла даст число 1.014512. Ошибки, содержавшиеся в значениях функции„слегка увеличились в результате. Коэффициент увеличения несуществен в данной формуле, однако для формул Ньютона— Котеса, основанных на полиномах более высокой степени, он становится неприемлемым. Мы выбрали формулу 8-го порядка с тем, чтобы длина каждого элементарного отрезка получалась из длины Ь вЂ” а основного интервала делением на степень двойки.
Поэтому на двоичной машине при вычислении узлов х; не будет ошибок округления. формула 8-го порядка используется для вычисления величины Р;, определенной в предыдущем параграфе. Величина Я; получается применением той же формулы к двум половинам интервала; тем самым используется 16 элементарных отрезков.
В программе величинам Р; и 1~, соответствуют переменные ЯРКЕЧ и ЯЫО%. Переменные Я(.ЕРТ и (,1м1ОНТ содержат результаты, полученные по двум половинам интервала. Поскольку основная формула точна для полиномов девятой степени, то приложйм анализ предыдущего параграфа прн р=!0. Ошибка в Яв приблизительно равна ошибке в Р;, взятой с мно- э численнОе интеггиРОвлние !!б жителем 1/(2" — 1)=1/1023, при условии, что подынтегральная функция гладкая и члены более высокого порядка можно игнорировать. Как было указано в 2 5.4, результат, полученный для данного подынтервала, приемлем, если ! Ь! В программе ( ()АХС8 этот тест выполняется сравнением двух величин ЕВТЕКК = АВ5(! 1(ч0% — ОРКЕЕВ)/1023 ТО1ЕКК =- АМАХ1(АВЬЕКК, КЕ) ЕКК АВБ(АКЕА))« (ЬТЕР/5ТОХЕ).
Отношение ЯТЕР / ЯТОМ)Е равно л!/(6 — а), а АКЕА есть оценка интеграла по всему интервалу, так что е есть максимум из АВЯЕКК и КЕ1.ЕКК» ~)/Е Подынтервал приемлем, если ЕБТЕКК.).Е.ТО1ЕКК. Если подынтервал принимается, то ! 1ЫОЮ добавляется к выходному параметру КЕМ).Т, а ЕВТЕКК добавляется к выходному параметру ЕККЕБТ. Кроме того, величина (ЯрО%в ЯРКЕЪ)/1023 прибавляется к внутренней переменной СОК11, поскольку 9!+ (1)! — Р;)/1023 обычно является более точной оценкой для /и чем само Я;. В действительности это не что иное, как один шаг экстраполяции Ричардсона.
В то время как формулы Р! и Я! по отдельности точны для полиномов степени 9, комбинация этих двух формул дает точное значение интегралов от полиномов степени 11. Окончательное значение переменной СОК11 перед самым выходом из программы прибавляется к КЕБ(/1Т. При каждом делении подынтервала узлы и значения функции для правой половины запоминаются для последуюшего использования. Поскольку максимальный объем памяти для хранения этих величин должен быть указан заранее, то устанавливается предел (.Е'»/МАХ=-30 на уровень деления пополам. Когда этот максимум достигнут, то подынтервал принимается, даже если оценка ошибки слишком велика; однако ведется счет числа таких подынтервалов, и это число содержится в выходной информации в качестве целой части параметра Г).АО. Длина каждого из этих подынтервалов очень мала — она равна (6 — а)/2»», так что для большинства подынтегральных функций несколько подобных подынтервалов могут быть отражены в конечном результате без серьезного ухудшения точности.
Однако если подынтегральная функция имеет разрывы или неограниченные производные или «зашумлена» ошибками округлений, то предельный уровень деления пополам может достигать- вБ. подпРОГРАммз оохксз 1!7 ся очень часто. Поэтому установлен другой предел ХОМАХ на количество вычислений функции. Когда оказывается, что этот предел может быть превышен, параметр (.ЕЧМАХ уменьшается до (.ЕЧО()Т стем, чтобыоставшаяся часть интервала могла быть обработана с меньшим, чем КОМАХ, количеством вычислений функции. Кроме того, точка ХО, где встретилась трудность, отмечается и величина ( — ХО)7( — А) содержится в выходной информации в качестве дробной части параметра Е АО. Это та часть интервала, которую предстоит обработать с меньшим пределом на уровень деления пополам.
Так как ХО==  — (дробная часть Р).АО)+( — А), то положение трудного места можно определить. Программа Я()АХС8 основана на кусочно-полиномиальной аппроксимации и потому не предназначена для вычисления некоторых типов интегралов. Грубо говоря, это интегралы от функций 7(х), для которых некоторая производная )'м(х), й' 10, не ограничена или не существует. Например, интегралы вида 1 при а) — 1 определены корректно, но если а не является целым числом, то Я1)АМС8 может потребовать большого количества вычислений функции и не получить хорошую оценку ошибки.
Любой из границ погрешностей мпжет быть задано нулевое значение; однако если обе равны нулю, то, за исключением очень простых случаев, будет достигнуто предельное количество вычислений функции. Приводимая ниже иллюстрирующая программа вычисляет 8!(2), где 5!(г) есть интегральный синус: 8! (2) =- ) — с)х. а Описание ЕХТЕКМА1. в основной программе, относящееся к подпрограмме Е()г), существенно. Поскольку основная программа транслируется отдельно от подпрограмм, то нет никакого другого указания на то, что Р1)!ч' не является простой вещественной переменной.
Данная задача очень проста, и программе ЯОАМС8 требуется лишь ЗЗ значения функции. Результат, полученный на !ВМ 360 с использованием удвоенной точности, таков: КЕБАБ).Т= !.6054129768 ЕККЕВТ=0.23 Т)-!6. 5. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ 1!8 БОВЙОпт1ме ООАмсз(РОМ, А, В, АВБейй, ке1ейй, йебиьт, ЕК Р ЕБТ, НОГИН, РГАЛИ) РЕАЛ Р())5 А, В, АВБЕЙК, ЙЕЬЕКК, КЕБ()ЬТ, ЕККЕБТ, ЕЬАО 1МТЕОЕК А.*ОРич) С С С Г. С С С С С ОЦЕНИТЬ ИНТЕГРАЛ ДЛЯ ГИЯ(Х) ОТ А ДО В С ЗАДАННОЙ ПОЛЬЗОВАТЕЛЕМ ТОЧНОСТЬЮ. АВТОМАТИЧЕСКАЯ АДАПТИВНАЯ ПРОГРАММА, ОСНОВАН- НАЯ НА ФОРМУЛЕ НЬЮТОНА — КОТЕСА 8-ГО ПОРЯДКА ВХОДНАЯ ИНФОРМАЦИЯ..
Р()м ИМЯ ПОДПРОГРАММЫ-ФУНКЦИИ Р())Ч(Х), РЕАЛИЗУЮЩЕЙ ПОДЫНТЕГРАЛЬНУЮ ФУНКЦИЮ ') То есть не вынолннлсн критерий точности.— Прим, нерее. Если мы изменим задачу на о то поведение программы будет совершенно иным. Подынтегральная функция имеет неинтегрируемую особенность при х=п/2ж ж1:57. (;Ц)А)ЧС8 выдает значение 91.21 для параметра Л.А6. Это показывает, что для 91 интервала не было сходимости" и оставалось обработать 0.21 от всего интервала интегрирования, когда было обнаружено трудное место. Так как 2 — 0.21и2=1.58, то особенность легко находится. С ИЛЛЮСТРИРУЮЩАЯ ПРОГРАММА ДЛЯ О()АМС8 С КЕА1.
РОМСТ(О)4 Р()К(Х) КЕАЬ Х 1Р (Х .ЕО. 0.0) Р(1М= 1.0 1Р (Х .!(Е. 0.0) Р()М =Б1М(Х)(Х ЙЕТ()КМ ЕЬО С ехтей!(Аь Рич КЕАЬ А, В, АВБЕКК, КЕЬЕКК, КЕБ()1Т, ЕККЕБТ, РЬАО 1МТЕОЕй )(ОР!)Г( А =0.0 В =2.0 РЕМЕЙК =- !.ОŠ— !О АВБЕйй =- 0.0 САНЬ О()А)(С8(РИМ, А, В, АВБЕКК, КЕЕЕЙЙ, КЕБ1)1Т, ЕййЕБТ, МОРОМ, РААС) ТУК!ТЕ(6, !) РЕЯЛ.Т, ЕККЕБТ 1Р (РЬАО .КЕ. 0 0) ТУР!ТЕ(0,2) Г1.АО 1 РОКМАТ(8Н РЕЯЗ! Т=, Г!5.10, !ОН ЕЙЙЕБТ=, Е!0.2) 2 ГОКМАТ (44Н ттАК!41)ЧО..
ЙЕБОЬТ МАТ ВЕ 1)МКЕПАВЕЕ. Р1.АО=, Рб.2) БТОР ЕИО З. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ 122 с с С с СОБРАТЬ ЭЛЕМЕНТЫ, НЕОБХОДИМЫЕ ДЛЯ СЛЕДУЮЩЕГО ИНТЕРВАЛА. ЯРК ЕН = ОЙ !ОНТ(Ь ЕН) ХО= Х(16) РО= Г(16) ОО 78 1=1, 8 Р(2*1)=Г5АНЩ), 1.ЕН) Х(2*1) = ХЗАНЕ(1, 1.ЕН) 78 СО?(Т1Н()Е ОО ТО ЗО С С ***ЭТАП 8**' ЗАКЛЮЧИТЕЛЬНЫЕ ОПЕРАЦИИ И ВЫХОД С 80 Й ВЯЛ.Т = КЕ51Л.7+ СОК!1 С С ОБЕСПЕЧИТЬ, ЧТОБЫ ЗНАЧЕНИЕ ПЕРЕМЕННОЙ ЕЙКЕ5Т С БЫЛО НЕ МЕНЬШЕ УРОВНЯ ОКРУГЛЕНИЙ. С 1Г(ЕЙЙЕЗТ .Е().
0.0)КЕТЦЙН 82 ТЕМР = АВ5(Й Е5()ЬТ)+ ЕК К ЕЗТ 1Р (ТЕМР .НЕ. АВ5(ЙЕ5()ЬТ))ЙЕТ()ЙН ЕЙКЕЗТ=2.0*ЕЙКЕЗТ ОО ТО 82 ЕНО УПРАЖНЕНИЯ 8.1. Интеграл, определяющий функцию ошибок ег1 (х)=- — ( е г г(Ь У я .1 4 н = — ох, ,) 1+ х' о найдите с помощью численного интегрирования приближения к числу и. а) Используйте формулу трапеций и формулу прямоугольников с злементарныл~и отрезками одинаковой длины Д= 1!н, взяв п=8, 32, 128. Отметьте, что ошибка приблизительна пропорциональна Де. б) Примените сплайн-квадратуры с теми же звачениями Д.
Будет ли ошибка приблизительно пропорциональна какой-либо степени й? очень просто находится с помощью численных квадратур. Напишите программу, которая, используя ОНАНС8, печатает таблицу значений функции ег((х) для х — — 0.0, 0.1, 0.2, ., !.9, 2.0. Сравните вашу таблицу с опубликованными значениями или со значениями, полученными по надежной подпрограмме, имеющейся на вашей машине. Вычисление каждого значения в вашей таблице производится независимо от прочих значений,и это не особенно эффективный способ генерярования подобных таблиц. Относительна более совершенного метода см. упр. 6.1. 6.2. Используя равенство ! УПРАЖН Е НН Я !23 в) Используйте программу О()АНС8 с различными значениями допустимой погрешности.