Метод монте карло для уравнений

Решение интегральных уравнений методом Монте-Карло

Южно-Российский государственный политехнический университет

NovaInfo53, с. 9-12
Опубликовано 26 октября 2016
Раздел: Физико-математические науки
Просмотров за месяц: 3
CC BY-NC

Аннотация

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

Ключевые слова

Текст научной работы

Метод Монте-Карло находит широкое применение в практике решения вычислительных задач, в том числе при решении интегральных уравнений 1. В то же время не все возможности данного метода используются полностью. Так при решении интегральных уравнений преимущественно используется вариант этого метода, основанный на суммировании резольвенты и использовании цепей Маркова [1]. Это обстоятельство значительно сужает класс решаемых задач, так как требуется ограничение нормы интегрального оператора единицей. В данной статье с целью восполнения отмеченного пробела рассматривается применение метода Монте-Карло в его классической форме к решению широкого круга интегральных уравнений типа Фредгольма.

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

Описание метода

Рассмотрим интегральное уравнение Фредгольма:

\mu u(x) — \lambda \int_V K(x,y) u(y) dy = f(x), x \in V

где u(x) — искомая функция, x=(x_1, \dots, x_m), y=(y_1, \dots, y_m)

— точки области V из m-мерного евклидова пространства, \mu

— некоторые вещественные или комплексные числа, K(x,y) — ядро интегрального оператора, f(x) — свободный член.

Предположим, что известны n точек области V: y^1=(y_1^1, \dots,y_m^1), \dots, y^n=(y_1^n, \dots, y_m^n)

, полученные из распределения с плотностью p(y), y \in V

Интеграл в (1) можно приближенно вычислять при помощи традиционной схемы вычисления интегралов методом Монте-Карло [1]:

\int_V K(x,y) u(y) dy \approx 1/n \sum_^n S_j (x), x \in V

где S_j (x) = K(x,y^j) u(y^j) /p(y^j)

Перепишем (1) в эквивалентном виде:

\mu u(x) — \lambda/n \sum_^n S_i (x) — \lambda R_n (x) = f(x), x \in V

— остаточный член формулы интегрирования Монте-Карло:

\int_V K(x,y) u(y) dy = 1/n \sum_^n S_j (x) + R_n (x)

Используем точки y^1=(y_1^1,\dots,y_m^1),\dots, y^n=(y_1^n,\dots,y_m^n)

как узлы коллокаций в известном вычислительном методе, при помощи которого получим из (2) соответствующую СЛАУ для нахождения приближенных значений решения в рассматриваемых точках:

\mu u_i — \lambda/n \sum_^n K(x^i,y^j)/p(y^j) u_j = f_i, i=1,\dots,n

Поскольку остаточный член квадратурной суммы метода Монте-Карло с любой наперед заданной вероятностью стремится к нулю при стремлении числа узлов к бесконечности, то обоснованно предполагать, что при достаточно гладком ядре и ограниченности оператора, обратного к оператору интегрального уравнения (1), решение СЛАУ (3) сходится к точному в одной из вероятностных мер. В литературе соответствующие вопросы сходимости детально рассмотрены применительно к задаче суммирования ряда Неймана [1].

Пример применения метода

Исходные данные модельной задачи:

K(x,y)= x_1 \dots x_m, y_1 \dots y_m

f(x)=x_1 \dots x_m + g (x_1 \dots x_m) ^2; g=10

Область интегрирования — m-мерный куб

u(x)=c_0 x_1 \dots x_m + g x_1 \dots x_m (x_1 \dots x_m + c_1)

c_0=0.5, c_1=c_0 (3/4)^m, \lambda =-3^m, \mu =1

Результаты трех последовательных вычислений решения при размерности области m=10 и числе узлов n=10:

• точное решение (0.02375, 0.01527, 0.00363, 0.04009, 0.04210, 0.08175, 0.00694, 0.03348, 0.03155, 0.01348);
• приближенное (0.02679, 0.01758, 0.00448, 0.04425, 0.04638, 0.08800, 0.00831, 0.03722, 0.03516, 0.01561).

Погрешность решения в норме l₁ около 11%.

На двух последующих вычислениях решения погрешность много меньше: около 1% и 0.4%.

Читайте также

Средства стохастической подготовки обучающихся на основе информационных технологий

Инструментальная реализация прикладной математической подготовки бакалавра экономики и менеджмента

1. Синчуков А.В.

NovaInfo59, с.24-28, 13 февраля 2017 , Физико-математические науки, CC BY-NC

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

1. Ермаков С.М. Метод Монте-Карло в вычислительной математике (вводный курс). Санкт-Петербург: Издательство Бином. 2011. 192 с.
2. Некрасов С.А., Ткачев А.Н. Теория вероятностей и ее приложения: Учеб. пособие/ Юж.-Рос. гос. техн. ун-т. – Новочеркасск: ЮРГТУ, 2007. 148 с.
3. Некрасов С.А. Методы ускоренного статистического моделирования и их применение в электротехнических задачах/ Изв. вузов. Электромеханика. — 2008. — No 5. — С. 13 — 19. http://elibrary.ru/item.asp?id=12159957

Цитировать

Некрасов, С.А. Решение интегральных уравнений методом Монте-Карло / С.А. Некрасов. — Текст : электронный // NovaInfo, 2016. — № 53. — С. 9-12. — URL: https://novainfo.ru/article/8242 (дата обращения: 21.02.2022).

Поделиться

Электронное периодическое издание зарегистрировано в Федеральной службе по надзору в сфере связи, информационных технологий и массовых коммуникаций (Роскомнадзор), свидетельство о регистрации СМИ — ЭЛ № ФС77-41429 от 23.07.2010 г.

Соучредители СМИ: Долганов А.А., Майоров Е.В.

Алгоритм Монте-Карло для решения систем линейных алгебраических уравнений методом Зейделя∗ Текст научной статьи по специальности « Математика»

Аннотация научной статьи по математике, автор научной работы — Товстик Татьяна Михайловна, Волосенко Ксения Сергеевна

Для решения системы линейных алгебраических уравнений методом Монте-Карло используется алгоритм последовательных приближений. Очередная итерация моделируется в виде случайного вектора, математическое ожидание которого совпадает с приближением процесса итерации в форме Зейделя. Выводится система линейных уравнений, которым удовлетворяют взаимные корреляции компонент предельного вектора и корреляции двух последовательных приближений. Доказывается существование и конечность предельных дисперсий случайного вектора решений системы. Библиогр. 7 назв. Табл. 1.

Похожие темы научных работ по математике , автор научной работы — Товстик Татьяна Михайловна, Волосенко Ксения Сергеевна

MONTE CARLO ALGORITHM FOR SOLUTION OF A SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS BY THE SEIDEL’S METHOD

For solution of a system of linear algebraic equations by the Monte Carlo method a method of the successive approximations is used. The next in turn iteration is modeled as a random vector, the expectation of which coincides with the corresponding Seidel’s iteration. The existence and the boundedness of dispersions of the limiting vector is proved. The system of linear equations, to which satisfy the cross-correlations of the limiting vector, and the limiting correlations of two successive approximations is delivered. All these correlations are unknowns of the common linear system and may be found only simultaneously. The boundedness of dispersions of the limiting vector leads to the stochastic stability of algorithm. The number of iterations giving the desirable exactness is estimated. As an examples, the systems of the orders of n = 3 and n = 100 are studied. Refs 7. Tables 1.

Текст научной работы на тему «Алгоритм Монте-Карло для решения систем линейных алгебраических уравнений методом Зейделя∗»

Вестник СПбГУ. Сер. 1. Т. 3(61). 2016. Вып. 3

ДЛЯ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ

АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ МЕТОДОМ ЗЕЙДЕЛЯ*

Т. М. Товстик, К. С. Волосенко

Санкт-Петербургский государственный университет,

Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9

Для решения системы линейных алгебраических уравнений методом Монте-Карло используется алгоритм последовательных приближений. Очередная итерация моделируется в виде случайного вектора, математическое ожидание которого совпадает с приближением процесса итерации в форме Зейделя. Выводится система линейных уравнений, которым удовлетворяют взаимные корреляции компонент предельного вектора и корреляции двух последовательных приближений. Доказывается существование и конечность предельных дисперсий случайного вектора решений системы. Библиогр. 7 назв. Табл. 1.

Ключевые слова: алгоритм Монте-Карло, метод Зейделя, система линейных алгебраических уравнений.

1. Введение. Метод Монте-Карло для решения системы линейных алгебраических уравнений рекомендуется использовать в том случае, если ее порядок достаточно велик [1]. В монографиях [1, 2] с помощью метода Монте-Карло вычисляется одна компонента вектора решения или скалярное произведение вектора решения и произвольного вектора. В работах 3 оценивается сразу весь вектор решений. В работе [3] используется стохастический метод итераций, причем математическое ожидание случайного вектора очередной итерации совпадает с суммой соответствующего отрезка ряда Неймана. В статьях [4, 5] представлены алгоритмы Монте-Карло, позволяющие получать итерационное решение системы линейных уравнений при условиях более слабых, чем условия обычного метода Монте-Карло. В данной работе, в отличие от алгоритма в [3], особенность построения случайного вектора заключается в том, что математическое ожидание очередной итерации совпадает с итерацией Зейделя [6]. В дополнение к работе [7] доказывается существование и конечность предельных дисперсий случайного вектора решений системы.

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

где А = — квадратная матрица, а X = (Хь . Хп)т и / = (/ь . fn)T —

В работе рассматривается алгоритм Монте-Карло, в основе которого лежит метод Зейделя. Будем предполагать, что для нормы матрицы A выполняется условие

IIAll = max ]Т\Aik| 1.

2. Алгоритм Монте-Карло для метода Зейделя. Рассмотрим метод, который имитирует алгоритм Зейделя. Будем моделировать случайные векторы то > 0 такие, что выполняется £(0) = /, а при то > 1 справедливо Е £(т) = X(т), причем компоненты X(т) удовлетворяют равенствам (3). В основе моделирования лежит имитация цепи Маркова, которая задается вектором начального распределения п и стохастической матрицей Р вероятностей перехода:

\\Pij||, Р. > ° 1 0 для тех пар (г,]), для которых справедливо Аг. = 0. На матрицу Р не накладывается никаких ограничений. В частности, она может состоять из отдельных блоков.

Моделирование вектора С(т) осуществляется следующим способом. Выбирается начальное приближение С(0) = /, далее при то > 1 компоненты случайного вектора обновляются последовательно, начиная с г = 1, причем при моделировании г-й компоненты учитывается результат обновления всех предыдущих компонент ((т\ ] 1 моделируются последовательно, начиная от г = 1 до п по формулам

(т) , 1 ^АЪ /РЪ пРи ^ г,

где случайные состояния ^ разыгрываются в соответствии с распределением (5). Формулы (6) допускают представление

Лг) Ai,ji Am) i Ai,j Am— 1) , x 1 г»

Q =2^X(k,ji)-j2-Cji ‘+2^X(k,ji)-j2-Cji > + fi, г = 1,2. п, (7)

где X(k,p) — символ Кронекера.

Вычисление математического ожидания от обеих частей равенства (7) приводит к равенствам

ЕСГ = £ Aij EZ(m) + ^ ECjm—1) + fi, г = 1, 2. n, m > 1. (8) j=1 j=i

Отсюда при m =1 получаем

EZ(1) = £ A1j fj + f1, EZ(1) = £ Aij Ej + E Aij fj + fi, г = 2. n. (9)

Из формул (9) следуют равенства EZ(1) = X(1), EZ(1) = X(1), в которых X(1) те же, что и в (3) метода Зейделя, если начальный вектор удовлетворяет X(0) = f. Нетрудно проверить, что при всех m > 1 выполнены равенства EZ(m) = X\m\ Теорема 1. Пусть для системы линейных алгебраических уравнений вида

k=1 Piji k=i Piji

Тогда справедливы равенства

ECr = ХГ = £ AijX(m) + ^ Aijx(m—1) + fi, 1 0 при то ^ ж для всех г. Так как имеет место равенство Е(г(т) = Х(т), выполняется \Хг — Е^™^ —> 0 при то ^ ж, тем самым теорема доказана.

Для оценки решения уравнения (1) следует достаточно большое число N

раз моделировать случайный вектор С(М), тогда будем иметь ес(м) « Не можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Теорема 2. Пусть выполнены условия теоремы 1 и норма матрицы В = [В. ] = [А. /ргз] удовлетворяет неравенству

\\В\\ = та? £Вш 1. Пользуясь равенствами (6), находим второй момент разности

/ Ат) с \ тл I ^(т)\ Лг(т)

(С /г) и с учетом Е I Сг ) = Х> получаем

2 г — 1 Л2 2 п л2 2

Е (с(т)) =Е—Е(Сзт)) +Е—Е(^1)) то > 1. (18)

Равенство (18) дает

КГ] — КТ—1) = £ Вгз (Ц) — К.— 1)) + £ Вгз — 1) — Я.^ +

+ 2fi (X(m) — X(m-1)) . (19)

Покажем, что для всех г будет выполняться \Я(т) — Я™ ^^ 0 при то ^ ж.

В [6] доказано, что в методе Зейделя для разностей итераций X(т), X(т—1) справедливы неравенства

||Х(т) _х(т-1)|| М 0 такие, что выполняется ат Не можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

Пользуясь равенствами (23) и симметрией ковариаций Д™1^, получаем неравенства, аналогичные (20),

Пример 1. Рассмотрим систему линейных алгебраических уравнений вида (1), где матрица А и вектор f определяются выражениями

0.3 —0.5 0.1 A = I —0.2 0.3 0.4 0.4 -0.3 0.2

Норма матрицы А равна ||А|| = ш&х^ | = 0.9. Элементы матрицы Р вычис-

лялись по формулам

1 М*, а при больших N N = 106 и больше), наоборот, брать М Надоели баннеры? Вы всегда можете отключить рекламу.

VMath

Инструменты сайта

Основное

Навигация

Информация

Действия

Вспомогательная страница к разделу СИСТЕМЫ ЛИНЕЙНЫХ УРАВНЕНИЙ
Для понимания материалов этого пункта полезно ознакомиться с идеологией метода Монте-Карло

Решение системы линейных уравнений методом Монте-Карло

Рассмотрим систему из $ n_<> $ линейных уравнений относительно $ n_<> $ неизвестных $$ \left\< \begin a_<11>x_1 &+a_<12>x_2&+ \ldots&+a_<1n>x_n &=b_1,\\ a_<21>x_1 &+a_<22>x_2&+ \ldots&+a_<2n>x_n &=b_2,\\ \dots & & & & \dots \\ a_x_1 &+a_x_2&+ \ldots&+a_x_n &=b_n, \end \right. $$ которую иногда будем представлять и в матричном виде $$ AX= \mathcal B \ . $$ Решение этой системы равносильно нахождению минимума квадратичной функции $$ F(X)=\sum_^n \alpha_j \left(a_x_1+\dots+a_x_n-b_j \right)^2= (AX-\mathcal B)^ <\top>\left( \begin \alpha_1 & 0 & \dots & 0 \\ 0 & \alpha_2 & \dots & 0 \\ \vdots & & \ddots & \vdots \\ 0 & 0 & \dots & \alpha_n \end \right) (AX-\mathcal B) \ , $$ где $ \< \alpha_j\>_^n $ — положительные числа, а $ <>^ <\top>$ означает транспонирование.

Если исходная система линейных уравнений имеет единственное решение $ X=X_<\ast>=(x_<1\ast>,\dots, x_) $, то в пространстве $ \mathbb R^ $ уравнение $$ F(X) = 1 $$ задает эллипсоид с центром в точке $ X_ <\ast>$. Каждая из $ (n_<>-1) $-мерных гиперплоскостей $ x_k=x_ $ ( линейных многообразий), проходящих через центр эллипсоида, делит его объем пополам.

Построим $ n_<> $-мерный параллелепипед $$ A_1 ☞ ЗДЕСЬ: $$ V_ <\mathrm E>\approx 3133.207748 \ . $$ $$ \begin N & 1000 & 5000 & 10000 & 20000 & 50000 \\ \hline M & 62 & 297 & 581 & 1181 & 2885 \\ \hline V_ <\Pi>M/N & 3339.6300 & 3199.581 & 3129.5565 & 3180.7282 & 3108.0105 \\ \hline \overline <\xi_1>& 3.2592 & 3.2745 & 3.1230 & 3.1020 & 3.1798 \\ \hline \overline <\xi_2>& 0.9406 & 1.5346 & 1.4669 & 1.4867 & 1.5141 \\ \hline \overline <\xi_3>& 0.3453 & 0.5163 & 0.3996 & 0.5001 & 0.4299 \\ \hline \overline <\xi_4>& -1.7029 & -2.2198 & -2.1676 & -2.1545 & -2.2413\\ \end $$ Решение системы $$ x_1=\frac<257> <84>\approx 3.05952,\ x_2=\frac<53> <36>\approx 1.47222,\ x_3=\frac<55> <126>\approx 0.43651 , x_4=-\frac<547> <252>\approx -2.17063 \ . $$

Источник

Бусленко Н.П., Шрейдер Ю.А. Метод статистических испытаний Монте-Карло и его реализация в цифровых машинах. М.: Физматгиз, 1961. 266 с.