close

Вход

Забыли?

вход по аккаунту

?

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

код для вставкиСкачать
Известия Тульского государственного университета
Естественные науки. 2013. Вып. 3. С. 193–201
Прикладная математика и информатика
УДК 511.3
Алгоритмы численного интегрирования
с правилом остановки
Н. К. Серегина
Аннотация. В работе рассматриваются алгоритмы численного
интегрирования с правилом остановки для параллелепипедальных
сеток. Приводятся программы в системе MATHCAD15, реализующие
эти алгоритмы.
Ключевые слова: параллелепипедальная сетка, алгоритм с
правилом остановки, квадратурные формулы, метод оптимальных
коэффициентов, количественная мера качества сетки.
Введение
Актуальность и цель работы. Актуальность темы состоит в том,
что построение и исследование алгоритмов численного интегрирования с
правилом остановки являются важным современным направлением в теории
теоретико-числовых методов приближенного анализа (см. [4]).
Цель работы — подробно рассмотреть классический случай теоретикочислового метода Коробова оптимальных коэффициентов интегрирования
периодических функций многих переменных, а также дать программную
реализацию алгоритмов численного интегрирования с правилом остановки
для последовательности концентрических параллелепипедальных сеток.
Как показали Н. Н. Ченцов и И. Ф. Шарыгин вопрос об интегрировании
произвольных функций с помощью периодизации задачи сводится к случаю
периодических функций. Поэтому рассматриваем только классический
случай теоретико-числового метода Коробова — интегрирование
периодических функций многих переменных. В работе предлагаются
соответствующие программы в системе компьютерной математики Mathcad15, реализующие данные алгоритмы численного интегрирования.
Параллелепипедальные сетки Н. М. Коробова. В 1959 году
профессор Н. М. Коробов предложил новый класс теоретико-числовых сеток
— параллелепипедальные сетки:
¾
½
¾¶
µ½
as k
a1 k
,...,
(k = 1, 2, . . . , N )
Mk =
N
N
Н. К. Серегина
194
и соответствующие квадратурные формулы с равными весами
Z1
Z1
...
0
0
¾
½
¾¶
N −1 µ½
1 X
a1 k
as k
f (~x) d~x =
f
,...,
− RN [f ],
N
N
N
k=0
где RN [f ] – погрешность квадратурной формулы.
На классе Esα периодических функций с быстро сходящимися рядами
Фурье были получены наилучшие результаты
lnα(s−1) N
(Н. С. Бахвалов, Н. М. Коробов).
Nα
Качество параллелепипедальной сетки. Количественной мерой
качества набора коэффициентов a0 , a1 ,. . . , as параллелепипедальной сетки
называется величина
½
¾¶2
p−1 s µ
aj · k
3s+1 X Y
1−2·
,
(1)
H(p, ~a) =
p
p
|RN [f ]| <<
k=0 j=0
которая равна приближенному значению интеграла от функции h(~x) =
s
s+1 Q
= 3p
(1 − 2xj )2 по квадратурной формуле с параллелепипедальной
j=0
сеткой
1=
Z
1
...
0
Z
0
1
½
¾¶2
p−1 s µ
aj · k
3s+1 X Y
h(~x)d~x =
1−2·
− Rp [h],
p
p
k=0 j=0
где Rp [h] — погрешность приближенного интегрирования.
Выбор функции h(~x) и величины H(p,³~a) связан
с тем, что функция h(~x)
´
π2
α
является граничной функцией класса Es ·, 6 (подробности см. в [10]).
Вычисление количественной меры качества набора коэффициентов. Следующая программа вычисляет по заданному модулю p и набору
коэффициентов a = (a0 , a1 , . . . , as )T значение функции H(p, a), которая
является количественным критерием качества набора коэффициентов
для использования в квадратурных формулах с параллелепипедальными
сетками:
HT (p, a) :=
s ← length(a) − 1
·p−1 · s h
³
³
´´i¸2 ¸
P Q
s+1
a ·k
a ·k
R← 3p ·
1 − 2 · jp − f loor jp
k=0
j=0
R
Например, при p = 101 и a = (1, 19, 85)T получим
HT (p, a) = 1.1030731532962952.
Ясно, что трудоемкость вычисления по этой программе равна O (s · p)
элементарных арифметических операций.
Алгоритмы численного интегрирования с правилом остановки
195
Вычисление точек сетки. При реализации произведения двух
параллелепипедальных сеток естественно иметь программу вычисления
точек сетки в виде двумерного массива, что и делает программа Setka(p, a):
Setka(p, a) :=
s ← length(a) − 1, Sp−1,s ← 0
f or k ∈ 0..p − 1
f or j ∈ 0..s
³
´
a ·k
a ·k
Sk,j ← jp − f loor jp
S
Эта программа также требует для формирования двумерного массива
координат точек сетки O (s · p) элементарных арифметических операций.
Следующая программа:
Setka1(p, a) :=
s ← length(a) − 1, Sp−1,s ← 0
f or k ∈ 1..p − 1
f or j ∈ 0..s
r ← Sk−1,j + aj
Sk,j ← if (r < p, r, r − p)
S
S←
p
Программа Setka1(p, a) при условии, что выполнены неравенства 1 6
6 aj < p (j = 0, . . . , s) обеспечивает более точное вычисление координат
точек при больших значениях p, так как операция деления выполняется
один раз над всем массивом, а при вычислении числителей координат все
вычисления ведутся с числами в диапазоне от 0 до 2p и не происходит потери
точности из-за вычитания двух плавающих чисел из разных диапазонов
значений, которое встречается в первой программе.
Вычисление количественной меры качества сетки. При наличии
массива координат точек параллелепипедальной сетки естественно иметь
программу, которая вычисляет значение функции H непосредственно по
координатам точек. Таким образом получаем программу HS(S):
HS(S) :=
p ← rows(S), s ← cols(S) − 1
·p−1 · s
¸2 ¸
P Q
s+1
R← 3p ·
(1 − 2 · Sk,j )
k=0
j=0
R
Отметим, что, несмотря на то, что MATHCAD15 является интерпретатором,
перечисленные выше программы работают достаточно эффективно и
требуют не более одной минуты времени на персональном компьютере
средней производительности при p = 5000000.
Н. К. Серегина
196
1. Оптимальные параллелепипедальные сетки
1.1. Оптимальные коэффициенты сетки. Покоординатный
алгоритм. Переходя к вопросу о вычислении оптимальных коэффициентов
по модулю p, прежде всего рассмотрим покоординатный алгоритм
вычисления, который требует для своей реализации O (s2 · p2 ) элементарных
арифметических операций:
HT 1(p, s) := as−1 ← 1, b0 ← 1, a0 ← 1
f or j ∈ 1..s − 1
R ← 3j+1
f or c ∈ 1..p − 1
bj ← c
"
#
· j h
³
´´i¸2
³
p−1
P Q
j+1
r← 3p
1 − 2 blp·k − f loor blp·k
k=0
aj ← c, R ← r
bj ← a j
l=0
if r < R
a
1.2. Обоснование
покоординатного
алгоритма. Необходимо
отметить, что первоначально Н. М. Коробов свой алгоритм обосновал
только для простых модулей. Позднее он его модифицировал для составных
модулей равных произведению двух простых, что позволило сократить
трудоемкость алгоритма вычисления оптимальных коэффициентов с O (N 2 )
до O (N ). Хотя в литературе не встречается обоснование этого алгоритма
для произвольного составного модуля, но не трудно понять, что такое
обоснование, но без оценок качества, нетрудно дать.
Действительно, из результатов разных авторов (например, см. [3])
следует существование оптимальных коэффициентов по любому составному
модулю. Кроме того, известно, что s-мерный набор оптимальных
коэффициентов всегда можно дополнить до (s + 1)-мерного набора. Поэтому,
если среди всех (s + 1)-мерных наборов с фиксированными первыми s
коэффициентами взять набор с минимальным значением функции H —
количественной меры качества оптимальных коэффициентов, то этот набор
и будет оптимальным.
1.3. Произведение сеток. Пусть ~x = (x1 , . . . , xs ) ∈ Rs , тогда {~x} =
= ({x1 }, . . . , {xs }) — дробная часть ~x.
Для любых двух сеток M1 и M2 их произведение есть сетка
M = M1 · M2 = {{~x + ~y }|~x ∈ M1 , ~y ∈ M2 }.
Различные аспекты теории произведения обобщенных параллепипедальных
сеток рассматриваются в работе [5].
При построении алгоритма Л. П. Добровольской поиска оптимальных
коэффициентов существенно используется представление параллелепипе-
Алгоритмы численного интегрирования с правилом остановки
197
дальной сетки по составному модулю как произведение сеток по взаимно
простым модулям.
1.4. Модифицированный алгоритм Л. П. Добровольской.
Следующая программа по заданной сетке S из P точек, представленной
двумерным массивом, находит оптимальную параллелепипедальную
сетку из p точек, такую, что их произведение будет оптимальной
параллелепипедальной сеткой из P · p точек. Важной особенностью
данной программной реализацией является тот факт, что массив с
результирующей сеткой представляет из себя p подмассивов, каждый из
которых суть модифицированная исходная параллелепипедальная сетка.
Эта особенность потом будет использована при реализации алгоритма
численного интегрирования с правилом остановки.
Данный алгоритм является модификацией алгоритма Л. П.
Добровольской.
Здесь вместо очень сложной системы целевых функций используется
количественная мера качества сетки H(S):
HT S(S, p, s) := as−1 ← 1, b0 ← 1, a0 ← 1, P ← rows(S)
f or j ∈ 1..s − 1
R ← 3j+1
f or c ∈ 1..p − 1
bj ← c
· P −1 p−1 · j h
³
P P Q
j+1
1 − 2 Sn,l +
r ← 3p·P
n=0 k=0 l=0
aj ← c, R ← r
bj ← aj
SP ·p−1,s−1 ← 0
f or k ∈ 1..p − 1
f or n ∈ 0..P − 1
f or j ∈ 0..s − 1
if r < R
SP ·k+n,j ← Sn,j +
S
aj ·k
p
bl ·k
−
p
³
−f loor Sn,l +
³
− f loor Sn,j +
aj ·k
p
bl ·k
p
´´i¸ 2
#
´
Программа HT S(S, p, s) требует для своего выполнения память объемом
O (s · P · p) и O (s2 · P · p2 ) элементарных арифметических операций.
1.5. Об алгоритме Л. П. Добровольской. Алгоритм Л. П. Добровольской основан на двух идеях.
Во-первых, этот алгоритм предназначен для поиска оптимальных
коэффициентов по составному модулю и параллелепипедальная сетка
рассматривается как произведение двух параллелепипедальных сеток по
взаимно простым модулям.
Во-вторых, хотя алгоритм работает для любого составного модуля, но
своих оптимальных характеристик он достигает, когда модуль является
произведением специальной последовательности простых чисел.
Н. К. Серегина
198
1.6. Допустимая последовательность простых. Последовательность p0 , . . . , pk называется допустимой последовательностью простых, если
для P = p0 · . . . · pk выполнены соотношения
pk · pk−1 · . . . · pj+1 · p2j 6
P
2j−2
(1 6 j 6 k).
Для этого достаточно, чтобы
pj−1 · . . . · p1 · p0
(1 6 j 6 k).
pj <
2j−2
Если
pj−1 · . . . · p1 · p0
pj−1 · . . . · p1 · p0
< pj <
(1 6 j 6 k),
j−1
2
2j−2
то получаем монотонную последовательность простых p0 < . . . < pk .
1.7. Программная реализация модифицированного алгоритма
Л. П. Добровольской. Следующая программа HT S1(pp, s) по заданному
вектору pp = (p0 , . . . , pk )T , где p0 , . . . , pk допустимая последовательность
простых, строит оптимальную параллелепипедальную сетку из P = p0 · . . . ·
pk точек:
HT S1(pp, s) :=
k ← length(pp) − 1
a ← HT 1(ppk , s), S ← Setka(ppk , a)
f or j ∈ k − 1, k − 2..0
S ← HT S(S, ppj , s)
S
1.8. Трудоемкость алгоритма. Нетрудно видеть, что количество
операций требуемых для выполнения этой программы для допустимой
последовательности простых есть
¡
¢
O s2 (p2k + pk · p2k−1 + . . . + pk · pk−1 · . . . · p1 · p20 ) =
!!
à à k
X P
¡
¢
+ P · p0
= O s2 · P · (4 + p0 ) .
= O s2
j−2
2
j=1
2. Алгоритм интегрирования с правилом остановки
2.1. Тестовая
функция
s+1
f (x) = 3
функция.
·
"
s
Y
(1 − 2 · xj )
j=0
которая реализуется программой
f (x) :=
s ← length(x) − 1
· s
¸2
Q
3s+1 ·
(1 − 2 · xj )
j=0
В качестве тестирующей предлагается
#2
,
Z
0
1
...
Z
0
1
f (x)dx0 . . . dxs = 1,
Алгоритмы численного интегрирования с правилом остановки
199
2.2. Мультипликативная дискретная дисперсия по произведению
параллелепипедальных сеток. Пусть M1 и M2 — две параллелепипедальные сетки и M = M1 · M2 — их произведение.
Мультипликативной
дискретной
дисперсией
по
произведению
параллелепипедальных сеток M для периодической функции f (x)
называется величина
Ã
!2
1 X
1 X
DM [f ] =
f ({x + y}) −
|M2 |
|M1 |
y∈M2
x∈M1
Ã
!2
1 X 1 X
−
f ({x + y}) .
|M2 |
|M1 |
y∈M2
x∈M1
2.3. Численное интегрирование с правилом остановки.
Следующая ниже программа реализует алгоритм численного интегрирования
с правилом остановки. Предполагается, что сетка S вычислена по алгоритму
Добровольской. Тогда пользуясь тем, что двумерный массив, задающий
сетку, специально организован, данная программа вычисляет приближенное
значение интеграла от функции f (x), наращивая количество точек от первой
сетки из pk точек, до последней сетки, если это будет необходимо, из P
точек.
В качестве правила остановки берется малость величины мультипликативной дискретной дисперсии:
IntP o(pp, S, f, ε) :=
s ← cols(S) − 1, k ← length(pp) − 1
x ← submatrix(S, 0, 0, 0, s)T , M ← f (x), m ← M, Int ← 0
f or n ∈ 0..ppk − 1
x ← submatrix(S, n, n, 0, s)T , r ← f (x)
M ← r if r > M
m ← r if r < m otherwise
Int ← Int + r
Int
Int ← pp
, D ← 2 · ε, j ← k − 1, P ← ppk
k
while (D > ε) · (j > 0)
D ← Int2
f or l ∈ 1..ppj − 1
Ir ← 0
f or n ∈ 0..ppk − 1
P n ← P · l + n, x ← submatrix(S, P n, P n, 0, s)T
r ← f (x)
M ← r if r > M
m ← r if r < m otherwise
Ir ← Ir + r
Ir ← Ir
, D ← D + Ir2 , Int ← Int + Ir
P
D
Int
, D ← pp
− Int2 , P ← P · ppj , j ← j − 1
Int ← pp
j
(Int
D
m
j
M)
Программа в качестве результата выдает приближенное значение
интеграла, величину мультипликативной дискретной дисперсии, а также
минимальное и максимальное значения функции на сетке.
200
Н. К. Серегина
Работа может быть использована в рамках учебного процесса при
изучении квадратурных формул в курсах вычислительной математики и
численных методов, а также в курсе по выбору, посвященном изучению
теоретико-числовых методов в приближенном анализе
Список литературы
1. Бабенко К.И. Основы численного анализа. М.: Наука, 1986.
2. Бахвалов Н.С. О приближенном вычислении кратных интегралов // Вестник
МГУ. Сер.1. Математика. 1959. №4. С.3–18.
3. Бочарова (Добровольская) Л.П. Алгоритмы поиска оптимальных коэффициентов // Чебышевский сборник. Тула: Изд-во ТГПУ им. Л. Н. Толстого. 2007. Т. 8.
Вып. 1(21). С. 4–109.
4. Добровольская Л.П., Добровольский Н.М., Симонов А.С. О погрешности
приближенного интегрирования по модифицированным сеткам // Чебышевский
сборник. Тула: Изд-во ТГПУ им. Л.Н. Толстого. 2008. Т. 9. Вып. 1(25).
С. 185–223.
5. Добровольский М.Н., Добровольский Н.М., Киселева О.В. О произведении
обобщенных параллелепипедальных сеток целочисленных решёток //
Чебышевский сборник Тула: Изд-во ТГПУ им. Л.Н. Толстого. 2002. Т. 3.
Вып. 2(4). С. 43–59.
6. Коробов Н.М. Вычисление кратных интегралов методом оптимальных
коэффициентов // Вестник МГУ. Сер. 1. Математика. 1959. № 4. С. 19–25.
7. Коробов Н.М. О приближенном вычислении кратных интегралов // ДАН СССР.
1959. Т. 124. № 6. С.1207–1210.
8. Коробов Н.М. Свойства и вычисление оптимальных коэффициентов // ДАН
СССР. 1960. Т. 132. № 5. С.1009–1012.
9. Коробов Н.М. Теоретико-числовые методы в приближенном анализе. М.:
Физматгиз, 1963.
10. Коробов Н.М. Теоретико-числовые методы в приближенном анализе. М.:
МЦНМО, 2004.
11. Локуциевский О.В., Гавриков М.Б. Начала численного анализа. М.: ТОО Янус,
1995.
12. Ребров Е.Д. Алгоритм Добровольской и численное интегрирование с правилом
остановки // Чебышевский сборник. Тула: Изд-во ТГПУ им. Л.Н. Толстого.
2009. Т. 10. Вып. 1(29). С. 65–77.
13. Огородничук Н.К, Ребров Е.Д. Об алгоритме численного интегрирования с
правилом остановки // Алгебра и теория чисел: современные проблемы и
приложения: матер. 7 Международной конф. Тула: Изд-во ТГПУ им. Л.Н.
Толстого, 2010. С. 153–158.
14. Огородничук Н.К., Ребров Е.Д. ПОИВС ТМК: Алгоритмы интегрирования
с правилом остановки // Многомасштабное моделирование структур
и нанотехнологии матер. Международной научно-практической конф.,
посвященной 190-летию со дня рождения академика Пафнутия Львовича
Чебышёва, 100-летию со дня рождения академика Сергея Васильевича
Вонсовского и 80-летию со дня рождения члена-корреспондента Виктора
Алгоритмы численного интегрирования с правилом остановки
201
Анатольевича Буравихина, Тула: Изд-во ТГПУ им. Л.Н. Толстого, 2011.
С. 153–158.
15. Algorithms for computing optimal coefficients / N.M. Dobrovolskiy [et al.] // Computer Algebra and Information Technology: book of abstracts of the International
scientific conf. Odessa, August 20–26. 2012. P. 22–24.
16. Некоторые вопросы теоретико-числового метода в приближенном анализе /
Л.П. Добровольская [и др.] // Алгебра и теория чисел: современные проблемы
и приложения: тр. X Международной конф. Ученые записки Орловского
государственного университета. 2012. № 6. Ч. 2. С. 90–98.
Серегина Надежда Константиновна (nadj_nadj@mail.ru), аспирант,
кафедра алгебры, математического анализа и геометрии, Тульский
государственный педагогический университет имени Л.Н. Толстого, Тула.
Algorithms for numerical integration with stopping rule
N. K. Seregina
Abstract. The paper deals with algorithms for numerical integration of the
stopping rule for Parallelepipedal nets. We present a program in MATHCAD15
that implement these algorithms.
Keywords: parallelepipedal net, algorithm with a stopping rule, quadrature
formula, method of optimal coefficients, quantitative measure of the quality of
the grid.
Seregina Nadegda (nadj_nadj@mail.ru), postgraduate ctudent, department of
mathematical analysis, algebra and geometry, Leo Tolstoy Tula State Pedagogical
University.
Поступила 23.08.2013
Документ
Категория
Без категории
Просмотров
4
Размер файла
160 Кб
Теги
алгоритм, остановке, интегрированный, правилом, численного
1/--страниц
Пожаловаться на содержимое документа