Сравнение различных вариантов рандомизации метода последовательных приближений Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 11, Специальный выпуск, 2006

СРАВНЕНИЕ РАЗЛИЧНЫХ ВАРИАНТОВ РАНДОМИЗАЦИИ МЕТОДА ПОСЛЕДОВАТЕЛЬНЫХ ПРИБЛИЖЕНИЙ*

А. В. Войтишек, Е.Г. Каблукова Институт вычислительной математики и математической геофизики СО РАН, Новосибирск, Россия

e-mail: vav@osmf. sscc. ru

О. С. Герасимова Новосибирский государственный университет, Россия e-mail: gos@gorodok.net

Numerical algorithms for approximation of a solution of an integral equation of the second kind are investigated. Randomization of finite and infinite segments of Newmann series is used. Possibilities of exploitation of the homogeneous Markov chain segments of finite and random length or effective discrete-stochastic methods of numerical integration are examined. A test example is provided for which the displaced deterministic-stochastic estimates of the solution are advantageous.

1. Метод последовательных приближений

Рассмотрим функцию f (x), x = (x(1),..., x(1)), x E X С R1, из нормированного функционального пространства В с нормой ||f ||B и расстоянием dB(f1, f2) = ||f1 — f2||B, а также оператор

действующий из В в В: К : В ^ В. Выражение (1) определяет интегральный оператор Фредгольма с ядром к, Рассмотрим также интегральное уравнение Фредгольма второго рода для неизвестной функции

функция f € В и ядро к заданы. В качестве пространства В далее будем рассматривать В = С(X). Для решения уравнения (2) рассмотрим метод последовательных приближе-

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00046) и Президентской программы “Ведущие научные школы” (грант № НШ-4774.2006.1). © Институт вычислительных технологий Сибирского отделения Российской академии наук, 2006.

(1)

или ^ + f,

(2)

нии:

7™+1 = / + к1т, т = 0,1, 2,...,

с начальным элементом т0 = / Полагаем ||К||в < д < 1, В этом случае решение уравнения

(2) существует и единственно и справедливо неравенство [2]

дт

\\^~1ш\\в<~------||Т1-7о||в- (3)

1 - д

Решение уравнения (2) представимо в виде ряда Неймана

ГО

р(х) = / (х) + К/(х) + К 2/(х) + ... = ^ К*/ (ж), (4)

¿=0

где

К */ (х) = f ..^ / (Уо)к(Уо, У1)к(У1,У2) х ... X к(Уг- 1,х) ¿Уо ... йуг-1. (5)

X X

В силу равенства (5) и линейности операторов К* соотношение (4) представляет собой

х

2. Рандомизация приближения линейного функционала

Для многих теоретических и прикладных проблем (см,, например, [1]) важной является возможность приближенного вычисления линейных функционалов вида

I. = (^, Ь,)= <р(х)Н(х) ¿х. (6)

X

Комбинируя соотношения (4)—(6), имеем

ОО п п

'•-5.II

/(Уо)к(Уо, У1)к(У1, У2) X .. х к(У*-1, У*)Ь(У*) ¿Уо..¿У*-1 ¿У*. (7)

*=0 X X

Таким образом, функционал (6) представляет собой сумму интегралов бесконечно возрастающей кратности.

Можно задаться некоторым 5 < той взять приближение

*=о

Заметим, что

г(*’

^ - I.*1 < тах ^(х)| х ||^ - 7*||с(х), 7* = ^ К*/.

х *=о

В свою очередь, из неравенства (3) следует, что

д*

\\ч> - ъ\\с(х) < ~

*

Тогда величину 5 следует выбирать так, чтобы было выполнено неравенство

тах|ВД| х у— \\Kf- ¡\\с(х) < ё£1}. (9)

жеХ 1 — д

Для приближения интегралов из суммы 1(в, имеющих максимальную кратность I = ¡(в + 1) можно построить серию кубатурных формул 1 с числом узлов п1, обеспечивающим выполнение неравенства |1^ 11 < ё^. Положительные константы и

~(2) ’ 1

£в в сумме должны составлять величину допустимои погрешности £3.

Достаточно содержательной (и непростой) является задача оптимального согласованного выбора параметров в и щ. Трудности, в частности, связаны с оцеп кой величин д и \\К/ — /||с(х) в левой части неравенства (9), Кроме того, требуется учитывать суммарную трудоемкость Б (в, п1) получаемого алгоритма. Здесь применим подход, используемый при оптимизации дискретно-стохастических алгоритмов глобальной аппроксимации функций [3]: из уравнения ё^ + ё(2')(п1) = £ выражаем п1 через в и подставляем полученное выражение п1 = #(в) в формулу для Б. Далее исследуем функцию одного переменного Б (в, $(в)) на минимум, находя в = вШщ, При проведении практических расчетов полагаем в ~ 5Шщ, п1 ~ д(вшш). Пример соответствующего рассуждения для конкретной тестовой задачи из разд. 3 приведен в разд. 5,

Для малых £ и ё^ величипа I может быть достаточно большой, поэтому в качестве 1 следует выбирать простейшую стохастическую кубатурную формулу (т. е, метод Монте-Карло), При этом возникает проблема разумного выбора плотности соответствующего ¡-мерного случайного вектора. Согласно методу выборки по важности (см,, например, [1, 3]), наименьшая дисперсия получается для плотности

Р(У) = С

в г-1

^(7 ымы + £ /(Уо)к(уо, у1)к(у1, У2)-к(Уз> У'+ОМУ'+О

г=0 з=0

(10)

где у = (у0, у1,... ,ув) и С — соответствующая нормирующая константа. Выражение (10) является, как правило, достаточно громоздким, не дающим возможность построить эффективный алгоритм получения выборочных значений х случайного вектора х1 = (х]_,..., х{) согласно плотности р(у).

Отметим, что если требуется вычислить лишь одно слагаемое сумм (7) или (8) (Кг/, И), то использование плотности

Р(уо, у1,...,уг) » С/(уо)к(уо, у1)к(у1,у2)..к(уг-1, уг)И(уг)

может дать эффективный алгоритм метода Монте-Карло (см, разд. 3),

В качестве альтернативы подходу (8) рассмотрим стандартную несмещенную оценку по столкновениям [1, 3]:

4 = ЕС^е1 + --- + Сга, £ = У,(2Мхг). (и)

п 7=0

В соотношении (11) х0,х1,... ,хм — однородная цепь Маркова, обрывающаяся с вероятностью единица в состоянии хм (номер N — случайный) [3], начальной плотностью п(х) и переходной функцией

р(х', х) = т(х', х)(1 — ра(х')).

Здесь r(x', x) — плотность перехода, a pa(x') — вероятность обрыва траектории;

Л _ /Ы г, _ г, HXi-uXi)

Ц/о / \ 1 Wi Wi—i ( \ '

n(Xo) P(Xi-i ,Xi)

3. Сравнение подходов (8) и (11)

Численное сравнение подходов (8) и (11) проводилось на примере решения известного тестового интегрального уравнения второго рода из [1, 3], соответствующего анизотропному рассеянию при переносе излучения.

Рассмотрим следующую модель переноса малых частиц. Частицы двигаются из точки x = 0 в положительном направлении оси x случайными пробегами, длины которых распределены с плотностью e-x, x > 0, В конце пробега с вероятностью p частица поглощается, а с вероятностью q = 1 — p совершает очередной пробег xm-l ^ xm. Поставим задачу оценки вероятности P(H) вылета частицы за точку x = H.

Переходная функция для цепи столкновений определяется выражением p(x', x) = = q exp(—(x — x'))x(x — x'), где x(w) = 1 при w > 0 и x(w) = 0 при w < 0, Плотность начального состояния x0 равнa n(x) = e-x, x > 0. Для функции плотности столкновений y(x) можно записать уравнение

X

ф) = q j e-{x-X)^(x')dx' + e-x, x > 0, (12)

0

точное решение которого равно y(x) = e-px.

Рассмотрим локальную оценку (см,, например, [1, 3]), учитывая, что в данной ситуации происходит прямое моделирование (/(x) = n(x), k(x',x) = p(x',x)):

Ç(H) + + AH)

P(H) = <p(H) = Е£(я) + e~H « ^— + ••• + £* + е-я (13)

n

N

Ç(H) = q exP( — (H — xm))X(H — xm).

m=0

xN xN < H

точка, для которой выполнено неравенство xN > H. Функция h(x), определяющая вычисляемый функционал

Ih = (y, h) = P (H ) — e-H = EÇ(H ), (14)

равна h(x) = qexp(—(H — x))x(H — x).

Рассмотрим также метод последовательных приближений (8) для подсчета вероятности P (H). Погрешность такого приближения в метрике C можно оценить следующим образом:

£(l) = |Ц, _ I(s)||C = s = Il1 h Ih l|C =

4i=s+l

<

C

E Ki /

i=s+l

C.

C

E Kг/, h

Вычислим

Кг/(x) = i ... [ qie-X0e-(x-xi-1 ) x ... x e-(xi-Xo)x(x — xi-l) x ... x

^* X *

хх(х1 ~ хо) <1хо ... с1х^1 = ~—е~х.

Далее имеем

£ к*ч

(=5 + 1

С

тах

х£[0,Н ]

£

*=5 + 1

~гГ

< тах

х€[0,Н]

е хд5+1ея

X

5 + 1

(в + 1)!

5 Е [0, х].

одесь использовано то обстоятельство, что выражение > —— является остатком разло-

А—/ Ц

жения Тейлора формулы е9х. Следовательно,

*=5 + 1

£ кч

*=5 + 1

<

С

д5+1н 5+1 (5 + 1)'

д5 + 2 Н5 + 1

с = тах \де~{н~х)\ = д, < —---------------------—(15)

х£[0,И ]

(в + 1)! '

Выражения (К*Ч, К) представляют собой интегралы увеличивающейся (с ростом г) размерности, Для их приближенного вычисления применим алгоритм выборки по важности из [4], Необходимое количество интегралов в в сумме (8) можно найти исходя из требуемой погрешности, преобразовав неравенство для е5^ в равенство и воспользовавшись формулой Стирлинга:

‘ £(1) = д3+2{еН)3+1

5 у/27г(^+Т)('5 + 1)5+1

Отметим, что такая оценка числа слагаемых завышенная. Минимальное количество слагаемых можно найти исходя из известного точного решения уравнения (12):

=■(1)

е?' <

е—рН- е~н

д*Н *

*=0

(и)

—х

е

о

эо

5

В проведенных численных экспериментах для приближения интегралов выбирались плотности, являющиеся кусочно-постоянными приближениями подынтегральных функций, с числом т промежутков постоянства вдоль одной координаты не больше десяти, В частности, для интегралов малой размерности (г < 6) выбиралось т = 10, для интегралов большей кратности т < 10, В табл. 1 приведены затраты па вычисление вероятности Р(Н) вылета частицы за границу Н (в секундах) при использовании локальной оценки (строки I.осКн!) и метода последовательных приближений (строки ЯН'рнрр) для разных

Таблица 1. Сравнение трудоемкости метода последовательных приближений и локальной оценки

Я д = 0.1 д = 0.3 д = 0.5

1, LocEst 1.77 6.0 24.0

1, Stepapp 0.215 1.51 19.5

10, LocEst 5.9 5.51 1.52

10, Stepapp 0.11 2.3 29

20, LocEst 618 4886 102

20, Stepapp 0.1 _ _

Н д Р(Н)

е1

е1 < 0.01 ■ Р(Н) = 0.01 (е—рН + е—н) .

Для локальной оценки достижение соответствующей величины контролировалось выражением €\ < у/Щя/\Д^ гДе п2 — число моделируемых траекторий; величина вычисляется точно [1], В методе последовательных приближений количество интегралов з

г

е = е1/(2з), г = 1,..., з.

Проведенные расчеты позволяют сделать следующий вывод, В модели переноса частиц с анизотропным рассеянием алгоритм, связанный с использованием локальной оценки для подсчета вероятности Р (Н), уступает по трудоемкости методу последовательных приближений лишь в случае малой вероятности Р (Н) выхода за границу Н и малых значений вероятности рассеяния д < 0.5, т, е, при вычислении небольшого количества интегралов (з < 10). Для д > 0.5 трудоемкость локальной оценки меньше. Отметим, что если требуется вычислить лишь одно слагаемое суммы (8) (К*4, И), то использование метода выборки по важности может быть очень эффективно.

4. Использование специального функционала

Следует заметить, что функционал (14) дает не самую убедительную иллюстрацию преимущества метода последовательных приближений (8) из-за отсутствия “длинных” тра-

Н

е

погрешность; п2 и — число моделируемых траекторий и время реализации для приближения (13) и соответственно щ, ^ — для приближения (8); з — номер обрыва траектории.

Более показательные результаты дает выбор функции К(х), распределенной вдоль всей положительной полуоси, например,

И(х) = е—Ах, х > 0 (17)

(однако при этом теряется “физический смысл” функционала 4), Для приближения со-

5

ответствующего функционала, (8) использовалась оценка £(5) = И(х*), причем отрезок

*=0

цепи Маркова х0, х1,... ,х5 получался с помощью прямого моделирования.

Из табл. 3 видно, что для функционала (6) с функцией (17) метод последовательных приближений более эффективен, чем стандартный алгоритм с оценкой по столкновениям. По аналогии с рассуждениями при получении соотношения (15) можно рассмотреть

Таблица 2. Сравнение приближений (8) и (13) для функционала (14)

Я д е «2 ¿2 «1 ¿1 в

30 0.999 0.0032 55 000 0 : 09 60 000 0 : 10 45

70 0.999 0.0074 20 000 0 : 13 21 000 0 : 14 91

30 0.99 0.006 20 000 0 : 05 25 000 0 : 06 43

70 0.99 0.005 25 000 0 : 12 27 000 0 : 13 87

погрешность \'|.гл функцію нала Ih с функцией (17), Вычислим

(Kif, h) = j ... J qie-x0 e-(xi-xi-l) x ... x e-(xi-x0 )x(xi - xi-1) x 0 0

+o / x 1

xx(xi ~ xo)e~ÂXi dxо ... dxi-i dxi = q% J e~Xi^A+1'> ( J dxi-i ] dxi =

00

СО

= £ fe-x<W xi dx, qi

г! ] * * (А +1)*+^

0

Здесь использовано известное свойство гамма-функции

СЮ

т(") = 1 и}и—1е—ш dw : Г(г + 1) = г!.

0

Используя формулу суммы членов бесконечно убывающей геометрической прогрессии, имеем

£(1) = |Ц, _ I(s) ||C = s = llJh J-h l|C =

О

£

i=s+1

K if,h

C

A + 1) A +1 - q'

Общую погрешность £s алгоритма (8) перепишем в виде суммы

= єі1} + єі2) ~ ( —і ~ г- I -г—:-----------Ь

s s s

q Y+1 1 d(s)

где d(s) = v' DC' При проведении численных экспериментов рассмотрены следующие предельные случаи,

1, При s ^ то подтверждено, что £s ведет себя как погрешность для метода Монте-Карло без обрыва.

Таблица 3. Приближения (8) и (13) для функционала (6) с функцией h(x) = e Ax

А Q є «2 ¿2 Пі il s

0.5 0.999 0.01 50 000 0 : 43 300 000 0 : 03 19

0.1 0.99 0.023 60 000 0 : 05 65 000 0 : 03 60

0.3 0.999 0.014 50 000 0 : 43 1 000 000 0 : 19 21

Таблица 4. Близость погрешностей es и es^ при щ >> 1

s £d)

3 0.392955 0.392698

6 0.116376 0.116006

10 0.023428 0.022823

12 0.010377 0.010123

16 0.002193 0.001991

2, При щ ^ то подтверждено, что погрешность є зависит только от з (т. е, £ = £

_ (!)' з ъ, £ — £.3 ,

Соответствующие результаты расчетов для д — 0.999 А — 0.5, п1 — 107 представлены в виде табл. 4.

5. Условная оптимизация алгоритма

Продемонстрируем на рассматриваемом тестовом примере работу сформулированной в разд. 2 процедуры согласованного выбора параметров з и п1 в алгоритме (8). Заметим, что трудоемкость этого алгоритма пропорциональна величине £ = вп1. Задаем уровень погрешности £ и рассматриваем равенство

I = »(1) + ¡¡т = ('+‘___________1____+ М (18)

* 1 + ■ и+1/ А+1-« \/йГ УЩ

Полагаем ¿(в) ~ й — еопв^ так как при малом значепии £ выполнено

3 + 1'

2А + 1

Выражая п1 через в из соотношения (18), имеем

=<е 3/(ё-- (атт) ' (19)

Последняя функция численно исследовалась на минимум по в. Поиск минимума сводится к решению уравнения

Это условие равенства нулю производной функции (19).

В частности, для £ — 0.01, д — 0.999, А — 0.^ и й — 2.1 удалось установить, что минимум функции £(,§) достигается примерно при в — 19. Заметим также, что для в — 19 получилось совпадение теоретического и практического результатов (т. е, при в — 19 для получения заданного уровня погрешности потребовалось минимальное время счета — табл. 5).

Таблица 5. Оптимизация алгоритма (8) по в

в «і і, с

13 469 223 5.55

16 77 651 1.10

18 59 919 0.94

19 56 221 0.93

20 53 944 0.94

22 51 577 0.99

30 49 867 1.26

40 49 801 1.70

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

Список литературы

[1] Ермаков С.М., Михайлов Г.А. Статистическое моделирование. М.: Наука, 1982.

[2] Канторович Л.В., Акилов Г.П. Функциональный анализ. М.: Наука, 1977.

[3] Михайлов Г.А., Войтишек A.B. Методы Монте-Карло. М.: Академия, 2006.

[4] Войтишек A.B., Каблукова Е.Г. Исследование адаптивных дискретно-стохастических алгоритмов численного интегрирования // Матер. VI Междунар. семинара-совещания “Кубатур-ные формулы и их приложения”. Уфа, июль 2001. С. 46-52.

Поступила в редакцию 15 сентября 2006 г.