О КРИТЕРИЯХ ПРИ ОЦЕНКЕ ОСТАТКА КУБАТУРНЫХ ФОРМУЛ*
Т. М. Товстик
С.-Петербургский государственный университет, канд. физ.-мат. наук, доцент, [email protected]
В статье дается алгоритм вычисления критериев, участвующих в оценке остатка кубатурных формул при вычислении п-мерных интегралов от функций, имеющих интегрируемые смешанные производные не выше порядка 2п, причем не выше второго порядка по каждой из переменных. Расматриваются функции в двухмерном и п-мерном единичном кубе. Критерии представляют собой константы, участвующие в оценках остатков кубатурных формул, и равные максимумам вспомогательных функций, зависящих от параметров кубатурной формулы. Частным случаем критериев является дискрепанс.
Введение. При вычислении многомерных интегралов часто используются [1, 2] кубатурные формулы
где п — размерность области интегрирования, а Ск и X (к) —параметры кубатурной формулы.
В статье [3] анализируются остатки кубатурных формул
в зависимости от свойств интегрируемых функций, заданных в п-мерном единичном кубе К” = [0,1]”. В частности, рассматриваются функции /(X), имеющие интегрируемые вторые производные по каждой из переменных. Введение критериев позволяет выразить остатки через интегралы, содержащие вторые производные функций / (X). Ниже этот вариант рассматривается более подробно.
Рассмотрим определение критериев в случае п = 2. Использование разложения в ряд Тейлора функции /(X) в точке х\ = Х2 = 1 после ряда преобразований дает возможность записать остаток кубатурной формулы в виде
N
/(X) ^ = ^3 с/^к)), X = (Х1,...,Х”), X(k) = (Хl(k),...,Хn(k)), (1)
к=1
(2)
* Работа выполнена при финансовой поддержке РФФИ (грант №08-01-00194). © Т.М.Товстик, 2009
П1
/(2Х 2 (иьи2)
.2^.2
N
Ь 1 ио X
12 Е
к=1
4
СкЬ(м1 - Х1 (^))Ь(м2 - *2(^))
содержащем интегралы от соответствующих производных [3]. Здесь
т/ \ | г, 2 > 0,
) По, 2 < 0,
а коэффициенты Ск и узлы X(&) удовлетворяют соотношениям
ЙМ1 ^М2, (3)
(4)
N
Е
к=1
N
N
N
Ск
1, ^СкХ1 (й) = ^ Ск Ж2(&) = 1/2, ^Ск Ж1(*)Ж2 (&) = 1/4,
(5)
к=1
к=1
к=1
которые имеют место, если кубатурная формула точна для многочлена вида /(х1, Х2) = ао + а1Ж1 + а2Ж2 + 03X1X2. Супремумы модулей выражений, стоящих в квадратных скобках формулы (3), называем критериями [3]. Оценка остатка Д[/] выражается в виде суммы произведений критериев на интегралы от модулей производных функции /(X), входящих в (3).
1. Двухмерный случай. Из формулы (3) следует, что принципиально различными являются три критерия:
С(1) = яир |у>(и)|, 0<«<1
С(2)= яир |^(и) 0<«<1
С(3) = яир |0(и, V)
0<и,^<1
N
¥(и) = у ~^2скЬ{и - х (к)), к=1
1/2 N N
Ф(и) = 2 ( У ~2^2скЦи-х(к))(1 — у (к))
(6)
к=1
N
(и, V)
- СкЬ(и - х(*))Ь(^ - у(&)),
к=1
где X = Х1, у = Х2, и = и1, V = и2.
Замечание 1. Критерий 6(2) является частным случаем критерия 6(3) при V = 1, поэтому 6(3) > 6(2). При выполнении условий (5) справедливы равенства
N N 1 N 1
Есй(1 -х(к)) = Ес*(1 -у(к)) = Ес^1-^))(1-^)) = 1’
к=1
к=1
к=1
поэтому, если ввести обозначения Ск = 2Ск(1 - у(^)), то ^=1 Ск = 1, и удвоенный
критерий 26(2) будет иметь вид 6(1). Если при этом Ск = Ск, то получаем 6(2) = 1/26(1).
Замечание 2. Как и в [3], разложение в ряд Тейлора можно также произвести в точке х = у = 0. Пусть С?(3) —аналог критерия С(3) в этом случае:
6(3) = ЯИр |0(и, V)!, С(и, V)
0<и,^<1
(1 — и)2{ 1 — уУ 4
N
Е
к=1
Ск£(х(&) - и)Ь(у(^) - V) (8)
Вычисление критерия С?(3) сводится к вычислению критерия С(3) в результате замен ми V на 1 -и и 1 - V и х(^) = 1 - х(&), у(^) = 1 - у(&), 1 < & < N. Если выполнены
22
и^2
4
соотношения (5), то аналогичные соотношения выполняются и для x(k), y(k) (см. (Т)). Критерии С(З) и С?(З) совпадают, т. е. С(З) = С?(З), если при ранжировании в порядке возрастания узлов x(k) и узлов y(k) при одинаковых x, для всех k выполняются условия симметрии
x(k) = 1 - x(N - k + 1), y(k) = 1 - y(N - k + 1), ck = cn-k+l, 1 < k < N. (9)
2. Вычисление критерия G(l). Здесь и далее считаем, что точки x(k) различны
и расположены в порядке возрастания: О = x(l) < x(2) < ••• < x(N) = 1. Если точек x(l) = О и/или x(N) = 1 не было среди исходных данных, то добавим их с весом сі = О и/или cn = О, причем величина G(l) от cn не зависит. Так как функция ^(u) не меняет своего вида в областях
A, = [x(i), x(i + 1)), 1 < i < N - 1, (1О)
будем искать sup |^(u)| последовательно в этих областях, начиная с i =1. Положим
u2
= sup \ifi{u)\, ipi(u) = ip(u) = — ~22ck(u - x(k)) при и Є Д*. (11)
k=1
Тогда
где
Yi
G(l) = max(^l,... ,^n-l), = max(^j(xj+l), Yi),
\\Т!к=іскФ) - \{Y!k=ick)2V если x(i)<EUlCk<x(i + i),
О, в противном случае.
3. Вычисление критерия С(3). Расположим все узлы в порядке возрастания первой компоненты. Узлы
(х(к),у(к)), к =1,...,Ж, (12)
таковы, что х(1) = у(1) = 0, х(Ж) = у(Ж) = 1 и
х(1) < х(2) < • • • < х(Ж). (13)
Пусть x*(1) < x*(2) < ••• < x*(mx) и y*(l) < y*(2) < ••• < y*(my) —все различные значения среди абсцисс и ординат в (12), а l(i) —количество x*(i) в ряду (1З).
Так как критерии являются супремумами непрерывных функций, можем искать их вместо области К2 в К2 = [0,1)2, которую представим в виде объединения областей
Ajj = [x*(i),x*(i + 1)) * [y*(j),y*(j + 1)), 1 < i < mx - 1, 1 < j < my - 1.
Имеем
G(3) = max0ij, 0^ = sup |#jj(u, v)|, (u,v) = 0(u,v), (u,v) Є Ajj. (14)
i,j (u,v)eAij
u2v2
6ij(u, v) = —-----AjjMw + BijU + CijV - Dij, (15)
где
= Е СР’ = Е СР = Е СР = Е СР ж(Р У(Р)' (16)
р^^гз Р^^гз Р^^гз Р^^гз
Для любой пары (г, _?’) введем в рассмотрение множества номеров узлов
^г; = ик=11;(к) ^(к) = {Р : У(Р) < У*С?), в(к - 1) + 1 < Р < ^(к)} , (17)
где
к
£(к) = в(к), в(к) = 1(г), в(0) = 0. (18)
Г = 1
Если при одинаковых абсциссах соответствующие ординаты расположены в порядке возрастания, то
£(к) = в(к — 1) + шт(^’, /(к)), (19)
что позволяет сократить объем вычислений.
Множество ^ содержит номера тех узлов из набора (12), которые участвуют в формировании функции 0; (и, V) на Д;, а множество 1; (г) представим в виде
1;(*) = {Р : у(р) < У*СЛ, ж(р) = ж*(г)} .
При г = 0 в (16) положим Ао,; = 0, Во,; = 0, Со,; = 0, ^о,;- = 0. При г > 1
Аг,; = Аг — 1,; + Ср. (20)
ре/3- (г)
Аналогичные формулы имеют место для остальных величин в (16).
Алгоритм построения множества ^ аналогичен алгоритму, использованному при подсчете дискрепанса (см. [4-6]), в котором все точки входят с одинаковым весом 1/Ж. В данном случае нужно не только определить количество слагаемых в суммах (19), но и учесть различие весов Ск.
Найдем максимум |0; (и, -у)| в области
а < и < Ь, а = ж*(г), Ь = ж*(г + 1), (21)
С < V < Л, С = у*^), Л = у*(/ + 1), (22)
который совпадает с супремумом этой функции в области Д;.
Пусть при каком-либо ] и г € [го + 1, го + к] множества 1;(г) пустые, 1;(г) = 0.
Следствием этого являются равенства Аг0,; = Аг0+1,; = • • • = Аг0+к,;, и так как при этих
(г, _?) не меняются и коэффициенты Вг,; Сг,; Дг,;, остается постоянным и вид функции 0г,; (и, V). В этом случае ШаХ |0г,; (и, v)| следует искать сразу в области и^0=+0 Д^,;, т. е. по переменной и в области
а < и < Ь, а = ж*(го), Ь = ж*(го + к + 1), (23)
а по переменной V в области (22).
Так как вид функции |0^- (и, V) и область ее рассмотрения определены, опустим нижние индексы и перейдем к поиску ее максимума в общем виде. Пусть
2 2 и^2
в(и, V) = —---Аиу + Ви + Су — В, (24)
а < и < 6, с < V < й. (25)
Здесь максимум вычисляется аналитически. Ниже в п-мерном случае при п > 3 используется метод прямого перебора точек в области (25).
Максимум функции |0(и, V) может достигаться внутри области или на границе. Точки экстремума функции 0(и, V) при ВС = 0 находятся из кубического уравнения
Ви3 - 2АСи + 2С2 = 0,
получающегося из соотношений
22 и^ и V
6^(г(, у) = —-Ап + В = 0, 0'у(и, у) = —--Аи + С = 0.
Если точки экстремума попадают в область (25), то значения экстремума в этих точках следует сравнить с более просто определяемыми значениями максимума функции
|0(и, V)! на границах области (25). Случай ВС = 0 также является более простым. Этот алгоритм использован в приводимых ниже примерах.
4. п-мерный случай. Пусть функция /(их,... ,ип) задана в кубе К”. Так же, как в двухмерном случае, мы хотим оценить остаток |Л[/]| кубатурной формулы.
Пусть ск, 1 < к < N, —веса и
(жх(к), Ж2(к), .. ., жп(к)), 1 < к < N (26)
— узлы кубатурной формулы, упорядоченные по возрастанию первой координаты
жх(1) < жх(2) < • • • < жх^). (27)
Если в ряду (26) нет точек с координатами (0,0,... ,0) и (1,1,... ,1), то добавим их с нулевыми весами. Общее число точек считаем равным N.
Пусть при всех * (1 < * < п)
ж*(1) < ж*(2) < ••• < ж*(т;) (28)
— все различные значения *-ой компоненты узлов (26) и при всех ], 1 < ] < тх, 1(^’) — число элементов первой компоненты, равных ж|(^’).
Если кубатурные формулы точны для полиномов, т. е. выполнены равенства
1 N I
21 = ^2Ск П^1 ~х°р(к))’ 1^п> (29)
к=1 р=1
то при условии существования соответствующих производных остаток Л[/] кубатурной формулы в п-мерном случае можно записать в виде
ж, Г я2г+г / (й(/ ) 1)
т = Е(-1)' Ук, «-.№)■ (30)
Суммирование по г и 1 распространяется на все различные непересекающиеся наборы вг, в; номеров аргументов функции /, а именно,
вг = {^1, . . . , }, в; = { в 1, ..., в ;}, 1 < г < п, 0 < 1 < п — г, (31)
й(?г) = (и^ , ...,и^), и^ П иЯт = 0, 1 < к < г, 1 < т < 1. (32)
В формуле (30) производная функции / берется в точке, у которой все компоненты,
кроме компонент и(вг), равны единицам. Функция Ф(и(вг); в;) = Ф(и41,. .., и4г; в;) имеет вид
1 / 1 г N ; г \
Ф(й(^);в,) = ^ I — П и1т ~ 21^2ск П(1 ~х°Лк)) П ~х^{к)) I , (33)
\ т=1 к=1 р=1 т=1 /
и при оценке остатка |Д[/]| кубатурной формулы в п-мерном случае нужно найти супремумы модулей этих функций
С(вг ;в;)= 8Ир |Ф(ип ,...,и4г ;в;)|, (34)
0<^1 ,...,иг<1
которые при выполнении равенств (29) могут быть записаны в виде
N г
С(и;в1)= эир 1 0<^1 < 1 2
1
^7 П П ~х^(к))
т=1 к=1 т=1
, 1 < г < п, (35)
где
; N
= 2; с*П(1 — ж«Р(к)), = 1 (36)
р=1 к=1
и Ск = Ск при 1 = 0.
При г = п(1 = 0) имеем в” = {1, 2, ...,п}, во = 0 и критерий С(вп; 0) является аналогом критерия С(3) в двухмерном случае:
С(в”; 0) = эир |Ф(и1,... ,и” 0)|, (и1,...,и”) е К”, (37)
и1,...,ип
где
1 ” N ”
ф(«ь■••.,1»;0) = рП,‘‘ -_ж»(*0)- (38)
®=1 к=1 ®=1
Рассмотрим вычисление критерия С(вп; 0), так как остальные случаи С(вг, в;), 1 < г < п — 1, получаются заменой п на г при фиксированных значениях соответствующих переменных (и* = 1) в количестве п — г штук. Если г + 1 = п, то С(вг; в;) < С(вп; 0).
Пусть КТ” = [0,1)” = иДй1 , где
”
Дк1 = П[ж* (к®), ж* (к* + 1)), 1 < к < т; — 1, 1 < * < п. (39)
®=1
Функция Ф(и1,..., ип; 0) в области Дд^,...,^ является многочленом относительно произведений и*. Покажем, как найти коэффициенты при соответствующих степенях.
Множества ^и (кі) (аналог множеств (17) в двухмерном случае) опре-
деляются как
= иг = і/^2 (г), ^0,Й2 = 1Й2,-..,Йп (0) = 0, (40)
Ґ Й1-1 Й1 Ї
Ік2,...,кп(кі) = < р : X(р) < ж*(к), 2 < і < п, ^ 1(г) + 1 < р < ^ ((гН . (41)
I Г=1 Г=1 )
Номера узлов, образующих множество /^2,...,&„(кі), таковы:
ІЙ2,...,кп(кі) = {р : Хі(р) < ж*(к*), 2 < і < п, жі(р) = жі(кі)} . (42)
Формулы (40)—(42) поясняют алгоритм и удобны для программирования.
В области (39) функция (38) записывается в виде
-і П П
^,...,*>1, 0) = - п«I - Е Е(-1)'“| П^к1^ (43)
і= і РЄ^1,...,Л:п а І=і
где введены мультииндексы
п
а = (аь «2,. .., ап), «і = {0,1}, |а| =Е «і, (44)
і=і
а суммирование ^ проводится по множеству всех 2п значений а.
Если воспользоваться обозначением
п
А1,...,к„ = Е срП ж“* (р), (45)
р^^кі,...,кп і і
то в формуле (43) коэффициент при ПП=і иі, получающийся при а = 0 = (0,..., 0) и являющийся аналогом в двухмерном случае, можно найти по формуле
^к?,...,*» = Е СР (46)
или с помощью рекуррентного уравнения
Ак0)...,кп = Ак0-і,...,й„ + Е ср (47)
РЄІк2,...,кп (кі)
при А00)2 к =0. Аналогично находятся коэффициенты при других значениях а. Как и в двухмерном случае, если при каком-либо І0
^й2 ,...,&„ (і0 + 1) = • • • = /й2 ,...,&„ (і0 + г) = 0
то в области Д = и^О=+0 Д^,к2,...,к„ функция Ф^1 ,...,&„(«і,..., «п; 0) не меняет своего вида, следовательно, ее максимум следует искать сразу в области и^О=+0 Д^,й2,...,йп.
Критерий С(іп; 0) представляет собой максимум супремумов модулей функций
Ф&1 ,...,&„(«і,... ,«п; 0) в областях Д^,...,^ или ^=0 Д^,к2,..,к„.
Ясно, что для модуля остатка кубатурной формулы справедливо неравенство
д2г+;/(и(в), 1)
|Д[/]| <ЕС(в ;*)
, О К г
г,1
ди2(£г )ди( в ;)
йи(£ г). (48)
Введем весовые коэффициенты 7Г; (см.[7]) и, используя взвешанные оценки С(£г; в ;), получим модуль остатка кубатурной фомулы в виде
|Д[/]| < шах ( —
г,; \7rUK г
\7rUKг ди2(^)ди( вг) у
Замечание 4. Свойства, описанные в замечании 3, распространяются и на п-мер-ный случай. Веса и узлы (26), входящие в кубатурную формулу, у которых хотя бы одна из координат равна единице, не участвуют в вычислении критерия С(вг; в;). При разложении функции /(X) в ряд в точке (0,..., 0) (см. [3]) соответствующий критерий С не содержит весов и узлов, имеющих хотя бы одну нулевую координату.
Замечание 5. Рассмотрим оценку остатка кубатурной формулы, когда не выполнены условия (29). В этом случае остаток Д[/] имеет вид
д[/] = Е(-!)г 1кг ^д¥ЩЩ^ф{й{и); *)<вд-
п /1 N \ п /1 N \
_Е (■*■) ( 2 _ Е _ Ж*(^)) ) + Е (Ч и” Е - ж*(к))(1 - Ху (к)) | .
*=1 V к=1 / *.з = 1 V к=1 /
г=3
(50)
Запишем выражение для функции Ф(и(вг); в;) в области :
1 1
"1
*=1 ..........................к„ а *=1
Фй1,...,й„(м(^);«0 = я(«г) ( “ Е срЕ(-1)1/3'П^^г'1 *
(51)
где
N1 1 ; N
я(^) = ЕСйП^1 с~* = ~н<Т)СкП^1 -х°р(к))> Ес~* = 1-
к=1 р=1 ( в1) р=1 к=1
Здесь в — аналог а в формуле (43), а |в| = $^Г=1 в*. Производные /^(1) и /Х^.(1) берутся в точках с единичными координатами. В правой части равенства (51) выражение, стоящее в скобках, отличается от соответствующего выражения в формуле (43) только первым слагаемым, а остальные слагаемые находятся по прежней схеме, поэтому алгоритм вычисления супремума функции Ф(м(1г );в;) с очевидными изменениями повторяет алгоритм, описанный выше.
Пример 1. Найдем величины критериев С(вг; в;) при п = 3, (и, V, ад) € К3, N = 40.
Пусть первые 8 узлов кубатурной формулы расположены в вершинах куба К3, а со-
ответствующие им веса равны С1 = С2 = • • • = С8 = 1/120. На каждом ребре помимо
угловых точек (которые зафиксированы раньше) располагаются два узла, делящие ребро на 3 равные части. Этих точек 24 и они имеют одинаковые веса сд = • • • = С32 = 1/60. Остальные 8 точек имеют веса С33 = • • • = С40 = 1/15, а сами узлы находятся в вершинах куба со стороной 1/3, расположенного симметрично внутри исходного куба.
Критерий С(вз; во) принимает вид
С({1, 2, 3}; 0) = яир
0<и,^,ад<1
2 2 2 N
и2^2^2
8
_ 53 с* ^(м — х(к))Ь(-у — у(к))Ь(ад — г (к)).
к=1
(53)
Соотношения (29) выполнены. Так как компоненты узлов равноправны в смысле их расположения и весов, приведем величины критериев С(вг, в;) при наборах векторов с ?1 = {1}, в2 = {1, 2}, вз = {1, 2, 3} и соответствующими в;. В этом и последующих примерах запишем величины критериев, опуская фигурные скобки аргументов функции С, например, С(1;2, 3) вместо С({1}; {2, 3}):
1) С(1, 2, 3; 0) = 0.00632;
2) С(1, 2; 0) = 2 С(1, 2;3) = 0.01064;
3) С(1; 0) = 2 С(1;2) =4 С(1;2, 3) = 0.01389.
При разложении функции /(X) в ряд Тейлора в окрестности точки (0, 0,0) приходим к оценке критериев С?(вг,в;) (см. (8)), для которых, ввиду выполнений равенств типа (9), при всех аргументах имеют место равенства С?(вг.; в;) = С(вг; в;).
Пример 2а. В [7] узлы кубатурных формул образуются по следующей схеме: ж* (к) = {кМ*/Т}, 1 < к < Т, где фигурные скобки обозначают дробную часть числа, М* —целые числа взаимно простые с Т и М* < Т, а веса с* = 1/Т. Полагаем Т = 40, М1 = 7, М2 = 11, М3 = 19, добавляем точку с единичными координатами х* (41) = 1, С41 =0, N = 41. Для исходных данных соотношения (29) не выполняются. Получены следующие значения критерия:
1) С(1, 2, 3; 0) = 0.00719.
2) С(1, 2; 0) = С(1, 3; 0) = 0.00433, С(2, 3; 0) = 0.02547.
3) С(1; 0) = С(2; 0) = С(3; 0) = 0.01250.
4) С(1, 2; 3) = 0.00719; С(1, 3; 2) = 0.00552; С(2, 3; 1) = 0.00719.
5) С(1; 2) = 0.00391; С(1;3) = 0.00223; С(1; 2, 3) = 0.00516.
6) С(2; 1) = 0.00223; с(2; 3) = 0.02594; с(2; 1, 3) = 0.00719.
7) С(3; 1) = 0.00391; с(3; 2) = 0.02594; с(3; 1, 2) = 0.00590.
Пример 2Ь. Алгоритм построения узлов и весов тот же, что и в примере 2а.
Если взять Т = 40, М1 = 7, М2 = 23, М3 = 29, N = 41, то величина многих соответствующих критериев отличается на порядок:
1) С(1, 2, 3; 0) = 0.04578,
2) С(1, 2; 0) = 0.02547; 3) С(1, 3; 0) = С(2, 3; 0) = 0.03859.
3) С(1, 2; 3) = С(1, 3; 2) = С(2, 3; 1) = 0.04473, С(1;2) = С(2; 1) = 0.02594.
4) С(1; 3) = С(2; 3) = С(3; 1) = С(3; 2) = 0.03906.
5) С(1; 2, 3) = С(2; 1, 3) = С(3; 1, 2) = 0.04578.
В последних двух примерах вне зависимости от содержния векторов 1 и 2 имеем Н( в 1) = 0.51250, Н(в2) = 0.28906.
Использование критериев дает возможность подобрать оптимальный набор узлов и весов кубатурных формул.
Автор благодарит С. М. Ермакова за постановку задачи и весьма ценное обсуждение.
1. Крылов В. И. Приближенное вычисление интегралов. М.: Наука, 1967.
2. Соболев С. Л. Введение в теорию кубатурных формул. М.: Наука, 1974.
3. Ермаков С. М., Бурнаева Э. Г., Сидоровская М. В. Многокритериальные методы анализа остатка кубатурных формул // Математические модели. Теория и приложения. Вып. 7. СПб., 2006.
4. Товстик Т. М. Вычисление дискрепанса конечного числа точек в n-мерном единичном кубе // Вестн. С.-Петерб. ун-та. Сер. 1, 2007. Вып. 3. С. 118-121.
5. Bundschuh P., Zhu Y. The formulas of exact calculation of the disgrepancy of low-dimensinal finite point sets(I) // Chinese Sci. Bul. 1993. Vol. 38. N15. P. 1318-1319.
6. Bundschuh P., Zhu Y. The formulas of exact calculation of the disgrepancy of low-dimensinal finite point sets(II) // Chinese Sci. Bul. 1995. Vol.40. N7. P. 610-612.
7. Sinescu V., Joe S. Good latice rules with a composite number of points based on the product weighted star discrepancy // Monte Carlo and Quasi-Monte Carlo Methods 2006. P. 645-658. Springer, 2007.
Статья поступила в редакцию 18 сентября 2008 г.