Гарантированная точность вычисления многомерных интегралов

Введение
Быстродействие и память электронных вычислительных машин, как это отмечается во многих современных работах, за каждую пятилетку возрастают примерно в десять раз, и вместе с этим ростом исследованиям при помощи и посредством массированных компьютерных вычислений подвергаются все более и более сложные задачи естествознания и техники.
Объем перерабатываемой при этом информации становится поистине гигантским и традиционные, хорошо себя зарекомендовавшие для относительно малых объемов обрабатываемых данных, способы контроля за происходящими в компьютере вычислениями в изменившихся условиях свою эффективность утрачивают.
По-видимому, именно по этой причине в научных изданиях регулярно появляются статьи, в заголовки которых снова и снова выносится вопрос о том, можно ли и насколько можно доверять результатам компьютерных вычислений?, (см., например, [10], [27], [353]). Ответ на этот вопрос интересен еще и тем, что позволяет избежать ненужных затрат и вероятных потерь, многократно увеличивающихся по мере усложнения рассматриваемой научно-технической проблемы.
Желание полностью контролировать сложный вычислительный процесс породило в современной математике ряд новых направлений и одно из них связано с выработкой неклассичсских критериев качества рассматриваемого процесса. Существо такого рода критериев состоит, как правило, в определении сопутствующих выбранному вычислительному процессу числовых параметров (одного или нескольких), эффективно определяемых по исходным данным задачи и позволяющих (в зависимости от собственной величины) давать гарантированные заключения о близости или удаленности получаемого машиной числа и истинного результата. Наибольшее развитие указанный подход получил во второй половине прошлого века в
Введение
3
применении к задачам линейной алгебры (см., например, [93], [79]).
Успехи, достигнутые в этом направлении, побуждают к проведению аналогичной точки зрения в областях вычислительной математики, непосредственно с линейной алгеброй не связанных, и, в частности, к выработке соответствующих критериев в теории приближенного многомерного интегрирования.
1°. Практика современных теории приближений и численных методов такова, что наиболее развиты и широко используемы способы вычисления интегралов на основе метода кубатурных (квадратурных) формул. В связи с компьютерной реализацией этого метода, т.е. в связи с его реализацией в арифметике с конечной точностью, возникает вопрос о правомочности применения оценок погрешности, получаемых в рамках функционального подхода теории кубатурных формул, к оценке погрешности реального вычислительного процесса.
Особенность реализации кубатурной формулы в арифметике с конечной точностью состоит прежде всего в том, что веса формулы и узловые значения подынтегральной функции приходится представлять в системе счисления с заданным основанием как вещественные числа с плавающей точкой, причем порядки этих представлений лежат в наперед заданном и зависящем от компьютера фиксированном отрезке числовой оси, а мантиссы этих же представлений содержат одно и то же также наперед заданное конечное число значащих цифр.
К погрешности аппроксимации, возникающей в результате замены интеграла конечной суммой взвешенных узловых значений подынтегральной функции, неминуемо добавляются при этом погрешности, обусловленные как неточным вводом в компьютер начальных данных задачи (весов формулы и узловых значений подынтегральной функции), так и неточным же выполнением сопутствующих формуле арифметических операций (сложений и умножений).
Верхняя граница возникающей таким образом полной ошибки растет пропорционально числу N узлов рассматриваемой формулы. Тем самым, начиная с некоторого значения Лго, эта верхняя граница заведомо превзойдет теоретическую погрешность кубатурной формулы, т.е. погрешность, вычисленную в предположении, что начальные данные задачи введены в
Введение
4
компьютер абсолютно точно и арифметические действия над ними совершаются также абсолютно точно. Следовательно, гарантировать точность в практических вычислениях интегралов без скрупулезного анализа сопутствующей этим вычислениям суммарной погрешности немыслимо.
Появление суммарной погрешности метода принято связывать с тремя основными видами ошибок:
ошибками, возникающими при замене операции взятия интеграла конечным числом арифметических операций (ошибки ограничения);
ошибками, содержащимися в исходной информации (ошибки ввода),
ошибками, возникающими в результате необходимости представлять результаты арифметических операций в виде цифрового вектора конечной длины (ошибки округления).
Суммарная погрешность может быть как абсолютной, так и относительной, но практика такова, что большинство пользователей предпочитают алгоритмы, в которых малость погрешности понимается по отношению к той или иной норме подынтегральной функции. Этому же предпочтению следует и настоящая диссертация.
Хорошо известно, что в классической теории ошибка кубатурной формулы на классе функций характеризуется с помощью нормы ее функционала погрешности. В случае арифметики с конечной точностью использование для оценки погрешности только лишь упомянутой нормы явно недостаточно: в ней никак не учтены ошибки ввода и округления. Чтобы учесть их, в диссертации вместо нормы функционала погрешности предлагается использовать родственное понятие “гарантированного радиуса” кубатурной формулы. В пределе, когда точность используемой в методе арифметики неограниченно возрастает, гарантированный радиус формулы переходит в норму ее же функционала погрешности.
2°. Квадратурные формулы для вычисления интеграла от непрерывной функции стали применяться и содержательно исследоваться задолго до того, как собственно появилось само определение интеграла. К примеру, еще Архимед, рассматривая задачу вычисления площади параболического сегмента, применил для этой цели формулу трапеций и, по существу, доказал сходимость соответствующего квадратурного процесса [97, с. 240]. В дальнейшем в создание и развитие теории приближенного интегрирования
Введение
5
свой заметный вклад внесли многие и многие видные математики. Геометрические традиции метода квадратур (Архимед, Кавальери, Торричелли) продолжил Б. Паскаль [97, с. 257]. П. Ферма изобрел квадратурную формулу для вычисления площади, заключенной между параболой и конечным отрезком оси абсцисс [97, с. 254]. Начиная с конца XVII века в работах основоположников и создателей анализа бесконечно малых интегрирование сформировалось как важнейшее понятие математического анализа.
Параллельно развитию теории интегрирования создавался и совершенствовался аналогичный по своему назначению аппарат численного анализа (аппарат именно аналогичный, но не тождественный). Основными объектами аналогичных интегрированию вычислительных технологий служили разного рода приближенные формулы, среди которых особое место занимали квадратурные: как результат эволюции этой математической ветви возникли и получили широкую известность квадратурные формулы Ньютона — Котеса, Эйлера — Маклорена, Грегори, Гаусса, Симпсона, Чебышева, Маркова и др. Однако все эти формулы служили для приближения исключительно одномерных интегралов, и только в 1887 г. Дж.К. Максвелл опубликовал, по-видимому, первую работу [350] по многомерному интегрированию, в которой он методом неопределенных параметров получил кубатурные формулы для квадрата и куба, точные на полиномах седьмой степени. Надо отметить, что вплоть до сороковых годов XX века, число работ по многомерному интегрированию оставалось весьма незначительным [338, с. 481], но затем ситуация в силу ряда причин резко изменилась.
К концу прошлого XX века теория кубатурных формул сложилась из нескольких главных ветвей. Вот что писал по этому поводу в 1972 г. академик С.Л. Соболев [261]:
“Работы по теории кубатурных формул относятся к различным направлениям. В одном из них, в продолжение классических исследований Радона, отыскиваются формулы, которые могли бы точно интегрировать возможно ббльшее число из данного набора функций, например, многочленов или тригонометрических функций. ...Второе направление — направление вероятностных методов Монте-Карло. Здесь проблема состоит в том, чтобы построить алгоритмы не вполне случайные, основанные на некотором регулярном выборе случайных точек.
Наконец, третье направление касается изучения формул, использующих
Введение
6
в качестве узлов точки правильной решетки. Здесь удается сравнительно далеко продвинуться в изучении нормы функционала погрешности в различных функциональных пространствах, и, в частности, для пространств с ортогонально инвариантной нормой найти зависимость погрешности от решетки.”
В обзорной статье [338, с. 481] всевозможные кубатурные формулы подразделяются на четыре больших класса: к первому из них автор относит “эффективные формулы”(в частности, декартовы произведения одномерных формул и формулы, имеющие при заданной полиномиальной точности наименьшее число узлов). Второй класс в его же терминологии образуют “формулы с минимальной нормой”; к третьему классу он относит “теоретико-числовые формулы”, а к четвертому — “вероятностные формулы”.
По поводу теоретико-числовых методов интегрирования Н.С. Бахвалов в своих комментариях [22] к разделу II книги [295] пишет:
“В работе [90] рассмотрен вопрос о повышении эффективности метода Монте-Карло при вычислении континуальных интегралов. Существенным моментом, обеспечивающим это повышение, было впоследствии ставшее классическим следующее преобразование исходного континуального интеграла. Вместо общепринятого в то время разыгрывания движения частиц, при котором разыгрывается поглощение и вылет, было предложено разыгрывать движение частиц с переменной массой, при котором поглощение и вылет учитываются при каждом столкновении изменением этой массы. Такое преобразование приводит к существенному сглаживанию подынтегральной функции и соответственно к уменьшению дисперсии.
Результатом исследования этой проблемы явилось появление так называемых теоретико-числовых методов интегрирования.”
Теоретико-числовым методам построения кубатурных формул посвящены многие работы Н.М. Коробова, а также его монография [133].
Вероятностным формулам посвящено, наверное, наибольшее в рамках рассматриваемой теории число работ. Познакомиться с вероятностным подходом в теории кубатурных формул можно, например, по работам [282], [104], [88], [89], [103], [295], [80], и это, конечно же, совсем не полный список. Но в настоящей диссертации вероятностные методы Монте-Карло не затронуты совершенно, основное же внимание автора сосредоточил на фор-
Введение
7
мулах, хорошие аппроксимационные свойства которых обосновываются методами функционального анализа и классической теории приближений.
3°. Функциональные методы стали широко применяться в теории приближенного интегрирования начиная с работ академика С.М. Никольского и первого издания его книги “Квадратурные формулы” [160]. Изложенные в ней результаты относятся, в основном, к квадратурным формулам на классах WmLp функций одной переменной. Создание же теории кубатур-ных формул заслуженно связывают с исследованиями академика С.Л. Соболева. В научном наследии С.Л. Соболева работы по теории приближенного интегрирования, выполненные им в “сибирский” период жизни, занимают весьма заметное место: первую работу по кубатурным формулам он опубликовал в 1961 г., последнюю — в 1986 г., всего же их более трех десятков, в том числе две фундаментальные монографии [263] и [276]. Вторая из упомянутых монографий в 1997 году переведена на английский язык [281].
Объектом исследования в теории кубатурных формул служат последовательности приближенных равенств, предназначенных для аппроксимации и вычисления многомерных интегралов по ограниченной области О, пространства Rn, п > 2. При этом предполагается, что границей области интегрирования Q служит кусочно гладкая поверхность конечной площади, а в остальном Q произвольна. Как уже отмечалось, практика современных теории приближений и численных методов такова, что для приближенного вычисления интеграла по области Q чаще всего используются кубатурные формулы, т. е. приближенные равенства вида
Г N
/ ip(x) dx~y^ ск<р(х{к)) = 0. n *=I
Точки х^к\ х^ € Kn, и параметры с* называют соответственно узлами и весами (коэффициентами) кубатурной формулы. На практике узловые точки хW внутри О, часто заданы внешними обстоятельствами или же используется множество “равноудаленных” узлов [293, с. 95].
Во втором случае мы имеем дело с результатами, относящимися к формулам с параллелепипедалъпой решеткой узлов. В этом случае векторы хW нумеруются с помощью мультииндексов ß = yßn) Е Zn с цело-
численными координатами, т.е. любой из узлов х^ можно найти по фор-
Введение
8
муле = хр = /гЯ/З, где Я — матрица размера п х п, с1еЬ Я = 1, а положительный параметр Л называется шагом решетки.
Функциональный подход к построению и исследованию формул для приближения многомерных интегралов предполагает, во-первых, использование выбранной или построенной формулы не столько для какой-то одной конкретной функции, сколько сразу для целого их семейства, представляющего собой шар в некотором наперед заданном банаховом пространстве X.
Во-вторых, разность между интегралом и приближающей его конечной комбинацией значений подынтегрального элемента рассматривается как результат действия на подынтегральную функцию некоторой обобщенной функции, полностью определяемой исходной кубатурной формулой и называемой по этой причине ее функционалом погрешности.
В-третьих, исходное банахово пространство X предполагается вложенным непрерывно в пространство С(Г2) функций, непрерывных в замыкании области интегрирования П. В этом случае функционал погрешности кубатурной формулы не только линеен, но и ограничен на X, причем знание численной мажоранты для его нормы в сопряженном пространстве X* позволяет получать для произвольной функции из единичной сферы пространства X гарантированные оценки близости истинного значения интеграла от этой функции к рассматриваемой на ней кубатурной сумме.
В рамках функционального подхода классической теории кубатурных формул традиционно ставится ряд задач общего характера, решение которых в частных случаях, собственно, и обеспечивает дальнейшее развитие и совершенствование метода кубатурных формул. В качестве иллюстрации приведем здесь формулировки некоторых задач.
1. Задавшись областью интегрирования, описать конструкцию каждой из последовательных кубатурных формул, т. е. указать для каждой из них множество узлов и весов, либо предложить эффективный алгоритм их нахождения. Указать полиномиальную (или тригонометрическую) степень точности рассматриваемой формулы.
2. Задавшись “пробным” банаховым пространством Х} найти асимптотическое разложение нормы функционала погрешности в X* по числу узлов N соответствующей кубатурной формулы — важнейшему параметру вычислительного процесса. Вывести эффективные оценки для упомяну-
Введение
9
той нормы как сверху, так и снизу в виде явных функций от N.
3. В качестве исходного банахова пространства X на практике выбираются самые разнообразные классы функций. Разные авторы рассматривали кубатурные формулы в пространствах Соболева (изотропных и неизотропных, конечного и бесконечного порядков), Бесова, Никольского, в классах Гельдера, Жевре, пространствах с доминирующей смешанной производной, весовых функциональных пространствах, в бесселевых шкалах пространств и др. Такое изобилие в выборе “пробных” банаховых пространств явно свидетельствует о необходимости доказательства утверждений общего характера (теорем о сходимости, об оптимальных формулах, об асимптотически оптимальных формулах и т.п.), верных одновременно для всех перечисленных выше классов, либо по крайней мере для большинства из них.
4. Рассматривая совокупность индикаторов ограниченных областей в Кп как достаточно представительное подмножество пространства финитных обобщенных функций, естественно заключить, что теорию кубатуриых формул помимо прочего возможно рассматривать как часть общей теории обобщенных функций. В этой связи стоит задача распространения результатов классического метода кубатурных формул на случай приближения произвольной финитной обобщенной функции конечными линейными комбинациями сдвигов дельта функций Дирака.
Таким образом, функциональный подход помимо описания конструкций рассматриваемых формул, т. е. указания их узлов и весов либо алгоритмов их нахождения, а также указания в ряде случаев полиномиальной степени рассматриваемой формулы, подразумевает и вывод эффективных двусторонних оценок для норм соответствующих функционалов погрешности. В случае параллелей и педальных решеток узлов асимптотические разложения этих норм по шагу решетки интегрирования найдены для широкого круга банаховых пространств.
Функциональный подход порождает также критерий качества кубатур-ной формулы отличный от критерия, порождаемого подходом алгебраическим: в первом случае предпочтительней та формула, функционал погрешности которой имеет меньшую норму, во втором лучшей считается формула, точная на возможно большем числе полиномов.
Введение
10
Важную часть рассматриваемой теории составляют результаты по ку-батурным формулам, обладающим высокой полиномиальной степенью и инвариантным относительно преобразований группы вращений некоторого правильного многогранника [244], [276, гл. И], [152], [233]—[241], [137]-[143]. Требование точности кубатурной формулы с заданными узлами на многочленах до определенной степени сводит задачу отыскания ее весов к решению линейной системы уравнений. Чем выше требуемая точность и чем больше узлов, тем большие размеры имеет эта система. Однако в случае, когда область интегрирования обладает определенного рода симметрией, а для приближенного интегрирования используется инвариантная кубатурная формула, размеры соответствующей линейной системы можно существенно уменьшить (см. [276, гл. II]; там же предложен алгоритм построения узлов инвариантной кубатурной формулы на трехмерной сфере).
Наиболее развитое направление теории состоит из результатов по асимптотически оптимальным решетчатым кубатурным формулам в пространствах функций конечной гладкости [276, гл. IV—VI].
С.Л. Соболев рассматривал в этой связи гильбертовы пространства Ь Предложенная им конструкция регулярного пограничного слоя позволяет при сколь угодно большом числе узлов находить веса кубатурной формулы, решая лишь несколько стандартных систем линейных уравнений с размерами, зависящими только от гладкости га. Центральный результат исследований по асимптотически оптимальным формулам — это вывод асимптотического представления Х^^-нормы функционала погрешности с регулярным пограничным слоем. Соответствующая формула включает в себя два слагаемых, первое из которых записано в явном виде через так называемые обобщенные числа Бернулли, а второе пренебрежимо мало по сравнению с первым при малом шаге Н решетки интегрирования. В частности, из этого представления следует, что норма функционала погрешности с регулярным пограничным слоем убывает при к —> 0 как степенная функция кт. Это весьма глубокий аналитический факт, позволяющий дать функциональное определение порядка точности кубатурной формулы на классе [274].
Формула, выражающая Ь^^-норму функционала погрешности с регулярным пограничным слоем через обобщенные числа Бернулли, дает серьезное основание в пользу выбора в качестве множества узлов интегри-
Введение
11
рования точек параллелеп и педальной решетки. В самом деле, мы можем поставить задачу отыскания для заданного числа узлов N наилучшей кубатурной формулы, т.е. такой кубатурной формулы, чей функционал погрешности имеет приданном N наименьшую Х/^'-иорму, причем минимум берется не только по весам, но и по узлам формулы. При этом отношение Х^'-нормы наилучшей кубатурной формулы к Х^^-норме функционала погрешности с регулярным пограничным слоем, имеющего то же самое число узлов, ограничено снизу положительной величиной от N не зависящей. Это заключение сразу следует из теоремы Бахвалова, формулировку и доказательство которой можно найти, например, в [276, глава IV, § 3]. Тем самым вряд ли при увеличении N можно получить большую выгоду от использования вместо формул с узлами в точках параллелеп и педальной решетки каких-либо других формул с произвольным распределением узлов, тем более, что оптимизация формулы по узлам связана с решением систем нелинейных уравнений высокого порядка, что само по себе является задачей весьма трудоемкой. В противоположность этому узлы кубатурных формул с регулярным пограничным слоем заданы явно.
Полезно отметить, что всякая формула с регулярным пограничным слоем, представляя собой многомерный вариант решения классической задачи суммирования функций, относящейся к исчислению конечных разностей [91], является, по существу, многомерным аналогом классической квадратурной формулы Грегори. Тем самым поведение подынтегральной функции вблизи границы области интегрирования при построении формулы учитывается посредством специального задания тех весов этой формулы, которые соответствуют узлам решетки, лежащим в некоторой приграничной полосе. Во всех остальных узлах веса формулы с регулярным пограничным слоем одинаковы.
Замечателен метод, предложенный С.Л. Соболевым для явного представления Х^^-нормы функционала погрешности 1(х) и использующий понятие экстремальной функции и(х). Эта функция рассматривается как обобщенное решение многомерного полигармонического уравнения со специальной правой частью:
Ати(х) = (—1)т/(х).
Решение этого уравнения на числовой прямой — это кусочно полиноми-
Введение
12
альная функция класса L^, т. е. сплайн. В многомерном же случае такой подход позволил привлечь к исследованию классической проблемы анализа прекрасно развитые методы решения дифференциальных уравнений с частными производными.
В случае, когда область интегрирования представляет собой рациональный многогранник, построение формул с регулярным пограничным слоем сопряжено с конструкцией формального пограничного слоя [276, гл. III], позволяющей создавать относительно простые алгоритмы для вычисления весов решетчатых кубатурных формул.
Важные результаты по асимптотически оптимальным формулам получил М.Д. Рамазанов [204]—[214]. В.И. Половинкин [165]—[194] обобщил теорию асимптотически оптимальных формул на пространства £рт^(Г2) при 1 < р < оо и на случай весового интегрирования.
О.В. Бесов вывел асимптотические формулы для норм дроблений данного функционала погрешности [29]—[32]. Каждое такое дробление — это множество функционалов погрешности, действующих в более общих, чем пространствах. К ним, в частности, относятся классы функций, обобщенные производные порядка т которых принадлежат пространству Марциикевича М^(П) (либо пространству Лоренца А^(П), либо пространству Орлича ь\,т-
Кубатурные формулы в пространствах Wp7n* и в неизотропных пространствах Соболева рассматривались Ц.Б. Шойнжуровым [300], [301].
Кубатурные формулы с регулярным и формальным пограничным слоем исследовали Н.И. Блинов и Л.В. Войтишек [33]—[37], [82]-[8б], Ф.Я. Загирова [112].
Асимптотическое поведение норм функционалов погрешности в случае анизотропных классов функций Wp^ М исследовали М.Д. Рамазанов [212] и В.Н. Темляков [285]. В.Н. Темляков [285] получил также асимптотические разложения норм функционалов погрешности для анизотропных пространств Никольского NHq([0, 2л-]п).
Особое направление теории составляют исследования по сходимости кубатурных формул на пространствах бесконечно дифференцируемых функций [276, гл. VII]. В качестве исходных классов с указанным свойством бесконечной гладкости рассматриваются, например, пространства перио-
Введение
13
дических функций многих переменных с заданным поведением интегральных £г2 -норм при неограниченном возрастании т [260], [274). Используемая в [276, гл. VII] классификация охватывает известные пространства целых функций, имеющих определенный тип и порядок роста, пространства аналитических функций, а также классы Жеврс, включающие в себя квазиаиалитические элементы.
Рассмотрев действие в периодических пространствах указанного типа решетчатого функционала погрешности с равными весами, С.Л. Соболев получил асимптотическое представление логарифма нормы этого функционала. Как и в случае пространств конечной гладкости, соответствующая формула состоит из двух слагаемых. Первое из них явно выражено через параметры исходного класса, а второе пренебрежимо мало по сравнению с первым при малом шаге Н решетки. Эти результаты обобщили работу П. Дэвиса [323], в которой рассматривались кубатуриые формулы в пространствах аналитических функций одной переменной с дополнительным условием периодичности.
В одномерном случае ряд авторов в качестве исходного класса элементов бесконечной гладкости рассматривал пространство функций, аналитических в окрестности единичного отрезка (точнее, в эллипсе с фокусами в точках ±1). В частности, для этого пространства В.И. Крылов доказал сходимость интерполяционного квадратурного процесса [134], а Н.С. Бахвалов оценил скорость сходимости квадратурного процесса Гаусса [17]. Асимптотические разложения норм функционалов погрешности квадратурных формул Грегори и Эйлера — Маклорена в пространствах бесконечнодифференцируемых функций, являющихся одномерными реализациями соответствующих пространств Соболева из [274], найдены в [43]—[48]. Как оказалось, в общем случае использование при данном Н формул высокой полиномиальной степени не только не приводит к уменьшению нормы функционала погрешности, но, напротив, может повлечь за собой прямо противоположные последствия.
Обзор свойств формул приближенного интегрирования в пространствах бесконечно дифференцируемых функций сделал К.И. Бабенко [3].
Квадратурные формулы в пространствах бесконечно дифференцируемых функций рассматривал также В.Н. Белых, предложивший алгоритмы без насыщения в задаче численного интегрирования [25], исследовавший
Введение
14
ненасыщаемые квадратурные формулы на отрезке [26], а также указавший пример приложения этих квадратур к численному решению уравнения Лапласа [28].
Переход от функций конечной гладкости к бесконечно дифференцируемым сопровождается принципиально важным для практики эффектом: норма функционала погрешности кубатурной формулы, в первом случае убывающая не быстрее некоторой степени шага решетки интегрирования, во втором случае убывает быстрее любой такой степени. В этой связи было предложено считать, что кубатурная формула имеет в данном банаховом пространстве бесконечный порядок, если норма соответствую и щ го функционала погрешности в сопряженном пространстве убывает к нулю быстрее любой степени шага сетки интегрирования [274]. Пример кубатурной формулы бесконечного порядка дает известная теорема о среднем для гармонических функций. Эрмитовы кубатурные формулы бесконечного порядка в пространствах гармонических и полигармонических функций рассмотрены соответственно в работах [69] и [325]. До сих пор кубатурных формул, обладающих бесконечным порядком, известно весьма немного.
Исследования по решетчатым кубатурным формулам в пространстве с оптимальными весами составляют еще один важнейший раздел теории (276, гл. IX]. Центральное место в полученных в этом направлении результатах занимает предложенный С.Л. Соболевым аналитический алгоритм отыскания весов оптимальных кубатур. Соответствующая формулировка использует специальный разностный оператор [254], действие которого на функцию дискретного аргумента [276, гл. VIII] — это свертка, аналогичная действию полигармонического оператора Дт на непрерывно дифференцируемую функцию.
Задача вычисления сверточного ядра при произвольном т оказалась весьма непростой. Частично ее удалось решить в одномерном случае: здесь имеется формула, выражающая искомые значения через корни многочлена Эйлера степени 2га. Веса оптимальных решетчатых кубатурных формул представимы значениями в точках Zn финитной функции дискретной переменной, удовлетворяющей линейному конечно-разностному уравнению со специальной правой частью. Действие на эту правую часть сверточного дискретного аналога полигармонического оператора приводит к аналитической формуле для искомых весов. Результаты по весам оптимальных квад-
Введение
15
ратурных формул, изложенные в (276, гл. IX, §§ 6,8], обобщили некоторые результаты А. Сарда [360], И. Мейерса, И. Шенберги, С. Силлимена [359]— [365]. Ряд вопросов, возникающих в процессе реализации предложенного С.Л. Соболевым алгоритма отыскания весов оптимальных кубатурных (и квадратурных) формул, решен З.Ж. Жамаловым [ 107]—[ 108], Ф.Я. Загировой [113]—[114], М.Д. Рамазановым и Х.М. Шадиметовым [221], [297], [296].
4°. Приведенный в предыдущих пунктах очень краткий и не претендующий на полноту тематический обзор включил в себя лишь теоретические результаты, относящиеся к анализу ошибок ограничения метода кубатурных формул, т. е. к анализу погрешности, выполненному в предположении, что численный алгоритм реализуется в обладающей бесконечной точностью арифметике. Результаты же о распространении ошибок в начальных данных алгоритма (т. е. ошибок в задании весов формулы и значений подынтегральной функции в узлах) и ошибок округления, в литературе представлены весьма скупо и фрагментарно. В этой связи можем упомянуть здесь лишь книгу Д. Мак-Кракена и У. Дорна [147], в главе 6 которой приведена оценка сверху общей суммарной погрешности квадратурных формул трапеций и Симпсона.
Указанный пробел в анализе суммарной ошибки метода кубатурных формул восполняет в какой-то степени настоящая диссертация. В ней приводятся относящиеся к проблематике метода кубатурных формул результаты исследований, начатых автором в 1978 году под научным руководством академика С.Л. Соболева. Вопросы приближенного многомерного интегрирования рассматриваются здесь для случая многомерных ограниченных областей интегрирования и классов непрерывных в этих областях подынтегральных функций.
5°. Диссертация носит теоретический характер, состоит из введения, основной части, включающей три главы, списка литературы и оглавления. Краткое содержание основной части приводится далее.
В главе I диссертации определяется ряд основных понятий, необходимых для организации вычисления многомерных интегралов с гарантированной точностью.
В § 1 главы I приводятся первоначальные и самые элементарные сведения теории: определяются кубатурная сумма, веса и узлы кубатурной
Введение
16
формулы, а также ее функционал погрешности
N
Л <w(x(<:)).
к=1
В §2 главы I приводятся простейшие сведения о специальных конечных числовых множествах, называемых форматами, элементы которых — суть машинные числа. Элементы теории машинных вычислений излагаются здесь в рамках общепринятой модели арифметики с плавающей точкой, разработанной во второй половине прошлого века.
Исторические и библиографические сведения об арифметике с плавающей и фиксированной точкой имеются в книге Д. Кнута (125, с. 239, 257). Как отмечается в [93, с. 443] со ссылкой на свидетельство Дж. Уилкинсона [287, с. 178], первый детальный анализ ошибок округления при арифметических операциях с фиксированной точкой был дан Дж. фон Нейманом и Голдстайном в 1947 г., а для операций с плавающей точкой аналогичный анализ дал Дж. Уилкинсон в 1957 г. Модель арифметики с плавающей точкой широко применялась, получая при этом дальнейшее развитие, в работах С.К. Годунова и его учеников (см., например, [93, глава 4]), в книгах Дж. Форсайта, М. Малькольма и К. Моулера [291] и Д. Каханера, К. Моулера, С. Нэша [122, глава 2], а также в книге Дж. Голуба и Ч. Ван Лоуна [94]. В настоящее время понятие машинных чисел активно используется в современных образовательных курсах. Из недавно вышедших следует отметить книгу Н.С. Бахвалова, A.B. Лапина и Е.В. Чижонкова [23] и учебник по анализу [305], а из изданных относительно давно — упомянутую книгу Д. Мак-Кракена и У. Дорна [147].
Для каждой конкретной машины формат определяется четырьмя целочисленными параметрами. Пусть 7 — целое число, 7 > 1, Р+ и к натуральные числа, а Р- отрицательное целое число. Форматом F = F(y, Р-} Р+, к)
Г N
~ / [хп{х) -£**(* “ xw)](p(x)dx
J k= 1
= J ip(x) dx —
fi
Введение
17
называется числовое множество, определяемое соотношением
к
г = Р('у,Р.,Р+,к) = {0} и {г е Л I 2 = ±7р^т37--’, Р- < р < Р+,
;=1
р — целое; т1,т2,..., га* — целые числа из [0,7 — 1]; тп\ > 0}.
Параметр 7 при этом называется основанием формата, к — его точностью, а Р- и Р+ — границами области значений [Р_,Р+] показателя р формата. Помимо числового вектора (7, Р_, Р+, &) формат Г = Г(7, Р_, Р+, к) принято описывать эквивалентным образом, с помощью вектора (7»е0>^ь^оо)> где бо, 61 и боо — так называемые форматные константы. Их связь с параметрами 7, Р_, Р+ и к задается следующим образом:
б0 = 7Р--1. е1=71_\ «оо = 7Р+(1 ~7-*)-
Округление вещественного числа г из отрезка [—Соо^), боо(Р)] представляется как результат действия на 2: специального оператора
А : [_ево(Е),€во(Г)] -> Г,
сопоставляющего числу 2: соответствующий ему в формате Р элемент [г]?
наилучшего приближения. Для заданного вектора д = (д\,д2, ••• ><7лг) при
шах |$| < бГоо(Р) определен вектор 1 <i<N
Я(р) = {ЯЫ.йЫ Й(Ы} = Ыг-
Погрешность, возникающая при замене вектора д вектором А(д) допускает следующую оценку
1|р - ЙЫ!Ь < £1(Г)|Ы|2 + %/Ме0(Г).
Следует учитывать, что на конкретной вычислительной машине помимо аппаратно реализованного в ней “максимально” доступного формата возможно использовать и любой другой формат с тем же основанием, но с меньшим числом элементов. Это позволяет при реализации конкретного вычислительного алгоритма производить отбор самого виртуального формата, в котором этот алгоритм и будет реализовываться. Виртуальный
Введение
18
формат вовсе не обязан совпадать с максимально доступным, а удачный его выбор может существенно упростить вычисления. Естественно, что при этом следует заранее убедиться, что требуемая точность результата в выбранном виртуальном формате будет достижима.
В §4 главы I проведен анализ погрешности кубатурной формулы Едг с учетом как ошибок ввода так и сделанных в процессе реализации вычислительного алгоритма округлений.
Задав для произвольного функционала погрешности Iдг его реализацию в формате Е как следующую обобщенную функцию
Мр(*) = хп(я) - ^[с*]г<5(х - х^),
*= 1
определим затем (^(^»^(Е^-окрестность функционала /дг. Именно, будем говорить, что функционал погрешности вида
Лг
(я) = хп(я) - - х{к))
*=1
принадлежит (ео(Р), £1(Е))-окрестности функционала /^, если имеет место равенство
Таким образом, при заданном стандартном формате Е с константами €о(¥) и е\(Е) все функционалы Iдг из (ео(Е),£л(Е))-окрестности данного функционала /дг, включая и его самого, вводятся в ЭВМ одинаково.
Следовательно, для любой непрерывной в замыкании области Г2 функции <р(х) результат вычисления на ЭВМ кубатурной суммы один и тот же сразу для всех функционалов /дг из упомянутой (ео(Е),б1(Е))-окрсстности функционала /дг, и этот результат мы обозначаем через Едг(у?). С учетом этого обозначения рассмотрим еще один, уже нелинейный, функционал
= |/п(*>) - Едг(у?)|.
Величина /2р(у>, Ед/) допускает простую геометрическую интерпретацию. Именно, отрезок 1° ЧИСЛОВОЙ ОСИ С центром С В точке Едг (у?) и длины /(V?, Едг) = 2Лр(у?, Е/у) обладает, очевидно, следующими свойствами: во-первых, числовое значение интеграла от у? принадлежит Iе, и, во-вторых, Iе наименьший среди всех отрезков с центром в с = Едг (у?), обладающих
Введение
19
первым свойством. По этой причине Iе = Ig(<Pj E;v) назван гарантированным отрезком.
Чтобы оценить длину l(ip, Едг) гарантированного отрезка или, что то же самое, величину Rp(<P, %n), сделаем некоторые дополнительные предположения.
Пусть функция <р — это элемент некоторого банахова пространства X, непрерывно вложенного в пространство С(П). Константой этого вложения назовем такое конечное положительное число Л, что
sup|y>(i)| < А\\<р \ Х\\.
Определение. Половину точной верхней грани длин f(<p, Е/\г), взятой по всевозможным функциям (р из единичной сферы пространства X, т. с. величину
RF = RF(Xy Е/у) = sup |/ц(у?) - E/v(y?)|,
назовем гарантированным радиусом кубатурной формулы Е/у.
Как следует из этого определения, все кубатурные формулы с функционалами погрешности из (eo(F),ei(F))-OKpecTiiocTH данного функционала /дг имеют одинаковый гарантированный радиус.
Для любой функции единичной сферы пространства X интеграл 1п(<р) всегда принадлежит отрезку [с — <т, с + <т], где с = E/vCv7) — машинное число, полученное в результате реализации выбранного алгоритма вычисления кубатурной суммы, а с = Rf(X, Ejv) — гарантированный радиус кубатурной формулы. Из двух кубатурных формул лучшей на пространстве X является та, гарантированный радиус которой меньше, а среди кубатурных формул с одинаковым гарантированным радиусом более экономичны те, которые имеют минимальное число узлов.
Как показывают сделанные относительно Rf(X, E/v) заключения, именно его естественно рассматривать в качестве параметра, гарантирующего точность вычисления интеграла с помощью выбранной кубатурной формулы. Радиус Rf(X, Е/у), как это вытекает из его определения, является помимо прочего функцией алгоритма, посредством которого введенные в машину данные перерабатываются в число Е/у(<р). Оценка радиуса Rf(X, Е/у)
Введение
20
проведена в § 4 главы I для случая, когда в качестве алгоритма вычисления скалярного произведения в арифметике с конечной точностью выбран простейший алгоритм 1 (подробнее см. § 4 главы I). Этот алгоритм применяется к машинным векторам (с]р и [<^]р» соответствующим вектору с весов кубатурной формулы и вектору ф значений подынтегральной функции, а получающееся в результате машинное число, как уже было оговорено, обозначается через Ел*(<р). Доказана следующая
Теорема 2. При заданных форматных константах ео(Г), сі(Р), ^(Г) и константе вложения А гарантированный радиус Яр(Х}^2^) кубатурной формулы Е/у с вектором с весов таким, что ||с|І2 < \Доо(^)/2, допускает следующую оценку
|яР(ЛГ, Елг) - ||/* | хщ < 4ТУ ( |с*| + €0(Р)
г<?е N < тіп 1^1) .
Если кубатурная формула Е/у имеет неотрицательные веса и при этом точна па любой тождественно постоянной функции, то для ее гарантированного радиуса при указанных выше значениях ТУ справедлива оценка
Лр(Х, Елг) - IIїм | Х*||| < 4ЛГ (у1сі(Р)|П| + ео(Р)).
где через |П| обозначен объем области ЄІ.
/ м \ 1/р
Пусть 1 < р < оо и гр(Х,1Ч,£1) = (||г | Х*\\р Н- (4ТУЛєі)р £ Мр) •
Тогда в условиях теоремы 2 справедливо неравенство
йр(ЛГ>Елг)<(^+1)1/%(ЛГ,Лг,Єі(Р))+4ЛГЄо(Г), где - + - = 1.
Я Р
Из теоремы 2 следует также, что при данном N и одновременном стремлении к нулю форматных констант ео(^) и є\(¥) радиус Яр(Х, Е/у) имеет пределом норму ||//у | Х*\\ функционала погрешности кубатурной формулы. Очень интересен вопрос о поведении Пр{Х, Е/у) при ТУ —» оо и фиксированных константах Єо(Р), £Гі(Г)), но ответ на него не представляется очень простым.
Введение
21
Гарантированный радиус Яр(Х,Е#) при вычислениях в арифметике с конечной точностью играет ту же роль, что и норма функционала погрешности кубатурной формулы в классических задачах теории кубатурных формул.
Вместе с тем при заданных формате F и функциональном пространстве X семейство всевозможных гарантированных радиусов Ry(X,Y^^) ограничено снизу некоторой строго положительной величиной. Это означает, что в арифметике с конечной точностью неограниченное увеличение числа узлов кубатурной формулы в принципе не приводит к неограниченному же уменьшению гарантированного радиуса.
В целом качество всей совокупности имеющих заданное множество узлов Д кубатурных формул предлагается характеризовать следующим набором числовых характеристик
Рв{ F,X,A)= inf Др(ЛГ,Елг).
C=(Ci,C2,—,C/v)
Здесь F — стандартный формат, X — банахово пространство, а инфимум берется по всевозможным весам с = (ci,C2,... , Cjv) кубатурных формул. Посредством величин Рв(F,X, Д) в пространстве непрерывных функций <7(Q) выделяется класс “практически интегрируемых” функций.
Определение. Функцию ip 6 C(Q) назовем “практически -интегрируемой” с точностью е > 0 в стандартном формате F, если найдутся такие банахово пространство X и множество узлов Д, что
tp е X и Рв(FyX,A)<e.
Как уже отмечалось, со временем аппаратные средства все более и более совершенствуются и точность машинной арифметики возрастает. Это означает, что для форматных констант eo(F) и 61(F), становятся допустимы все меньшие, хотя и строго положительные, значения. Но при одновременном и по возможности согласованном стремлении параметров co(F) и 61(F) к нулю, а параметра 6oo(F) к бесконечности, гарантированный радиус R?(X, Еуу) стремится к норме в X* соответствующего функционала погрешности. Это замечание устанавливает естественную преемственность задач приближенного интегрирования в арифметике с конечной точностью с задачами классической теории кубатурных формул: последние являются,
Введение
22
по существу, предельными для задач, возникающих при анализе погрешности с учетом операций арифметических округлений.
В главе II диссертации исследуются инвариантные кубатурные формулы типа Грегори для простейшей области пространства Rn — его единичного куба
Q = {х £ Rn | х = (xltx2i... ,xn), 0 < xj < 1, j = l,2,...,n} .
Соответствующие кубатурные формулы имеют при этом параллелепипе-дальную решетку узлов которые нумеруются с помощью мультииндск-са р = (/?i,... ,0п) с целочисленными координатами и задаются формулой х/3 = hp, где hp £ Q, a h — положительный параметр, называемый шагом решетки и выбираемый таким образом, чтобы 1//*было натуральным. Тем самым компоненты узлов инвариантных кубатурных формул главы II явно заданы как рациональные числа.
Построения главы II начинаются с рассматрения в §1 многомерных аналогов квадратурных формул прямоугольников и трапеций. Для заданного шага h решетки узлов определены следующие обобщенные функции
Pi,h(z) = hn^2 5{х - /17), Т„(х) = - cq) + CQ).
hyeQ Ле£>(п)
Здесь и всюду далее £(х) — дельта функция Дирака, 7 = (7ь«-*»7п) ” мультииндекс с неотрицательными компонентами, сq — центр куба Q, D(n) — группа, образуемая диагональными матрицами размеров п х тг с ±1 на главной диагонали. Равенство
задает для куба <3 аналог функционала погрешности хорошо известной в случае отрезка квадратурной формулы прямоугольников, аналог же функционала погрешности хорошо известной квадратурной формулы трапеций задается равенством