Алгоритм Добровольской и численное интегрирование с правилом остановки Текст научной статьи по специальности «Математика»

ЧЕБЫШЕВСКИЙ СБОРНИК Том 10 Выпуск 1 (2009)

УДК 511.9

АЛГОРИТМ ДОБРОВОЛЬСКОЙ И ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ С ПРАВИЛОМ ОСТАНОВКИ 1

Е. Д. Ребров (г. Тула)

Аннотация

В работе рассматриваются вопросы численного интегрирования периодических функций многих переменных по системе квадратурных формул с парралелепипедальными сетками, построенных по алгоритму Л. П. Добровольской.

Введение

В 1957 - 1960 годах ([7] — [9]) при создании теоретико-числового метода в приближенном анализе Н. М. Коробов ввёл в рассмотрение широкий класс периодических функций Е<а (а > 1) с быстро убывающими коэффициентами Фурье, состоящий из функций f (х\, ..., х3), имеющих по каждой из переменных Х\, ... ,х3 период, равный единице, и для которых их ряды Фурье

удовлетворяют условиям

где константа С те зависит от ш\,..., т3, и для веществе иных ш полагаем ш = тах(1, |т|). Ясно, что такие ряды Фурье сходятся абсолютно, а поэтому для любого (а > 1) они представляют непрерывные функции.

Относительно нормы

пространство является несепарабельным банаховым пространством изоморфным пространству — всех ограниченных комплексно значных последовательностей (см. [5]).

ГО

||/(жь ... ,ж5)||е? = вир\С(т)\(ті. ..т3)а

1Работа выполнена по гранту РФФИ 08-01-00790

Усеченной нормой называется величина q{x) х\ ... где для вещественного х обозначаем х = тах(1, |ж|). Тогда усеченной норменной поверхностью с параметром t называется множество N(t) = {x|q(x) = t}. Для натурального

t на усеченной норменной поверхности имеется r*(t) целых ненулевых точек, 2

ГДР

Ts(t) = ^ 1 (3)

m eN (t)

— число представлений натурального числа t в виде t = тГ... т,

Используя новые обозначения, можно написать другое выражение для нормы функции Ilf (x)||Ea. Справедливо равенство

Ilf(x)||Ea = supta max |C(m)|.

s t^i meN (t)

Через Ea(C) обозначается множество функций из Ef с нормой, не превосходящей C, то есть шар в банаховом пространстве Ef радиус a C с центром в нуле.

Нетрудно видеть, что произвольная периодическая функция f (x) из E^(C) по модулю ограничена величиной C (1 + 2((a))s, при этом данная оценка достижима на функции

°° c

f(x)= V т=----------- _ч e2m(tn’g)

^ (mi ■ ■ ms)a

m =—оо

в точке X = 0, Здесь и далее, как обычно, ((а) — дзета-функция Римана.

Очевидно, что Ea(C) С Ee(C) при а ^ в- Для любой периодической функции f (X) Є E3a(C) С Ee(C) справедливо неравенство для норм

llf(X)|k“ ^ llf(X)||Ef-

Равенство достигается только для конечных тригонометрических многочленов вида,

f (X)= c(m) e2m(m’g).

m gN (i)

Рассмотрим квадратурную формулу с весами

i i N

\ \ f(xi,... ,xs)dx1.. .dxs = ^ ^2 pkf [Ci(k),-■■ ,Cs(k)] - RN[f]. (4)

N

о о к=1

Здесь через ] обозначена погрешность, получающаяся при замене интеграла

1 1

!-!!

оо

23десь и далее ^ означает суммирование по системам (mi,..., ms) = (0,..., 0).

средним взвешенным значением функции f (х1,... ,х3), вычисленным в точках

Мк = (Ш,...,Ш) (к = 1 ..^).

Совокупность М точек Мк называется сеткой М, а сами точки — узлами квадратурной формулы. Величины рк = р(Мк) называются весами квадратурной формулы. В этой работе будем везде предполагать, что все веса вещественнозначные.

Для произвольных целых т1,...,т8 сумм ы Бм^(т1,.. .,т3), определённые равенством

N

Бм,Р<т1,.. .т) = 5]Ркв2пг[т^(к)+-+т(к)\ (5)

к=1

называются тригонометрическим,и суммами сетки с весам,и.

Будем, также, рассматривать нормированные тригонометрические суммы сетки с весам,и

Зм,р(ть- • •’т«) = ■^13мАть- ■ -,т3).

N

Положим р(М) = ^2 \рз |, тогда для всех нормированных тригонометриче-

3 = 1

ских сумм сетки с весами справедлива тривиальная оценка

|5^й)| < ^Р(М).

1

ма сетки и писать Бм (т) и нормированная тригонометрическая сумма сетки Б*м(т).

Обобщая работу [6], дадим следующее определение.

Определение 1. Дзета,-функцией сетки М с весам,и р и параметром р ^ 1 назовем, функцию ((а,р\М, р) заданную в правой полуплоскости а = а + й (а > 1) рядом Дирихле

= — го у 7 п=1

где

3"(р,М,р,и)= £ \^м)|р. (7)

т (п)

Непосредственно из определения следует неравенство

с(ра,р\М,р) ^ Ср(а.1|М,р). (8)

1М параметром р и писать ((а,р\М) .

Справедливы следующие две обобщенные теоремы Коробова о погрешности квадратурных формул.

Теорема 1. Пусть ряд Фурье функции / (X) сходится абсолютно, С(т) — ее коэффициенты Фурье и 3р(т) — тригонометрические суммы сетки с весами, тогда справедливо равенство

/ \ ^

Як [/] = С(0) ( ^,,„(6) - 1) + 1 С(т)Ям,Дт) =

' ' т-1 т_ — — с-*-',

т1,...,тв = -х сю

= С(0) (31,,Д0) - 1) + £’ С(т)ЙМ,ДЙ)

Ш1,...,Шв = — ^

(9)

и при N ^ ж погреитость Ям [/] будет стрем,и,ться, к нулю тогда, и только тогда, когда, взвещенные узлы квадратурной формулы равномерно распределены в единичном в-мерном, кубе.

Теорема 2. Если, /(хг,... ,х3) € Е<а(С), то для погрешности квадратурной формулы справедлива, оценка

\Ям[/]\< С

0) -1

с

N

£'

\3м,р(т )\

С

Ш1,...,Ше= — <^

3мА0) - 1 + С ■ с (а, 1\М,р)

(т 1... т3)с

(10)

где сумма 3м,з(т) определена, равенством (5). На, классе Е<а(С) эту оценку нельзя, улучшить.

Модифицированной сеткой М (в) будем называть сетку, состоящую из модифицированных узлов

Мк (в) = Ш*) + вг},..., (&(к) + вЛ) (к =1 ).

Линейный функционал погрешности приближенного интегрирования периодической функции / (X) из класеа Е^ то квадратурной формуле с весами р и модифицированной сеткой М (в) будем обозначать через Ям(в) д[/ (X)], а его

Я„

норму — через

ьм (3),р

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

Я

м (3),3

Е?

~^8м,р{0) - 1

1

+ N

£'

Ш1,...,тз=—<ж

I £м,/>(гаь... ,т3 (тГ... т1)а

3м,з(0) — 1 + С(а, 1\М,Р)'

(н)

Из этого равенства следует, что для всех модифицированных сеток М (в)

М

ции. Это просто объяснить через понятие граничной функции класса, которое впервые встречается в работе Н. М. Коробова [11], а более подробно в его монографии [12] (см. также [1]).

Функции /(X) с единичной нормой, для которых абсолютная погрешность приближенного интегрирования равна норме линейного функционала погрешности, следуя Коробову, будем называть граничными функциями класса, Е^. Легко видеть, что граничной функцией класса Е^ для сетки М с весами р будет функция с коэффициентами Фурье

Ясно, что если при некоторых значениях ш1; ,,,, ш3 тригонометрическая сумма Бм,р(Ш) = 0, то граничная функция класса сеткой определена неоднозначно.

Нетрудно видеть, что если f (X) — граничная функция класса Еа для сетки

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

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

3

по алгоритму Добровольской.

Вопросы численного интегрирования с правилом остановки рассмотрены в работе [3]. Следуя этой работе дадим частное определение мультипликативной дискретной дисперсии для случая произведения двух параллелепипедальных сеток с равными весами.

М

дением, двух параллелепипедальных сеток 3

при Бм,р(ш) = 0, при Бм,р(ш) = 0.

(12)

М, то д(Х) = f (X — 3) — граничная функция класса Еа для модифицированной сетки М (33),

М = Мі ■ М2 = {{X + у} IX Є Мі, у Є М2} ,

3Для любого г Є М8 полагаем {г} = ({гі},..., {^я}) Є [0; 1)я.

то мультипликативной дискретной дисперсией погрешности приближенного

М

риодической функции / (X) из пространства, Еа будем называть величину 0*м м [/(X)], заданную равенством

Из определения видно, что, вообще говоря, 0*м м [/(X)] = В*М2.м1 [/(X)],

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

1 Программная реализация алгоритма численного интегрирования с правилом остановки

В этом разделе мы рассмотрим программную реализацию алгоритмов численного интегрирования периодических функций многих переменных. Как показали И, И, Чепцов и И, Ф, Шарыгин вопрос об интегрировании произвольных функций с помощью периодизации задачи сводится к случаю периодических функций. Поэтому, для простоты изложения мы рассматриваем только классический случай теоретико-числового метода Коробова, то есть интегрирование периодических функций многих переменных.

Система компьютерной математики МаЛсас114 выбрана по нескольким причинам.

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

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

1.1 Алгоритм Коробова вычисления оптимальных коэффициентов

Следующая программа вычисляет по заданному модулю р и набору коэффициентов а = (а0, аг,..., а3)Т значение функции Н(р,а), которая является количественным критерием качества набора коэффициентов с точки зрения

2 2

при \М\ ^ ж будем иметь 0*м м [/(X)] ^ 0,

применимости для использования в квадратурных формулах с параллелепп-педальными сетками.

HT(p, a) :=

p

p-i s s

E 1—2

k=0 j=0

a,j - к P

a,j - к P

Например, при p = 101 и a = (1,19, 85)T получим

HT (p,a) = 1.1030731532962952.

Ясно, что трудоемкость вычисления по этой программе равна O (s • p) элементарных арифметических операций.

При реализации произведения двух параллелепипедальных сеток естественно иметь программу вычисления точек сетки в виде двумерного массива, что и делает программа Setka(p, а).

Setka(p, а) := s •*— length(a) — 1, S'p-i,. for к £ 0..p — 1 forj £ 0..s

Sktj - ^ ~ floor (^)

S'

Эта программа также требует для формирования двумерного массива координат точек сетки O (s • p) элементарных арифметических операций.

При наличии массива координат точек параллелепипедальной сетки есте-

H

HS(S)

HS(S) :=

р

p-i s 2

E П (1 — 2 • Sk,j)

k=0 j=0

R

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

HT 1(p,s) :=

as-i ^ 1, b0 ^ 1, a0 ^ 1 forj Є 1..s — 1 R <- 3J+1 for c Є 1 ..p — 1 bj <- с

~ , 3i+1

p-i j

E П

k=0 .1=0

-i 2‘

bi-k

P

bi-k

P

2

p

a

Необходимо отметить, что первоначально Н, М. Коробов свой алгоритм обосновал только для простых модулей. Позднее он его модифицировал для составных модулей равных произведению двух простых, что позволило сократить трудоемкость алгоритма вычисления оптимальных коэффициентов с О (М2) до

О Хотя в литературе не встречается обоснование этого алгоритма для

произвольного составного модуля, но не трудно понять, что такое обоснование, но без оценок качества, нетрудно дать. Действительно, из результатов разных авторов следует существование оптимальных коэффициентов по любому составному модулю. Кроме того, известно, что ^-мерный набор оптимальных коэффициентов всегда можно дополнить до 5 + 1-мерного набора. Поэтому, если среди всех 5 + 1-мерных наборов с фиксированными первыми 5 коэффициентами взять набор с минимальным значением функции Н — количественной меры качества оптимальных коэффициентов, то этот набор и будет оптимальным.

1.2 Алгоритм Добровольской вычисления оптимальных коэффициентов по составному модулю

Алгоритм Л, П, Добровольской [2] основан на двух идеях.

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

Во-вторых, хотя алгоритм работает для любого составного модуля, но своих оптимальных характеристик он достигает, когда модуль является произведением специальной последовательности простых чисел.

Следующая программа по заданной сетке Б из Р точек, представленной двумерным массивом, находит оптимальную параллелепипедальную сетку из р точек, такую, что их произведение будет оптимальной параллелепипедальной сеткой из Р ■ р точек.

HTS(S,p,s) :=

a,s—i — 1, bo — 1, ao — 1, P — rows(S) for j £ 1..s — 1 ' R<- 3J+1 for c £ 1..p — 1

bj <- с

oi+1 г <— p-P P—1 p—1 E E n = 0 k = 0 Д [l - 2 (sn,i + ^ - floor (sn>1 + ^))] 2'

aj — c, R « — r if r < R

°j aj SP-p — 1,s-1 — 0

for k £ 1..p — 1

forn £ 0..P — 1 for j £ 0..s — 1

Sp.k+nJ <- SnJ + ^ - floor (snJ +

S

Нетрудно видеть, что программа HTS(S,p,s) требует для своего выполнения память объёмом O (s • P • p) и O (s2 • P • p2) элементарных арифметических

операций.

Важной особенностью данной программной реализацией является тот факт,

р

. гы 1**1 из которых суть модифицированная исходная параллелепипедальная сетка. Эта особенность потом будет использована при реализации алгоритма численного интегрирования с правилом остановки.

Перейдем к описанию второй составляющей алгоритма Добровольской. р

простое число непревосходящее 100, Будем говорить, что р1,р2,... ,рк — допустимая последовательность простых чисел первого типа длины к для простого р

р

существует допустимая последовательность простых чисел первого типа прок

Заметим, что р2 ^ 5, поэтому любая допустимая последовательность простых р1, р2,... ,рк первого типа строго монотонно возрастает.

Непосредственными вычислениями, используя таблицу простых, нетрудно проверить, что для простого р =3 допустимой последовательностью простых первого типа длиной 5 будет последовательность 3, 5,11, 61,1847, которая определяет модуль N = 18582667.

р р1, р2, . . . , рк к

последовательность р1,р2,... ,рк и одну макспмальную р1 ,р2,,рк, которые обладают свойством

р1, р2, . . . , рк

и тем же р. Тем самым определяется минимальный модуль Мрк = р1 ■ р2 ■... ■ рк

вого типа длины 5 будет последовательность 3, 5,11, 43, 919, а максимальной

- 3,5,13,97, 4463. Соответственно, М 5 = 3 ■ 5 ■ 11 ■ 43 ■ 919 = 6704105 и М;5 = 3 ■ 5 ■ 13 ■ 97 ■ 4463 = 84701370.

Пусть р > 3 — фиксированное нечетное простое число, например, любое нечетное простое число отличное от 3 и непревосходящее 100. Будем говорить, р1, р2, . . . , рк кр выполнены условия

(2 $ і $ к).

(14)

р) < Рі <р) (3 $ і $ к)

(15)

р

существует допустимая последовательность простых чисел второго типа прок

Следующая программа НТБ 1(рр, 5) то заданному вектору рр = (р0,... ,рк)т, где р0, ,,,, рк допустимая последовательность простых первого или второго типа, строит оптимальную параллелепипедальную сетку из Р = р0 ■ ... ■ рк точек,

НТБ1(рр, в) := к *— 1епдЬк{рр) — 1

а ^ НТ 1(рр&, ®), Б ^ БвЬка^рри, а)

/ог ] € к — 1, к — 2..0

Б ^ НТБ(Б,ррз ,в)

Б

Нетрудно видеть, что количество операций требуемых для выполнения этой программы для последовательности простых первого типа есть

О Ш + рк ■ р\-1 + ... +Рк ■ рк-1 ■ ... ■ р1 ■ р2о)) = к

O

Y.^k + p-p°) 1 = О (s'2 ■ (Р ■ (4+р„)) .

j=1

2j

1.3 Алгоритм интегрирования с правилом остановки

Следующая ниже программа реализует алгоритм численного интегрирования с правилом остановки. Предполагается, что сетка Б вычислена по алгоритму Бочаровой, Тогда пользуясь тем, что двумерный массив, задающий сетку, специально организован, данная программа вычисляет приближенное значение интеграла от функции / (х), наращивая количество точек от пер вой сетки из рк

Р

IntPo(pp, S, f, s) :=

s — cols(S) — 1, k — length(pp) — 1 x — submatrix(S, 0, 0,0, s)T, M — f (x), m — forn £ 0..ppk — 1

x 4— submatrix(S, n, n, 0, s)^, r f(x)

M — r if r > M m — r if r < m otherwise Int — Int + r Int <— D <— 2 • e, j <— к — 1, P <— рр^

PPk k

while (D > e) ■ (j > 0)

D <— Int2 for l £ 1..ppj — 1 ' In-0 forn £ 0..ppk — 1

Pn — P • l + n,x ■

M, Int 0

submatrix(S, Pn, Pn, 0,s)T ,r — f (x)

Ir

(Int

Int

D

M

m— Ir _ Ir_

P

r if r > M r if r < m otherwise - Ir + r

D — D + Ir2, Int — Int + Ir

_ 1 nt j-y

ppj >

m M)

-2- -Int2,P ■ ppj

P • ppj ,j — j — 1

Программа в качестве результата выдает приближенное значение интеграла, величину мультипликативной дисперсии, а также минимальное и максимальное

значения функции на сетке.

В качестве тестирующей предлагается функция

__ 0^+1

f (х) = 3

которая реализуется программой

Ї(х) :=

П(1 — 2

3=0

Ц (1 — 2 ■ хі) .3=0

СПИСОК ЦИТИРОВАННОМ ЛИТЕРАТУРЫ

[1] Бочарова (Добровольская) Л. П. О граничных функциях некоторых классов // Наукоемкое образование. Традиции. Иновации, Перспективы. Сборник межвузовских научных статей. Тула, АНОВО "ТИНО 2006. С. 198 - 202

[2] Бочарова (Добровольская) Л. П. Алгоритмы поиска оптимальных коэффициентов // Чебышевекий сборник 2007 Т. 8. Вып. 1(21). Тула, Из-во ТГПУ им. Л. Н. Толстого. С. 4 — 109.

[3] Добровольская Л. П., Добровольский Н. М., Симонов А. С. О погрешности приближенного интегрирования по модифицированным сеткам / / Чебышевекий сборник 2008 Т. 9. Вып. 1(25). Тула, Из-во ТГПУ им. Л. Н. Толстого. С. 185 - 223.

[4] Добровольский М. Н,, Добровольский Н. М., Киселева О. В. О произведении обобщенных параллелепипедальных сеток целочисленных решёток // Чебышевекий сборник 2002 Т. 3. Вып. 2(4). Тула, Из-во ТГПУ им. Л.Н.Толстого. С. 43 — 59.

[5] Добровольский Н. М., Манохин Е. В. Банаховы пространства периодических функций // Изв, ТулГУ, Сер. Механика. Математика. Информатика. Т. 4. Вып. 3. Тула, 1998. С. 56-67.

[6] Добровольский Н. М., Манохин Е. В., Реброва И. Ю,, Рощеня А. Л. О непрерывности дзета-функции сетки с весами // Известия ТулГУ. Сер. Математика. Механика. Информатика. Т. 7. Вып. 1. Тула, 2001. С. 82-86.

[7] Коробов Н. М. Приближенное вычисление кратных интегралов с помощью методов теории чисел // ДАН СССР. 1957. N 6. С. 1062 - 1065.

2

[8] Коробов Н. М. Вычисление кратных интегралов методом оптимальных коэффициентов // Вестн. Моск. ун-та, 1959. N 4. С. 19 — 25.

[9] Коробов Н, М. О приближенном вычислении кратных интегралов // ДАН СССР. 1959. Т. 124, N 6. С. 1207 - 1210.

[10] Коробов Н. М. Теоретико-числовые методы в приближенном анализе. М.: Физматгиз, 1963.

[11] Коробов Н. М. Квадратурные формулы с комбинированными сетками // Математические заметки. 1994. Т. 55. Вып. 2. С. 83 — 90.

[12] Коробов Н. М. Теоретико-числовые методы в приближенном анализе, (второе издание) М.: МЦНМО, 2004.

Тульский государственный университет Получено 10.04.2009