§ 17. Принципы построения стандартных программ с автоматическим выбором шага

Как отмечалось в § 12, интегралы от функций с особенностями типа хорошо вычисляются методами интегрирования с переменным шагом, если узлы интегрирования распределены оптимальным образом. По-видимому, столь же хорошо будут интегрироваться функции с особенностями других типов. Поэтому представляется заманчивым строить стандартные программы численного интегрирования так, чтобы для любой функции распределение узлов являлось оптимальным или близким к нему.

В § 11 указана возможность распределения узлов, близкого к оптимальному, после исследования поведения подынтегральной функции на редкой сетке. Однако в случае резко меняющихся функций, например функций типа , реализация этой возможности приводит к недостаточно удовлетворительным результатам.

Поэтому при разработке стандартных программ интегрирования приняты несколько другие процедуры распределения узлов интегрирования, обеспечивающие лучшее приближение к оптимальному распределению узлов для функций с особенностями. Рассмотрим некоторые из них.

Для вычисления интегралов по элементарным отрезкам разбиения выбираются: квадратурная формула

и мера погрешности

Пусть вычисляется

Первая процедура, которую естественно назвать горизонтальной, определяется заданием параметров , и . Полагаем . Предположим, что каким-то образом уже вычислено приближенное чение интеграла . Программа располагает в каждый момент времени некоторым значением , с которым надо начинать считать оставшуюся часть интеграла. Вычисляем величину , соответствующую отрезку . Если оказалось, что , то вычисляем приближенное значение по формуле (1) и полагаем . Мы получили приближенное значение величины . В случае, когда , полагаем , в противном случае полагаем . Мы готовы к следующему шагу. Если оказалось, что , то принимаем за новое значение величины и возвращаемся к исходной позиции: вычислено приближенное значение интеграла и задан шаг .

Начальные условия для применения процедуры:

Процедура должна также иметь блок окончания работы: если оказалось, что , то следует положить . Установилась практика брать .

Другая процедура, которую можно назвать вертикальной, определяется заданием числа и заключается в следующем. Пусть на каком-то шаге возникает необходимость вычисления интеграла по отрезку разбиения , т.е. . Вычисляется величина , соответствующая этому отрезку. Если она оказалась меньше , то этот интеграл вычисляется по соответствующей формуле (1) и программа переходит к следующему справа отрезку разбиения. В противном случае отрезки и объявляются отрезками разбиения и программа сначала обращается к вычислению интеграла по левому из этих отрезков. В начале работы программа обращается к вычислению исходного интеграла .

Распределение узлов, осуществляемое этими процедурами, не является асимптотически оптимальным в смысле, определяемом в § 11, по следующим причинам: величина не является главным членом погрешности квадратуры (1), а, как правило, есть некоторая грубая, завышенная по порядку оценка для него; отрезки для него могут принимать липхь довольно редкий ряд значений (при этом для первой процедуры они имеют вид , а для второй .

Рассмотрим некоторые моменты, связанные с практическим использованием описанных выше процедур.

Чтобы в отдельных случаях сделать распределение узлов более близким к оптимальному, иногда величины , соответствующие отрезкам разбиения, сравнивают не с , а с , где специально подбирается. На первоначальном этапе применения ЭВМ сложилась традиция производить сравнение величины с величиной , т. е. брать . Оказалось, что для резко меняющихся функций получающееся распределение узлов было далеко от оптимального, а в случае функций с разрывами программа не могла вычислить интеграл. В самом деле, пусть, например,

Если при каком-то q величина содержит одновременно значения в узлах, где и , то имеет порядок , поэтому нет оснований надеяться что неравенство будет выполняться при малых программа будет осуществлять дробление тага до машинного нуля. Если функция меняется на малом участке очень резко, то шаг будет дробиться неоправданно сильно. Это обстоятельство было отмечено пользователями, и после теоретического анализа, в результате которого и была решена задача об оптимальном распределении узлов, было принято полагать если нет особых причин против этого. В случае резко меняющихся функций следует иметь в виду, что программа может не заметить участка резкого изменения функции. Пусть когда степени и пусть реальная подынтегральная функция есть

где практически отлична от нуля лишь на малом отрезке, например

Рис. 3.17.1.

(рис. 3.17.1). Для определенности обратимся к первой из описанных выше процедур. Если отрезок разбиения удален от точки С, то и программа будет увеличивать шаг. При подходе к точке С шаг может оказаться настолько большим, что окрестность, где функция существенно больше 0, окажется заключенной между узлами. Тогда будет близко к 0 и на отрезке, содержащем точку С. В результате в качестве значения мы получим значение погрешность этого приближенного значения близка к 1. Для контроля над точностью можно было бы попытаться произвести интегрирование с другим значением , однако с большой вероятностью получился бы тот же результат.

Не следует думать, что этот недостаток свойствен только методам интегрирования с автоматическим выбором шага. Если производить вычисление этого интеграла по формулам с постоянным шагом , то при все может оказаться, что во всех узлах интегрирования, и мы получим приближенное значение . Производя численное интегрирование при нескольких шагах и оценивая погрешность по правилу Рунге, можно прийти к неправильному выводу, что приближенное значение интеграла . В целом можно сказать, что в случае использования алгоритмов интегрирования с автоматическим выбором шага возможность получения подобных неправильных выводов несколько больше.

При составлении стандартных программ всегда приходится балансировать между двумя крайностями: гарантией требуемой точности для любой подынтетральной функции или быстрым вычислением интеграла с нужной точностью для большинства предъявляемых к решению задач.

По-видимому, можно говорить о каком-то принципе неопределенности по отношению к методам высокой эффективности: дальнейшее повышение скорости работы должно сопровождаться уменьшением надежности. Чтобы избежать «проскакивания» областей резкого изменения функции, можно предусмотреть в программе наличие отрезков, где шаг интегрирования дробится принудительно. Например, в описание стандартной программы можно включить концы некоторого отрезка . Если отрезок разбиения , пересекается с , то следующий отрезок выбирается по более сложному правилу. В качестве следует задавать некоторый отрезок, где подынтегральная функция резко меняется.

Обратим внимание на ряд дополнительных моментов. Пусть для гладкой функции

Для понимания механизма изменения шага предположим, что — кусочно-постоянная функция. В области постоянства имеем

Зададимся некоторым числом . Рассмотрим случай . Пусть для шага , с которым мы приступаем к интегрированию в данной области постоянства производной, выполняется соотношение . Тогда шаг будет увеличиваться каждый раз в , пока не станет больше . Поскольку при шаге в раз меньшем, выполнялось соотношение , то интегрирование будет осуществляться с наименьшим шагом , для которого . Если же для начального в этой области шага выполняется соотношение , то шаг станет наибольшим из шагов , для которых выполняется условие . Если, как мы предположили, , то может оказаться, что в зависимости от поведения функции при меньших в рассматриваемой области выбираются различные шаги. Мы видели, что шаг интегрирования желательно распределять так, чтобы оценка погрешности на всех отрезках разбиения интервала интегрирования была примерно одинаковой. Исходя из сказанного мы делаем вывод о желательности употребления . Рассмотрим случай . Можно указать шаг такой, что

Если оказалось, что

то при интегрировании по рассматриваемому отрезку будет выбран именно такой шаг. Однако при

программа выработает следующий шаг, равный . Для этого шага, вследствие (6), окажется, что , поэтому шаг будет признан непригодным. Программа возьмет шаг, равный поскольку для него выполняется условие (6), то интегрирование по текущему отрезку будет закончено и за исходный шаг для дальнейшего счета будет принят шаг . Таким образом, фактическое интегрирование происходит с шагом h и на каждом шаге делается попытка произвести интегрирование с шагом . Если вычисления с шагами h и независимы, то на каждом шаге объем работы будет удваиваться. Известен следующий экспериментальный факт: если мы зададимся произвольным и из различных независимых источников получим некоторые числа , то величины будут асимптотически равномерно распределены на отрезке [0, 1].

Задача 1. Считая, что - числа, происходящие из независимых случайных источников, показать, что при математическое ожидание затрат работ пропорционально . (Отсюда следует целесообразность выбора ).

Наши рассуждения проведены в предположении, что , и поэтому требуют экспериментальной проверки. Исходя из результатов такой проверки на практике обычно берется равным или несколько большим, чем , например, в одной из распространенных стандартных программ .

Как и в случае интегрирования с постоянным шагом, возникает вопрос практической оценки погрешности. Обратимся ко второй процедуре, имеющей более простое описание. Как и ранее, проведем исследование на примере кусочно-постоянной производной . Шаг интегрирования в области постоянства для этой процедуры определяется условием: это наибольший шаг h, вида , для которого выполнено условие

Очевидно, что для такого максимального шага справедливо соотношение

Произведем интегрирование с некоторым . В области, где

«старый» шаг интегрирования будет уже непригоден, поэтому произойдет дробление шага не менее чем вдвое. В области, где выполняется условие

дробления шага не произойдет. Таким образом, при может случиться, что дробление шага интегрирования произойдет лишь на части отрезка интегрирования или его вообще не будет.

Из сравнения результатов расчетов с и мы не получим информации о величине погрешности интегрирования в той области, где дробления шага не произошло. При условии

шаг дробится всюду, однако отсюда не следует, что сравнение результатов расчетов с и любым , удовлетворяющим (8), дает надежную гарантию точности результата.

Задача 2. Пусть . Подобрать постоянные с такие, что при

приближенные значения и интеграла, соответствующие и , будут одинаковыми, а погрешность отлична от нуля.

Если непрерывна, то при и при введении в алгоритм условия уменьшения максимального шага, пропорционального можно предложить и обосновать аналог правила Рунге оценки погрешности

здесь — порядок главного (при ) члена погрешности квадратуры (1).

Литература

1. Крылов В. П., Бобков В. В., Монастырный П. И. Начала теории вычислительных методов. Интерполирование и интегрирование. — Минск: Наука и техника, 1983.

2. Крылов В. И., Шульгина А. Т. Справочная книга по численному интегрированию. — М.: Наука, 1966.

3. Крылов В. И., Бобков В. В., Монастырный П. И. Вычислительные методы. Т.1. — М.: Наука, 1976.

4. Мысовских И. П. Интерполяционные кубатурные формулы. — М.: Наука, 1981.

5. Никифоров А. Ф., Суслов С. К., Уваров В. Б. Классические ортогональные полиномы дискретной переменной. — М.: Наука, 1985.

6. Никифоров А. Ф., Уваров В. Б. Специальные функции. — М.: Наука, 1979.

7. Никольский С. М. Квадратурные формулы. — М.: Наука, 1979.

8. Stroud А. Н. and Secrest D. Gaussian Quadrature Formulas. — Englewood Cliffs, N. Y.: Prentice-Hall, 1966.