CC BY
9
Вычислительные технологии
Том 22, № 1, 2017
О выборе плотностей распределения узлов адаптивных сеток в стохастических алгоритмах многомерного численного интегрирования и приближения функций
А. В. ВойтишЕК*, Д. А. Прасол
Институт вычислительной математики и математической геофизики СО РАН,
Новосибирск, Россия
*Контактный e-mail: vav@osmf.sscc.ru
Сформулирован аналог принципа эквираспределения для стохастических адаптивных сеток, используемых при решении многомерных задач численного интегрирования. Рассмотрена проблема выбора плотности распределения узлов при реализации алгоритма самоорганизации Т. Кохонена для построения адаптивных сеток при решении задачи численного приближения функций. Для гладких приближаемых функций в одномерном случае предложено использовать в качестве плотности адаптивных узлов корень квадратный из модуля второй производной приближаемой функции. Приведен пример функции, для которой данная рекомендация приводит к улучшению аппроксимационных свойств сеточного приближения.
Ключевые слова: адаптивная сетка, принцип эквираспределения, вероятностная плотность распределения, метод выборки по важности, схема Т. Кохонена, выбор плотности узлов при приближении функций.
Введение
При численном решении ряда прикладных задач целесообразно использовать адаптивные сетки, в которых сгущения узлов в подобластях позволяют более эффективно получать приближенные решения.
В данной работе обсуждаются применение и численная реализация адаптивных сеток в алгоритмах численного интегрирования и аппроксимации функций. Особое внимание уделено решению многомерных задач, в которых целесообразно применять стохастические (вероятностные) алгоритмы построения адаптивных сеток.
В разд. 1 приведены рассуждения о том, что в основном принципе построения адаптивных сеток — принципе эквираспределения [1] — так называемая управляющая функция имеет основные свойства вероятностной плотности распределения [2]. Сформулирован вероятностный аналог принципа эквираспределения (в интегральной форме) для случая, когда управляющая функция в точности является вероятностной плотностью, индуцирующей неупорядоченный адаптивный набор точек (как в алгоритме выборки по важности [3]).
© ИВТ СО РАН, 2017
В разд. 2 приведен краткий анализ использования стохастических адаптивных сеток при построении квадратурных и кубатурных формул для приближенного вычисления интегралов [3, 4], в том числе рассмотрены особенности компьютерной реализации упомянутого выше алгоритма выборки по важности.
В разд. 3 приведены соображения об использовании стохастического алгоритма построения адаптивных сеток — схемы Т. Кохонена [5, 6] — в задаче приближения функций, а также сформулированы рекомендации по выбору плотности распределения "узлов притяжения" в этом алгоритме.
Отметим особенность представления материала в данной работе. Хотя речь идет об алгоритмах, особенно эффективных для многомерных задач, ряд принципиальных аналитических рассуждений и рекомендаций в тексте сделан для одномерного случая. Перенос соответствующих суждений на многомерный случай технически не сложен (соответствующие примеры приведены в заключении данной статьи), однако он требует использования "тяжеловесных" многомерных обозначений и рассуждений, что приводит к затруднениям при восприятии основных идей данной статьи.
1. Вероятностный аналог принципа эквираспределения
При построении адаптивных сеток, как правило, реализуется тот или иной вариант принципа эквираспределения [1]. Этот принцип связан с построением невырожденного отображения x = x(y), которое переводит зафиксированную (чаще всего равномерную) сетку на "вычислительной" области Q в адаптивную сетку на заданной "физической" области X, и удовлетворяет соотношению
/ (x(y ))J (y) = Hi = const, (1.1)
где J(y) — якобиан отображения; Hi — положительная константа; f (x) — специально выбранная управляющая функция, такая, что
/ (x) > 0 при x е X и / (x) = 0 при x/X. (1.2)
На практике обычно рассматриваются случаи, когда X С Rd, а Q С Rs, где 1 < s < d.
В дальнейшем будем полагать, что функция f (x) непрерывна или кусочно-непрерывна, и выполнено соотношение
Jf (x) dx =1. (1.3)
Свойствами (1.2), (1.3) обладает, в частности, вероятностная плотность распределения (x) некоторого случайного вектора (d-мерной случайной величины) ^ [2]. Но более существенным свойством вероятностной плотности является то, что она определяет "сгущения" выборочных значений (реализаций) ^ случайного вектора
Управляющая функция f (x), обладающая свойствами (1.1) — (1.3), также определяет "сгущения" узлов реализуемой с ее помощью адаптивной сетки.
Покажем это для простейшего одномерного случая, когда d =1 и
x = Ж е X с R и y = у е Q =[0,1] с R.
< Хм-I < Хм}
(1.4)
где
х
3 = х(уj), у, = jAy, Ay = 1/М, j = 0,1,..., М, (1.5)
в этом случае связано с рассмотрением частного случая соотношения (1.1)
f (х(у))х'у (у) = Н2 = const. (L6)
Здесь х'у (у) обозначает соответствующую производную по переменной у (это в данном случае якобиан отображения х = х(у)).
Интегрируя это соотношение на отрезке [yj ,yj+i] и предполагая непрерывность функции f (х) на этом отрезке, получаем
xj+1
i f (z) dz = H2(yJ+1 - y3) (1.7)
или
f (Хз)(ж,-+1 - х,) = H3, (1.8)
Н2
где Ху е ух,,Хз+1\, Н3 = Н2Ау = —, ] = 0,1,...,М - 1.
Таким образом, при заданном количестве узлов М для подобластей, в которых функция /(ж) велика, длины отрезков [х^ ,Хз+\\ малы.
Учитывая соотношения (1.2), (1.3), (1.8), функцию /(х) можно назвать плотностью распределения узлов адаптивной сетки.
Теперь рассмотрим одномерную неупорядоченную стохастическую адаптивную сетку для метода Монте-Карло (см. [3], а также разд. 2 данной работы)
= {^1,...,^»} ,
где ^г — независимые выборочные значения случайной величины £, имеющей плотность распределения Д(г), а также соответствующий аналог левой части соотношения (1.7) (которое можно назвать интегральной формой принципа эквираспределения):
€¿+1
1(г) = ( Ш ¿г = ^(6+1) - Р^г), г = 1,...,п - 1.
Здесь ^(ж) = Р{£ < х} = Д(г) ¿х — функция распределения случайной величины
и —<х
£ [2].
Хорошо известно [3], что случайная величина
« = Рц (О
имеет равномерное распределение в интервале (0,1). Таким образом, случайная величина I(г) равна
€¿+1
[ & (г) йг = а(1) - а(2) (1.9)
X
для любой плотности распределения Д (г) и любых случайных узлов ^ и £¿+1- Здесь а(1) и а(2) — независимые стандартные случайные величины [3], равномерно распределенные в интервале (0,1).
Используем следующее утверждение.
Лемма [2]. О'умма ^(1) + ^(2) двух независимых случайных величин ^(1) и ^(2) с плотностями распределения /(1)(^) и /(2)(^) имеет плотность распределения в виде свертки
(г) = (/^(1) * /ч(2)) (г) = J /(1)(м)/(2)(г - м) ¿и.
Из леммы следует, что случайная величина I(г) имеет "универсальное" (для различных г = 1, ...,п — 1) распределение с плотностью
(1 + г для — 1 < г < 0, 1 — г для 0 < г < 1, (1.10)
0 иначе.
Здесь р(1)(г) — В -сплайн первого порядка (или "функция-крышка" [7]).
Таким образом, соотношение (1.9) можно назвать вероятностным аналогом принципа эквираспределения.
2. Использование адаптивных сеток в стохастических алгоритмах численного интегрирования
Рассмотрим задачу приближенного вычисления многократного интеграла
I = J д(х) ¿х = J д(х) ¿х = J д(х) ¿х.
X
Известно [3, 4], что для достаточно больших размерностей й (конкретнее, для d > 4 [3]) наиболее эффективным алгоритмом вычисления интеграла I является весовой метод Монте-Карло, основанный на представлении интеграла в виде математического ожидания
'=/*=<=т
где случайный вектор ^ распределен согласно выбираемой плотности (х) с последующим численным приближением интеграла I в соответствии с законом больших чисел [2, 3]
1 , 9(Ь)
п
1 - Я < = ^
г=1
Здесь ^1,...,— реализуемые на компьютере выборочные значения случайного вектора
Одно из принципиальных преимуществ метода (2.1) состоит в том, что он допускает корректную аналитическую оптимизацию на основе принципа "тот алгоритм эффективнее, который позволяет получить заданный уровень погрешности за меньшее
время счета", в то время как "детерминированные" методы численного интегрирования — квадратурные и кубатурные формулы — сравниваются (не вполне корректно) только по зависимости погрешности от шага сетки [4].
Затраты (время счета) алгоритма (2.1) равны величине
s = nt, (2.2)
где t — среднее время для получения выборочного значения Q = g(^i)/f^(£j). Как несложно показать [3], при фиксированном уровне погрешности величина (2.2) пропорциональна величине трудоемкости
S = DC ■ t, (2.3)
где D^ — дисперсия случайной величины ("весовой оценки") (.
В свою очередь, при оптимизации алгоритма (2.1) на основании минимизации величины (2.3) широко используется принцип выборки по важности, состоящий в том, что среди алгоритмов (2.1), имеющих примерно одинаковое среднее время t реализации выборочных значений Q, тот имеет меньшую трудоемкость S из (2.3), для которого
(x) « #4 |0(x)| , Н4 = (U)| du. (2.4)
Термин "выборка по важности" соответствует английскому термину importance sampling. Такое название объясняется тем, что если (x) пропорциональна модулю подынтегральной функции g(x), то в частях области X, где |#(x)| больше и вклад которых в интеграл I более существен, будет выбираться больше случайных точек
Таким образом, в алгоритме выборки по важности (2.1), (2.4) реализуется адаптивная стохастическая сетка
= {&,..., U . (2.5)
Что касается конкретных алгоритмов компьютерной реализации выборки (2.5), то здесь можно воспользоваться соображениями из разд. 3.2.4 учебника [3] (более детально они представлены в [8]). Эти соображения состоят в том, что в качестве плотности (x) целесообразно использовать приближения модуля подынтегральной функции g(x) вида
м
f£(x) = H5Lm|y(x)| = #5 ^ Wj(g) Xj(x), (2.6)
j=1
где {Xj(x)} — заданные базисные функции, согласованные с заданными узловыми точками {xi,..., x м}, коэффициенты {wj(g)} зависят от значений g = (|g(xi)|,... , |^(xM)|), а Н5 — нормирующая константа.
Здесь возникают требования моделируемости приближения Lm|#(x)| [3, 8]. Из соотношения (2.6) следует, что для численного моделирования адаптивных случайных узлов (2.5) можно использовать алгоритм метода суперпозиции [3].
Алгоритм 1.
1. Используя один из алгоритмов моделирования целочисленной случайной величины rq (см., например, [3]) с распределением
P{V = j} = Pj, Pj = ^g* #6j) = /ь»du, i = 1,...,M,
щл J
реализуем выборочное значение щ = т, г =1,..., п.
2. Реализуем выборочное значение ^ согласно плотности
/т(х) = н6т)Хт(х), х = (х(1),...,х(а)) е И• (2.7)
Для корректной реализации алгоритма 1 требуется неотрицательность коэффициентов {wj^)}. Кроме того, базисные функции (х)} должны быть пропорциональными плотностям распределения случайных векторов вида (2.7), для которых имеются эффективные алгоритмы реализации выборочных значений. Последние два требования определяют моделируемость приближения Ьм|#(х)| из соотношения (2.6).
Проведенные в работе [8] исследования известных приближений вида (2.6) показали, что требованиям моделируемости наилучшим образом удовлетворяет конечно-элементная аппроксимация Стренга — Фикса [7]. Базисные функции строятся следующим образом.
Предполагается, что в И задана равномерная прямоугольная сетка с шагом Н и каждый узел х^, используемый в аппроксимации (2.6), можно описать как
х.- , Л1)
где к^) — целые числа. Соответствующая узлу х^ базисная функция ^ (х) имеет специальный вид
* (х) = е (г){ х - " м( х - • (2.8)
Здесь [3(г)(г) — Б-сплайн г-го порядка.
Для базиса Стренга — Фикса (2.8) плотность (2.7), используемая в алгоритме 1, имеет вид
/т(х) = /т^ )) • ••• • /т^ , = Н^^ ^~Н , в =
т. е. в данном случае Н6т) = НЛ. Это означает, что при моделировании вектора ^ = ^1), следует учесть, что компоненты независимы и моделируются по фор-
муле
^ = Н4Г) + к$)Н, (2.9)
где — выборочное значение случайной величины ф(г), имеющей плотность распределения [3(г) (г).
Учитывая рекуррентное определение 5-сплайна г-го порядка
/3(г)(г)= (/3(г-1) * $(0)) (г),
где
р (0)м = {
Н 0 Для " 1/2 " 2 - 1/2' (2.Ю)
1 0 иначе,
а также сформулированную выше лемму и то обстоятельство, что функция [3(0)(г) из (2.10) является плотностью распределения случайной величины а(0) — 1/2 (а(0) — стандартная случайная величина), несложно обосновать моделирующую формулу
Ф() = ^«о - Т^ + уа1 - + ... + ^ - Т^ = ао + «1 + ... + аг - , (2.11)
где а0,а1, ...,аг € и(0,1) — стандартные случайные числа [3].
Что касается практического применения алгоритма 1 с формулами (2.9), (2.11), то случай г = 0 (кусочно-постоянное приближение (2.6)) подробно рассмотрен в работе [9].
Трудности использования "гладких" приближений (2.6), (2.8) (для достоточно больших г) описаны в диссертации [10], где подробно рассмотрен случай кубических сплайнов для г = 3. Они связаны, в частности, с вычислением коэффициентов {wjдающих высокий порядок аппроксимации Ьм |#(х)| (здесь может нарушаться требование положительности этих коэффициентов).
В работе [8] показано, что наиболее эффективно применение так называемого муль-тилинейного приближения (2.6) с базисными функциями (2.8) и сплайнами первого порядка (1.10). Особо отметим, что в практических вычислениях вместо формулы (2.11), которая для случая г =1 имеет вид
фо = «1 + «2 - 1,
целесообразно использовать обоснованную в разд. 1 более экономичную формулу
Фо = «1 - «2 (2.12)
(см. соотношение (1.10)).
3. Стохастический алгоритм самообучения (схема Т. Кохонена) для построения адаптивных сеток в задаче приближения функций
Рассмотрим задачу численной аппроксимации функции ^(х) на некоторой компактной области X из И. Соответствующее численное приближение описывается формулой, аналогичной (2.6):
м
^(х) и Ьм^(х) = ^ Wj(х). (3.1)
j=l
Здесь, как и ранее, {xj (х)} — заданные базисные функции, согласованные с узлами сетки
X(м) = {х1,..., хм}, (3.2)
а коэффициенты {wj(^)} зависят от значений ^ = (^(х1),... , ^(хм)).
Рассмотрим вопросы построения адаптивной сетки (3.2) на основе принципа экви-распределения (1.1) для задачи (3.1).
Для простоты рассмотрим сначала одномерный случай (д = 1), когда Q = X = [0,1]. Здесь для построения приближений точек адаптивной сетки (3.2) (или (1.4)) можно применить следующую численную схему [1].
Дифференцируя соотношение (1.6) по у, получаем уравнение
{
(/(х(у))х'у(у)Уу = ° (3 3)
ж(0) = 0, ж(1) = 1. (3.3)
Аппроксимация этой задачи на равномерной сетке
{у3 = з/М; j = 0,1,..., М} (3.4)
в вычислительной области Q = [0,1] (см. соотношение (1.5)) приводит к системе уравнений для Xj = x(yj):
fj+1/2xj+1 - (fj+1/2 + fj-1/2)xj + Л-1/2Х-1 = 0, j = 1 М - (3.5)
W = = f (X(V')+2x(V'+1)) , Л-1/2 = f(x1-1/2) = f (х(v'-1)+x(у')) ,
которую с учетом граничных условий и того факта, что значения функции f(x) заданы во всех точках X = [0,1], можно решить методом прогонки [1].
Несложно показать, что для решений {xj} системы (3.5) выполнен принцип экви-распределения в разностной форме (это аналог соотношения (1.8)):
W-W = f (З^) ^^^ = = COnSt. (3.6)
Описанная методика построения приближений точек (1.4), (3.2) весьма чувствительна к размерности задачи. Когда эта размерность равна двум и более, построение соотношений типа (3.3)-(3.6) сопряжено с выполнением специальных критериев качества сетки, связанных с требованиями выпуклости и вытянутости ячеек, ортогональности линий сетки и т. д. Это обстоятельство вместе с упомянутой выше проблемой выбора управляющей функции /(x) представляет собой трудно разрешимую проблему даже для приближения решений ^(x) простейших классических задач математической физики [1].
Многих указанных недостатков лишен стохастический алгоритм построения адаптивных сеток — схема Т. Кохонена [5, 6]. Его главным преимуществом является свойство самообучения, обеспечивающее автоматизированное построение адаптивных сеток с нужными свойствами.
Кратко опишем стохастический алгоритм в общем (многомерном) случае [5, 6]. Пусть в "физической" области X (или на ее поверхности Gx) требуется построить сетку (3.2), распределение узлов которой соответствует заданной плотности /(x), x Е Rd. Структура этой сетки (порядок и структура расположения узлов) задается "картой" Q(M) = {q1,..., Ям} и системой так называемых нейронов Е(м) = {e1,..., вм} (где ег = (qi, x^)), определяющих соответствие между сетками X(м) и Q(M).
Приближение нейронов осуществляется с помощью процедуры самообучения, представляющей собой итерационный процесс, основанный на последовательном формировании обучающего множества S(T) = {^0(1),...,^0(Т)} в виде выборки из вероятностного распределения случайного вектора имеющего плотность /(x). Здесь Т — количество итераций, а £0(t) Е X. Кроме того, с помощью значений специального коэффициента обучения d4j (t, qj) Е [0,1) на каждом шаге итерации устанавливаются латеральные связи между нейронами е^ и ej. В результате получается последовательность сеток X(м\t) = {x1(i),... , xm(t)}, при этом требуется выполнение соотношения
X(м) « X(м) = lim X(м\t) « X(м)(Т).
(3.7)
Знаки приближенных равенств в соотношении (3.7) означают не только малость расстояний p(xj(то),x), p(x(го), Xj(Т)) и p(x(Т),x), где
p(x, y) = - у(1))2 + ... + (ж№ - y(d))2, x = (x(1),...,x(d)), y = (y(1),...,y(d)),
но и воспроизведение (с хорошей точностью) требуемых свойств сетки X(м) (например, свойств прямоугольности, отсутствия граничного эффекта и др.) при реализации следующего итерационного процесса. Алгоритм 2 [5, 6].
1. Устанавливаются начальные положения узлов сетки X(м)(0) = {x1(0),..., хм (0)}.
2. На каждой итерации с номером t = 1,... ,Т выполняются следующие действия:
а) выбирается очередной элемент £0(i) выборки S(T);
б) вычисляются расстояния p(£0(i), x(t — 1)) от точки £0(i) до всех узлов x(t — 1) и выбирается ближайший к (t) узел xm(t — 1) в соответствии с условием
т = arg min р ($o(t), x(t — 1));
г=1,...,М
такой узел xm(t — 1) называют победителем;
в) проводится корректировка положений всех узлов в соответствии с формулой
x(t) = x(t — 1) + вЧт (t, q) (£o(i) — x(t — 1)) (3.8)
для всех г = 1,..., M.
На каждой итерации (3.8) основного алгоритма узлы сетки сдвигаются к случайной точке ^0(i). Поэтому в местах высокой концентрации элементов выборки постепенно скапливается все больше узлов, за счет чего достигаются сгущения сетки. В работе [6] показано, что при Т ^ то и специальном выборе коэффициента обучения 84j (t, q^) основной алгоритм ведет к выполнению принципа эквираспределения, т. е. к получению требуемой плотности сетки, определяемой функцией f (x) (см. соотношения (1.1), (1.7)). Отметим также, что в работах [11, 12] предложены конструктивные подходы к аналитическому исследованию алгоритма 2.
Интересен и принципиален вопрос о выборе плотности f (x). Здесь определенные трудности возникают даже в одномерном случае. Например, в пособии [1] без особой аргументации для достаточно ограниченного класса задач предлагается выбирать управляющую функцию в виде
f (х) = 1 + С0Иж)|60 + W (ж)|61 + С2 (x)lb2,
где неотрицательные константы С к ,Ьк, к = 0,1, 2, подбираются экспериментально, а вместо производных берутся их сеточные приближения.
Авторами статьи предлагается следующий конструктивный подход к выбору плотности f (х). Рассмотрим кусочно-линейное приближение на равномерной сетке {yi} из соотношения (3.4):
Ьм Ф) = ipfa) + (х — у,) ^(Уз+1) — ^), j = 0,1,...,М — 1. (3.9)
У3+1 — Уз
Используя рассуждения из параграфа 3 гл. 2 монографии [4], можно показать, что для этого приближения выполнено соотношение
Ро[0,1](^,Ьм= sup 1р(х) — Ьмр(х)1 < 1 sup |<р"(x)lh2, (3.10)
же[о,1] 8 же[о,1]
1 = Vj+1 - Уз = 1/М, j = 0,1,...,М - 1.
Введем допустимый уровень погрешности и приравняем его правой части неравенства (3.10):
- sup |(^"(x)|12 = 8 или I sup '(x)| ■ 1 = Н8, (3.11)
8 же[о,1] у хе[о,1]
где Н8 = л/88. Сравнительный анализ соотношений (1.8) и (3.11) приводит к рекомендации использовать в качестве плотности ( x) сеточное приближение функции
Пх) = Н9у/ К(х)\, (3.12)
где Нд — соответствующая нормировочная константа.
Приведем пример, показывающий целесообразность выбора плотности в соответствии с рекомендацией (3.12). Рассмотрим функцию
р(х) = е-2х, х е [0, 2].
Введем равномерную сетку х^пг^отга^ = (^ — 1)Д, ] = 1,.., 9, с шагом А = 0.25. Затем, следуя рекомендации (3.12), определим плотность /(х) в виде
е2
/(х) = х-г, х е [0, 2]. (3.13)
е 2 — 1
Здесь использовано обстоятельство, что <р"(х) = 4е-2ж.
Заметим, что для плотности (3.13) можно построить моделирующую формулу метода обратной функции распределения [3] для соответствующего выборочного значения
е 2 - 1
ее = — ь(1 —е
(1 - V1 ао),
где а0 Е U(0,1) — стандартное случайное число. Эта формула может быть применена при реализации соответствующей одномерной версии алгоритма 2.
Применим принцип эквираспределения и рассмотрим соответствующее соотношение (1.6) с граничными условиями
f(x(y))x' (у) = Ню, Н10 = const,
x(0) = 0, x(1) = 2.
Решаем получившееся дифференциальное уравнение и находим x(y). Общее решение имеет вид
/ р2-х \
dx = Н10 dy
,е 2 - 1 или
-е2-ж = Ню(е2 - 1)у + Нп, Нп = const. Определяем константы из граничных условий
Н10 = 1, Нц = -е2.
Окончательно получаем
х(у) = 2 - 1п (е2 - (е2 - 1) у) , у е [0,1]. Результаты пересчета узлов адаптивной сетки
х((] - 1)к), э = 1,...,9, к = 0.125,
(adaptive)
X
г и 1 о (uni form) „ (adaptive)
запишем в табличном форме. Здесь xj — узлы равномерной сетки, xj —
узлы адаптивной сетки.
Далее построим кусочно-линейное приближение функции <р(х) = е-2х на равномерной и адаптивной сетках. Для этого вычислим значения этой функции в узлах равномерной и адаптивной сеток.
(wni form) X.j 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
(adaptiv e) Xj 0 0.11 0.24 0.39 0.57 0.78 1.05 1.41 2
2x(uniform) e i 1 0.6 0.368 0.223 0.135 0.082 0.05 0.03 0.018
(adaptive) e-2xi 1 0.8 0.618 0.458 0.320 0.210 0.122 0.06 0.018
Найдем погрешности
6(adaptive) = Pc[0,2] UbÎdaptive)<р) = max ^(х) - L^daptive) ф)
gWonn) = fL(umform)\ = max - L(rform)/Ax)l
\ / x€[0,2]
T (adaptive) / \ m r\\ r (adaptive)
где L*~M ф(х) — кусочно-линеиная аппроксимация (3.9) с узлами {xj },
т (uni form) / \ m r\\ r (uniform) л
a LM y(x) — кусочно-линеиная аппроксимация (3.9) с узлами {xj }.
На равномерной сетке максимальное отклонение §(adapUvе) достигается на интервале [0, 0.25] и равно 0.024. Соответствующие вычисления на адаптивной сетке дают максимальное отклонение 5(umform) на интервале [0.57, 0.78], равное 0.01. Получаем, что погрешность аппроксимации на адаптивной сетке заметно меньше.
Отметим, что если функция <р(х) не задана явно, а лишь известны ее значения в узлах сетки, то для вычисления второй производной в узлах сетки можно использовать ее сеточные приближения (по аналогии с работой [1]).
Заключение
Перечислим основные результаты данной работы.
1. Получен вероятностный аналог принципа эквираспределения (в интегральной форме, одномерный случай) — см. соотношения (1.9), (1.10).
2. Показано, что метод выборки по важности, широко используемый в задачах многомерного численного стохастического интегрирования, связан с формированием сетки (2.5), адаптивной в том смысле, что плотность распределения узлов пропорциональна модулю подынтегральной функции. Приведен также конструктивный алгоритм метода
суперпозиции (алгоритм 1) с конкретным аппроксимационным базисом (2.8) и эффективными (экономичными) моделирующими формулами (2.9), (2.12).
3. Для задачи приближения функций многих переменных рассмотрен стохастический итерационный алгоритм построения адаптивных сеток (схема Т. Кохонена — алгоритм 2). Предложен подход к выбору плотности распределения адаптивных узлов (одномерный случай — соотношение (3.12)) и на простом примере показана эффективность такого подхода.
Отметим следующие возможности развития и применения перечисленных результатов.
1. Относительно несложно получить многомерные аналоги полученных здесь результатов (например, соотношений (1.9), (1.10) и (3.12)). В частности, индукцией по размерности , используя разложение совместной плотности распределения случайного вектора в произведение условных плотностей [2, 3], можно получить следующий аналог фор -мулы (1.9) для одинаково распределенных (согласно плотности /(х) = / ( г(1),..., )
(-мерных случайных векторов = (£(1), ) и = (£(1),..., :
г(1)
?2
/(
Д1)
(г(1)..А г(л) = а(1)
а
(а) _ а^+1)
а
(2Л)
(4.1)
М) <1
где а(в), 5 = 1,..., 2(1, — независимые стандартные (равномерно распределенные в интервале (0,1)) случайные величины.
Заметим, что формула (4.1) и ее обоснование подтверждают высказанный во введении статьи тезис о том, что многомерные аналоги полученных в данной работе результатов и их обоснований имеют относительно простой, но достаточно громоздкий вид. Такой вид будет иметь и многомерный аналог формулы (1.10). Здесь для получения плотности распределения случайной величины из правой части соотношения (4.1) следует использовать известные формулы для произведения и суммы независимых случайных величин с известными плотностями распределения [2].
2. Для развития теории и приложений конструктивных алгоритмов выборки по важности целесообразны дополнительные исследования моделируемых базисов, отличных от конечных элементов (2.8). Особую роль здесь играет построение экономичных реализаций соответствующих алгоритмов с целью их включения в геометрические вычислительные модули (подобно тому, как это сделано в работе [9]).
3. Возможно обобщение подхода к выбору плотности узлов в задаче аппроксимации функций из разд. 3 на многомерный случай и (или) на случай использования более гладких аппроксимаций Ьм'-р(х) (см. соотношение (3.1)). Здесь можно предложить (по аналогии с формулой (3.12)) конструирование управляющей функции (плотности) /(х) в виде комбинации модулей производных приближаемой функции, входящих в выражения для оценок погрешностей рВ(х)(<р^м<р) = ||^(х) — Ьм<р(х)\\в(х) в норме соответствующего функционального пространства В(Х) [3, 4, 7, 10]. Важно также приумножать примеры эффективного использования многомерных адаптивных сеток при решении прикладных задач (подобно тому, как это сделано в пособии [1]).
Список литературы / References
[1] Хакимзянов Г.С., ШЪкин Ю.И. Разностные схемы на адаптивных сетках. Ч. 1,2. Новосибирск: НГУ, 2009.
Khakimzyanov, G.S., Shokin, Yu.I. Difference schemes on adaptive meshes. Vol. 1,2. Novosibirsk: NGU, 2009 (In Russ.)
[2] Боровков А.А. Теория вероятностей. М.: Наука, 1986. 432 с. Borovkov, A.A. Probability theory. Moscow: Nauka, 1986. 432 p. (In Russ.)
[3] Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло. М.: Изд. центр "Академия", 2006. 368 с.
Mikhaylov, G.A., Voytishek, A.V. Numerical statistical modelling. Monte Carlo methods. Moscow: Izd. Tsentr "Akademiya", 2006. 368 p. (In Russ.)
[4] Березин И.С., Жидков Н.П. Методы вычислений. Т. 1. М.: Гос. изд-во физ.-мат. лит. 1959. 464 с.
Berezin, I.S., Zhidkov, N.P. Numerical methods. Vol. 1. Moscow: Gos. izd-vo fiz.-mat. lit., 1959. 464 p. (In Russ.)
[5] Kohonen, T. Self-organizing maps. Springer-Verlag, 2001. 502 p.
[6] Нечаева О.И. Нейросетевые модели, алгоритмы и комплекс программ для построения адаптивных сеток: Дис. ... канд. физ.-мат. наук. Новосибирск: НГУ, 2007. 120 с. Nechaeva, O.I. Neural network models, algorithms and program complex for constructing adaptive meshes: Dis. ... kand. fiz.-mat. nauk. Novosibirsk: NGU, 2007. 120 p. (In Russ.)
[7] Марчук Г.И., Агошков В.И. Введение в проекционно-сеточные методы. М.: Наука, 1981. 416 с.
Marchuk, G.I., Agoshkov, V.I. Introduction in projective-mesh methods. Moscow: Nauka, 1981. 416 p. (In Russ.)
[8] Voytishek, A.V., Kablukova, E.G. Using the approximation functional bases in Monte Carlo methods // Russ. J. of Numer. Analysis and Mathem. Modelling. 2003. Vol. 18, No. 6. P. 521-542.
[9] Бессмельцев М.В., Войтишек А.В. Модификации геометрических программных модулей, связанные с построением моделируемых вероятностных плотностей // Вычисл. технологии. 2011. Т. 16, № 4. С. 19-36.
Bessmeltsev, M.V., Voytishek, A.V. Modification of geometrical program modules, tied with construction of modelled probabilistic densities // Comput. Technologies. 2011. Vol. 16, No. 4. P. 19-36. (In Russ.)
[10] Милосердов В.В. Дискретно-стохастические численные алгоритмы со сплайн-восполнениями: Дис. ... канд. физ.-матем. наук. Новосибирск: НГУ, 2006. 83 с. Miloserdov, V.V. Discrete-stochastic numerical algorithms with spline approximations: Dis. ... kand. fiz.-mat. nauk. Novosibirsk: NGU, 2007. 83 p. (In Russ.)
[11] Войтишек А.В., Хмель Д.С. Аналитическое описание применения одномерной схемы Т. Кохонена для построения адаптивных сеток // Сиб. журн. вычисл. математики. 2011. Т. 14, № 2. С. 131-140.
Voytishek, A.V., Khmel, D.S. Analytical description for application of the 1D Kohonen scheme for constructing adaptive grids // Numerical Analysis and Applications. 2011. Vol. 4, No. 2. P. 105-113.
[12] Войтишек А.В., Хмель Д.С. Аналитический подход к изучению граничного эффекта в одномерном рандомизированном численном алгоритме построения адаптивных сеток // Журн. вычисл. математики и матем. физики. 2013. Т. 53, № 2. C. 195-208. Voytishek, A.V., Khmel, D.S. Analytical approach to study the boundary effect in the 1D randomized numerical algorithm for constructing adaptive grids // Zurn. Vychisl. Matematiki i Matem. Fiziki. 2013. Vol. 53, No. 2. P. 195-208. (In Russ.)
Поступила в 'редакцию 24 августа 2016 г., с доработки —13 января 2017 г.
On the choice of distribution densities in adaptive mesh nodes for stochastic algorithms of numerical integration and function approximation
Voytishek, Anton V.*, Prasol, Denis A.
Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk, 630090, Russia
* Corresponding author: Voytishek, Anton V., e-mail: vav@osmf.sscc.ru
In this paper we formulate the analog of the equidistributional principle for stochastic adaptive meshes which are used for solving the multi-dimensional problems of numerical integration. It means that the integral for density of continuous 1D random variable £ over the interval (£i, £2) between two sample values of the variable £ has the "universal" distribution with the density equal to the so called "hat-function".
The importance sampling technique is treated as a method for constructing an adaptive mesh for the problem of numerical integration. The discrete-stochastic approach with the "modelled" basis of the Strang — Fix approximation for numerical simulation of the corresponding adaptive mesh is proposed.
We also consider the problem on the choice of the distribution density of nodes in implementation of self-organization in the T. Kohonen's algorithm. The research is aimed at construction of adaptive meshes for solving the problem of numerical function approximation. For the smooth approximated functions we propose to use the square root of the second derivative module of the approximated function as density of adaptive nodes. The example of function for which the formulated recommendation leads to improvement of approximating properties of mesh approach is given.
The most part of results is formulated in the simplest ID-version. The survey of possible multi-dimensional versions of these results and directions of further investigations is presented in the conclusion.
Keywords: adaptive mesh, equidistributional principle, probabilistic distribution density, Monte Carlo method, importance sampling method, T. Kohonen's scheme, choice of nodes density for function approximation.
Received 24 August 2016 Received in revised form 13 January 2017
© ICT SB RAS, 2017
