ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.
TECHNICAL SCIENCE. 2021. No 1
УДК 517 DOI: 10.17213/1560-3644-2021-1-27-37
МЕТОДЫ РЕШЕНИЯ ИНТЕГРАЛЬНЫХ УРАВНЕНИЙ НА СЛУЧАЙНОЙ И ПСЕВДОСЛУЧАЙНОЙ СЕТКЕ И ИХ ПРИМЕНЕНИЕ В ПРИКЛАДНЫХ ЗАДАЧАХ
© 2021 г. Абас Висам Махди Абас1, Р.В. Арутюнян2
1Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия, 2Московский государственный технический университет имени Н.Э. Баумана, г. Москва, Россия
THE CALCULATION OF SOLUTIONS OF INTEGRAL EQUATIONS ON RANDOM AND PSEUDORANDOM GRID AND ITS APPLICATION
Abas Wisam Mahdi Abas1, R. V. Harutyunyan2
1Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia, 2Bauman Moscow State Technical University, Moscow, Russia
Абас Висам Махди Абас - аспирант, кафедра «Прикладная математика», Южно-Российский государственный политехнический университет (НПИ) имени М.И. Платова, г. Новочеркасск, Россия. E-mail: [email protected]
Арутюнян Роберт Владимирович - канд. физ.-мат. наук, доцент, кафедра «Вычислительная математика и математическая физика», Московский государственный технический университет имени Н.Э. Баумана, г. Москва, Россия. E-mail: [email protected]
Рассматривается подход на основе метода случайных кубатур для решения как одно-, так и многомерных сингулярных интегральных уравнений, уравнений Вольтерра и Фредгольма 1 рода, для некорректных задач теории интегральных уравнений и т.д. Изучен вариант квази-Монте-Карло для рассматриваемого метода. Интеграл в интегральном уравнении приближенно вычисляется при помощи традиционной схемы вычисления интегралов методом Монте-Карло. Применяется многомерная интерполяция на произвольном множестве точек. Рассмотрены примеры применения метода к одномерному интегральному уравнению с гладким ядром с использованием как случайных, так и низкодисперсных псевдослучайных узлов. Решено с помощью метода Ньютона многомерное линейное интегральное уравнение с полиномиальным ядром, многомерная нелинейная задача - интегральное уравнение Гаммерштейна. Показано существование нескольких решений. Рассмотрены многомерные интегральные уравнения первого рода и их решение с использованием регуляризации. Решение методами Монте-Карло и квази-Монте-Карло подобных задач в изученной литературе не проводилось. Был использован метод регуляризации Лаврентьева, а также случайные и псевдослучайные узлы, полученные при помощи последовательности Хальтона. Решена проблема собственных значений. Установлено, лучшим из рассмотренных методов является метод Леверье-Фаддеева. Результаты решения задачи для различного числа квадратурных узлов представлены в таблице. Исследован подход на основе параметрической регуляризации ядра, интерполяционно-проекционный метод, усредненные адаптивные плотности. Решены пространственные краевые задачи Дирихле для уравнения Лапласа для шаровой и тороидальной областей.
Рассматриваемые подходы позволяют расширить круг задач теории интегральных уравнений, решаемых методами Монте-Карло и квази-Монте-Карло, поскольку отсутствуют ограничения на величину нормы интегрального оператора. Рассмотрена серия примеров, демонстрирующих степень эффективности исследуемого метода.
Ключевые слова: интегральные уравнения; высокая размерность; метод Монте-Карло.
Abas Wisam Mahdi Abas - Graduate Student, Department «Applied mathematics», Platov South-Russian State Polytechnic University (NPI), Novocherkassk, Russia. E-mail: [email protected]
Harutyunyan Robert V. - Candidate of Physical and Mathematical Sciences, Associate Professor, Department «Computational mathematics and mathematical physics», Bauman Moscow State Technical University, Moscow, Russia. E-mail: [email protected]
ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION. TECHNICAL SCIENCE. 2021. No 1
The article considers an approach based on the random cubature method for solving both single and multidimensional singular integral equations, Volterra and Fredholm equations of the 1st kind, for ill-posed problems in the theory of integral equations, etc. A variant of the quasi-Monte-Carlo method is studied. The integral in an integral equation is approximated using the traditional Monte-Carlo method for calculating integrals. Multidimensional interpolation is applied on an arbitrary set of points. Examples of applying the method to a one-dimensional integral equation with a smooth kernel using both random and low-dispersed pseudo-random nodes are considered. A multidimensional linear integral equation with a polynomial kernel and a multidimensional nonlinear problem - the Hammerstein integral equation-are solved using the Newton method. The existence of several solutions is shown. Multidimensional integral equations of the first kind and their solution using regular-ization are considered. The Monte Carlo and quasi-Monte-Carlo methods have not been used to solve such problems in the studied literature. The Lavrentiev regularization method was used, as well as random andpseu-do-random nodes obtained using the Halton sequence. The problem of eigenvalues is solved. It is established that one of the best methods considered is the Leverrier-Faddeev method. The results of solving the problem for a different number of quadrature nodes are presented in the table. An approach based on parametric regularization of the core, an interpolation-projection method, and averaged adaptive densities are studied. Spatial Dirichlet boundary value problems for the Laplace equation for spherical and toroidal regions are solved.
These approaches allow us to expand the range of problems in the theory of integral equations solved by Monte Carlo and quasi-Monte-Carlo methods, since there are no restrictions on the value of the norm of the integral operator. A series of examples demonstrating the effectiveness of the method under study is considered.
Keywords: integral equation; high dimension; Monte-Carlo method.
Введение
Метод Монте-Карло находит широкое применение в практике решения вычислительных задач, в том числе при решении интегральных уравнений (ИУ) [1 - 3]. В то же время не все возможности данного метода используются полностью. Так при решении интегральных уравнений преимущественно используется вариант этого метода, основанный на суммировании резольвенты и использовании цепей Маркова [1]. Это обстоятельство значительно сужает класс решаемых задач, так как требуется ограничение нормы интегрального оператора единицей.
В работе [4] предложен подход, названный полустатистическим и основанный на применении квадратурной формулы со случайными узлами. Модернизация этого метода на основе последовательного усреднения результатов нескольких серий вычислений описана в [5]. В [2] рассматривается метод решения регулярных ИУ Фредгольма II рода, сочетающий метод Монте-Карло с проекционными методами. Решение ИУ находится в виде разложения по ортонормиро-ванному базису с учетом погрешности разложения. Интегралы, содержащие неизвестную погрешность, аппроксимируются формулами метода Монте-Карло, в результате чего образуется СЛАУ относительно коэффициентов разложения и значений погрешности в узлах случайной сетки интегрирования.
Рассматриваемый подход является еще малоизученным, так как практически не применялся для сингулярных ИУ, уравнений Вольтерра и ИУ
Фредгольма I рода, для некорректных задач теории ИУ и т.д. Круг решенных задач также является относительно узким (так в основном решены одно- и двумерные задачи). Практически не изучен вариант квази-Монте-Карло для данного метода. Эти и другие вопросы рассматриваются в данной статье.
Описание метода
Рассмотрим ИУ Фредгольма:
цы ( х)-А К ( х, у)ы (у) йу = / ( х), х еК, (1)
где и(х) - искомая функция, точки области V из т - мерного евклидова пространства; ц и X -некоторые вещественные или комплексные числа; К(х, у) - ядро интегрального оператора; Дх) - свободный член.
Предположим, что известны п случайных точек области V:у1 = (у1,...,ут)уп = (у",...,упт), полученные из распределения с плотностью р(у), yеV. Условие нормировки:
| Р (У ) йу =1.
V
Интеграл в (1) можно приближенно вычислить при помощи традиционной схемы вычисления интегралов методом Монте-Карло [1]:
1 "
IК(х,у)ы(у)йу (х), хе V,
V п}=1
где Sj (x) =
K (x, У )u( У)
p(yJ) '
ISSN 15б0-3б44 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.
TECHNICAL SCIENCE. 2021. No 1
Перепишем (1) в эквивалентном виде:
- п
ци (х) - - Х$ (х) - ХИп (х) — / (х), х 6 V, (2) П 7—1
где Яп(х) - остаточный член формулы интегрирования Монте-Карло:
1 п
IК(х,у)и (у)йу — - ££/ (х) + Ип (х).
V п ] —1
Используем точки у1 — (у-,...,у1т) ,...,
п / п п \
у — (у1 ,..., ут) как узлы коллокаций в известном вычислительном методе, при помощи которого получим из (2) соответствующую СЛАУ для нахождения приближенных значений решения в рассматриваемых точках:
Ни (х)--±К±}1 и) — ,,
п7—1 Р(у]) } (3)
и ~ и(у), 7 — 1,...,п.
Поскольку остаточный член квадратурной суммы метода Монте-Карло с любой наперед заданной вероятностью стремится к нулю при стремлении числа узлов к бесконечности, то обоснованно предполагать, что при достаточно гладком ядре и ограниченности оператора, обратного к оператору интегрального уравнения (1), решение СЛАУ (3) сходится к точному в одной из вероятностных мер. В литературе соответствующие вопросы сходимости детально рассмотрены применительно к задаче суммирования ряда Неймана [1] и для так называемого полустатистического метода [2, 4, 5].
Многомерная интерполяция
Интерполирование на случайной сетке используется в методе случайных квадратур [1], при вычислении адаптивной плотности классическим методом или методом отбора, при сглаживании, для графического отображения результатов и т.д.
В многомерном случае соотношения интерполяционного метода для аппроксимации функции п переменных и(х) на произвольном наборе точек имеют вид:
4 (х;ul,...,ип)="=^]У(х -у],И),
Ьп (хг;и1,...,ип) — и,7 — 1,.,п
"=У]У (Х7 - у] ,И ) — и7,7 — 1,.,п
, n, , n,
или финитные «пирамидальные» кусочно-линейные функции на основе функции
Ф (x) = max (0,1 - 2 Jx| / h), Wj - неизвестные коэффициенты.
Примеры применения метода
Одномерное ИУ с гладким ядром. Исходные данные модельной задачи:
K(x,y)=sin(x)- sin(y); f (x) = 1. Промежуток
интегрирования u ( x ) = ( nsin ( x ) -1) /
[0, я].
П
-3
Точное решение ; X = 1, ц = 1. Реше-
ние осуществлялось методами Монте-Карло для обычной равномерно распределенной псевдослучайной сетки и для соответствующей модифицированной низкодисперсной последовательности, квази-Монте-Карло (узлы сетки - точки низкодисперсной последовательности Хальтона НО, 1 = 1, ..., п; 5 - простое число), а также квадратурным методом центральных прямоугольников. Результаты представлены на рис. 1, 2, где сплошная линия соответствует точному решению.
1,2
0,96 0,72 0,48 0,24
0
-0,24 -0,48
1,2
0,96 0,72 0,48 0,24
1,3 1,9 2,5 3,1- 0 -0,24
-0,48
Рис. 1. Решение ИУ (1) методом Монте-Карло. Слева: n = 10, псевдослучайные узлы; справа: n = 10, модифицированная псевдослучайная последовательность
/ Fig. 1. Solution of IE(1) by the Monte-Carlo method. Left: n = 10, pseudo-random nodes; right: n = 10, modified pseudo-random sequence
1,2
0,96
0,72
0,48
0,24
0
-0,24 -0,48
1,2
0,96
0,72
0,48
0,24
3,1 у 0 -0,24
-0,48
0,6 1,3 1,9 2,5 3,1 • \
где v(x,h) - базисная функция, например, функции Гаусса v(x,h) = exp(-x2 I (2h2 )), h да n-m, m > 0,
Рис. 2. Численное решение ИУ (1). Слева: n = 10, последовательность Хальтона Hi (2); справа n = 10, формула центральных прямоугольников / Fig. 2. Numerical solution of the IE (1). Left: n = 10, the sequence of Halton H(2); on the right, n = 10, the formula for the central rectangles
ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.
TECHNICAL SCIENCE. 2021. No 1
Многомерное линейное ИУ с полиномиальным ядром. Исходные данные модельной задачи:
К ( х, у )=х1 •••хту1 - ут,
/ (х)= х1 — хт + g (х1 — хт )2; g = 10
Область интегрирования - т-мерный куб:
От ={0 < x1,•.., хт, у^- ут < 1}. Точное решение
ы(х) = С0х1 • • • хт + gx1■■■xm (х1 —хт + С1);
со = 0,5, С1 = (3/4)т, X = -3т, ц = 1.
Результаты трех последовательных вычислений решения при размерности области т = 10 и числе узлов п = 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).
Погрешность решения в норме А около 11 %. На двух последующих вычислениях решения погрешность много меньше: около 1 и 0,4 %.
Многомерная нелинейная задача - ИУ Гаммерштейна. Рассмотрим нелинейное ИУ Фредгольма вида [10]:
цы ( х)-А | К ( х, у) Ь (ы (у ), у ) йу = / (х ), х еК, (4)
V
соответствующая ИУ (4) квадратурная формула:
А "К (у', у}),, Л -цы, , )'Ь(и],у) = /, ...
=1 р( у]) ] (5)
ы' ~ ы(у, ), ' = 1, — , п.
В качестве тестовой задачи рассмотрим ИУ:
цы ( х)- АIК (х, у ) ы 2 (у)йу = / ( х), х е От, (6)
V
К(х,у) =х1 —хту1 —ут, /(х) = х1 • • • хт• Имеется два точных решения:
[±41 - А4- т+1
( x) = C{k) xi ... xm; c(l,2) = -
ы (х
2•4-тА
и(х) = е(кх ...Хт; X = -г3т, г = 1 ... 100, ц = 1.
Результаты решения ИУ (6) при помощи (5) и метода Монте-Карло при следующих ис-
ходных данных: m = 10, r = 100; начальное приближение м(0 = f (y), i — 1,..., n .
Набор координат узлов случайной сетки моделировался для плотности распределения:
p (х) — 4т (х1.. ,xm )3 , x е Dm, поэтому значения
координат находились по соответствующей
формуле метода обратных функций: yik — (y'k )1/4,
у'к — случайное число для соответствующих узла и его координаты. Значения координат представлены в табл. 1. Результаты итераций по методу Ньютона даны в табл. 2 и 3.
Таблица 1 / Table 1
Координаты узлов псевдослучайной сетки / Coordinates of pseudo-random grid nodes
Номер
узла
1 2
3
4
5
6
7
8
9
10
Координаты узлов
0,92 0,86 0,73
0,89 0,80 0,72
0,73 0,73 0,64
0,95 0,93 0,98
1,00 0,84 0,84
0,67 0,69 0,87
0,59 0,80 0,60
0,48 0,83 0,93
0,72 0,91 0,76
0,73 0,88 0,33
0,94 0,78
0,60 0,65
0,91 0,83
0,72 0,53
0,94 0,99
0,85 0,45
0,40 0,85
0,88 0,74
0,61 0,98
0,89 0,95
0,94 0,57
0,71 0,94
0,75 0,90
0,90 0,67
0,95 0,49
0,80 0,79
0,97 0,83
0,94 0,94
0,88 0,63
0,87 0,78
0,80 0,63
0,81 0,86
0,91 0,73
0,75 0,91
0,99 0,78
0,88 0,75
0,28 0,47
0,79 0,93
0,81 0,69
0,83 0,71
0,77 0,73 0,60 0,98 0,62 0,45 0,84 0,90 0,90 0,51
Таблица 2 / Table 2
Результаты для первой итерации / Results for the first iteration
Номер узла Погрешности итераций Приближенное решение Точное решение
1 -0,04038 0,04755 0,03006
2 -0,03057 0,03600 0,02276
3 -0,03223 0,03796 0,02400
4 -0,06108 0,07193 0,04547
5 -0,06733 0,07928 0,05012
6 -0,01310 0,01543 0,00975
7 -0,00409 0,00481 0,00304
8 -0,06549 0,07712 0,04875
9 -0,03773 0,04443 0,02809
10 -0,01630 0,01919 0,01213
После первой итерации погрешность решения составила 58 % для среднеарифметической нормы. После четвертой итерации погрешность оказалась меньше 0,01 %.
ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.
TECHNICAL SCIENCE. 2021. No 1
Таблица 3 / Table 3
Результаты для четвертой итерации / Results for the fourth iteration
Номер узла Погрешности итераций Приближенное решение Точное решение
1 -0,00009 0,03006 0,03006
2 -0,00007 0,02276 0,02276
3 -0,00007 0,02400 0,02400
4 -0,00014 0,04547 0,04547
5 -0,00016 0,05012 0,05012
6 -0,00003 0,00975 0,00975
7 -0,00001 0,00304 0,00304
8 -0,00015 0,04875 0,04875
9 -0,00009 0,02809 0,02809
10 -0,00004 0,01213 0,01213
Соответствующий выбор начального приближения (например, м(0) — -50 / (у- ), I — 1,..., п)
позволяет найти второе решение данного ИУ. Результаты вычислений аналогичны рассмотренным и потому не приводятся.
Многомерные ИУ первого рода и их решение с использованием регуляризации. Решение методами Монте-Карло и квази-Монте-Карло подобных задач в изученной литературе не описано [1 - 16]. Исходные данные модельной задачи:
К ( х у )=Х1... хту1... уп,
/(х) = х1...хп + я(х1...хп)2; я=1; область интегрирования
оп ={о < х1,..., хп,У1,...,Уп < 1}, п —10.
Рассматриваемая задача относится к некорректным, что является типичным для ИУ Фредгольма первого рода. Нелинейное слагаемое в Дх) играет роль возмущения. Для решения использован известный метод регуляризации Лаврентьева, в соответствии с которым ИУ преобразовывалось в ИУ второго рода с малым коэффициентом при функции решения ц=0,01; точное решение ИУ приобретает вид
и(х) = СоХ1 ... Хт + g(Xl ... Хт)2/ц; его упрощенное решение:
1 + ^ • 4-п /ц
uq(x) = coxl ...хт, c0 = -
X = -3m.
ц - X•3-п
Набор координат узлов случайной сетки моделировался для плотности распределения: р(х) — ст(х1 ...хт)с 1, хеОт, значения координат находились по формуле метода обрат-
у'к - случайное число
ных функций: Yk = (Yk f °,
для соответствующих узла и его координаты в методе Монте-Карло, либо число последовательности Хальтона H(2) в методе квази-Монте-Карло. Результаты вычислений при n = 100 для 10 испытаний для каждого метода, даны в табл. 4, 5. Параметр «степень c» является характеристикой плотности распределения узлов интегрирования.
Таблица 4 / Table 4
Результаты для псевдослучайной последовательности / The results for the pseudorandom sequence
Упрощенное решение Численное решение
Степень Средняя Стандарт Средняя Стандарт
с погреш- погрешности, погреш- погрешности,
ность, % % ность, % %
1 130 3 125 47
2 111 10 52 33
3 93 1 6 3
4 93 1 3 3
5 94 0,4 6 4
Таблица 5 / Table 5
Результаты для последовательности Хальтона Н(2) / Results for the Halton sequence Hi(2)
Степень с Упрощенное решение Численное решение
Средняя погрешность, % Стандарт погрешности, % Средняя погрешность, % Стандарт погрешности, %
1 128 3 166 64
2 106 5 26 16
3 92 2 9 6
4 92 1 3 2
5 94 1 6 4
Данные табл. 4, 5 показывают, что точность обоих методов приблизительно одинакова. Плотность распределения узлов подбиралась в соответствии с видом упрощенного решения, получающегося, если параметр g положить равным нулю. Степень с = 4 является оптимальной, это объясняется тем, что в подынтегральной функции основной вклад при интегрировании вносят кубические члены.
Проблема собственных значений. Одним из преимуществ данного метода является возможность численного решения полной проблемы собственных значений для многомерного интегрального оператора. Отыскание приближенных собственных значений и собственных функций осуществляется на основе решения полной проблемы собственных значений для соответствующей матрицы квадратурной формулы (3). Согласно известной теории предельное множество собственных значений аппроксимирующих интегральный оператор матриц включает в себя множество собственных значений этого инте-
ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.
TECHNICAL SCIENCE. 2021. No 1
грального оператора. Существует множество методов решения полной проблемы собственных значений для матриц: методы Якоби (вращений), Крылова, Леверье-Фаддеева, Данилевского и ряд других. В результате исследования установлено, что одним из лучших из рассмотренных методов является метод Леверье-Фаддеева: метод Якоби предполагает симметричную матрицу, метод Крылова экономичен, но его область применимости ограничена (для рассматриваемых задач не удалось найти решение данным методом), метод Данилевского относительно сложен. Результаты решения задачи для различного числа квадратурных узлов представлены в табл. 6. Они показывают, что, несмотря на небольшую точность аппроксимации функции - решения ИУ, точность аппроксимации собственных значений удовлетворительная даже для небольших п.
В качестве модельной рассмотрена задача (1) при
K (х, y) = sin (5х + 5 y), f (х) = х + 5; 0 < х < 1; А = 1,5, ц = 1.
Ядро ИУ вырожденное, точное решение имеет вид u(х) = C1cos(5х) + C2sin(5х) + С0 (х),
где С\, С2 - некоторые константы; Co(x) - один из нулей интегрального оператора (отвечает собственному значению zo = \). ИУ редуцируется к СЛАУ относительно констант С\, С2 с матрицей
лп
= 1 -15 (1 - cos 20 v
f
a 21 = 1,5
0,5 + -
sin
(10)), (10)
(
-<12
= -1,5
0,5 -■
sin (10) 20
A
20
22 =1 - M (1 - cos (10)) •
N Метод 1 Метод 2 Метод 3
Z1 Z2 Z1 Z2 Z1 Z2
5 1,58452 0,08765 1,75172 0,25856 1,50610 0,04589
10 1,60494 0,10736 1,51632 0,04268 1,63468 0,13484
15 1,60833 0,11063 1,78386 0,28386 1,59268 0,09551
20 1,60950 0,11176 1,55148 0,07998 1,56529 0,10672
Наилучшие результаты показал метод прямоугольников, несколько хуже точность у метода квази-Монте-Карло.
Результаты вычислений методом Монте-Карло менее надежны. В данном примере обнаружилось интересное свойство: сходимость для собственных значений много лучше (в сотни раз), чем для функции решения ИУ.
Сингулярные ИУ. Метод выделения особенности. Рассмотрим ИУ Фредгольма (1) с сингулярным ядром К(х,у), неограниченным при х = у. Для решения применим метод выделения особой точки:
цм (У)- X j K (У, у) u (y) dy -
Vi
-X j K(y1,y)u (y)dy = f (y1), i = 1,..., n.
(7)
Коэффициенты характеристического многочлена данной матрицы det (A — zE) — z2 — pi z — p2
равны р1 = 1,72414, р2 = - 0,18233; собственные значения (нули многочлена) равны z\ = 1,61096, Z2 = 0,11318; результаты, полученные тремя методами, даны в табл. 6 (метод 1 - метод прямоугольников, метод 2 - Монте-Карло, метод 3 -квази-Монте-Карло).
Таблица 6 / Table 6
Результаты вычисления собственных значений / Results of calculating eigenvalues
Ядро первого интеграла в (7) является неограниченным. Размеры области Vявляются относительно малой величиной порядка
О (п_1/т ), при п ^ да, поэтому вариация и(у) в
этой области стремится к нулю и интеграл может быть вычислен приближенно по формуле
{ К(у', у)ы (у)йу К (у, у)йуы (у' ) = Згы (у') ,
у'' еV .
Точку у, целесообразно размещать в центре области V-.
Ядро второго интеграла в (7) является регулярным, поэтому интеграл может быть вычислен при помощи метода Монте-Карло:
| К (у', у) ы (у) йу «—^ X $] (у') ,
Í л K(yi,У' ) / л
где Sj (y)= \ У u (yJ),
iy)
п - количество точек
из выборкИ, у1 =(у{, —, у1т ), уп =(у" у"т ),
принадлежащих замкнутой области V , плотность р(у) нормирована на единицу в области V \ .
В результате получаем из (7) соответствующую СЛАУ для нахождения приближенных значений решения ыi « ы (у1) в узлах случайной сетки:
(ц - %Ji)м -i = 1, ..., n.
V KÍMu = f,
n - ni jeVV p (У ) j ' (8)
ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.
TECHNICAL SCIENCE. 2021. No 1
Сингулярное ИУ с ядром типа Коши.
1
ци (x)-Xj( x - y) u (y) dy — x, x e( 0,1), (9) 0
Vi — jx e(0,1): |x -y,-| < kh /2}, i —1,...,n; h — —,
где с, к - некоторые положительные константы, X = 0,2, ц = 1.
Для решения (9) на основе (8) были рассмотрены:
1) метод Монте-Карло. Случайные равномерно распределенные на интервале (0,1) узлы (также модифицированное распределение с отсевом чрезмерно близких точек);
2) метод квази-Монте-Карло. Узлы сетки -точки низкодисперсной последовательности Хальтона И(), I = 1,...,и; 5 - простое число;
3) методы центральных прямоугольников и трапеций.
На рис. 3, а даны графики численного решения ИУ (9) при псевдослучайном равномерном распределении точек интегрирования.
и(х) 2,5
u(X 2,5
2,0
1,5
1,0
0,5
0
ЩУ
0,2 0,4 0,6 0,8 1,0
и(х) 2,5
2,0
1,5
1,0
0,5
0
0,2 0,4 0,6 0,8 1,0
Рис. 3. Приближенное решение ИУ (9): а - сплошная линия -метод центральных прямоугольников (п = 100), кружки -метод Монте-Карло для 100 псевдослучайных точек (слева к = 1; справа к = 5); б - сплошная линия - метод центральных прямоугольников (п = 100), кружки - метод
квази-Монте-Карло для 100 точек последовательности Хальтона (слева к = 1; справа к = 5) / Fig. 3. An approximate solution of the integral equations (9): а - solid line - the Central method of rectangles (n = 100), circles - Monte-Carlo method
for 100 pseudo-random points (left к = 1; right к = 5); б - solid line - the Central method of rectangles (n = 100), cups - method quasi-Monte-Carlo for 100 points in the sequence of Halton (left к = 1; right к = 5)
При чрезмерном сближении псевдослучайных точек точность, как правило, уменьшается, а число обусловленности матрицы СЛАУ увеличивается. Для большей равномерности распределения последовательность точек была модифицирована посредством исключения чрезмерно близких точек.
На рис. 3, б даны графики численного решения ИУ (9) методом квази-Монте-Карло. В качестве узлов сетки были выбраны точки последовательностей Хальтона [1].
Результаты вычислений показывают, что увеличение в некоторых пределах размера элемента разбиения, содержащего особенность ядра, приводит к улучшению решения. Однако дальнейшее увеличение размера рассматриваемого элемента приводит к снижению точности из-за погрешности аппроксимации решения в этом элементе.
Метод на основе параметрической регуляризации ядра. В уравнении (9) сингулярное ядро К(х, у) = (х - у)-1 заменяется на асимптотически близкое регулярное выражение
к (*, у) = (* - у) / (* - y)2 +
•0.
Рассмотрим результаты численного решения (9) при X = 1, ц = 1, представленные на рис. 4.
Регуляризация ядра и модификация псевдослучайной равномерной последовательности за счет исключения чрезмерно близких точек приводит к значительному улучшению результата. Недостатком метода является сильная зависимость от значений малого параметра в.
Интерполяционно-проекционный метод. В этом методе приближенное решение отыскивается в виде суммы
п I \ и (х)«£С/ф (х - у]), (10)
М
где ф(х) - базисная функция, например, широко употребительная финитная кусочно-линейная функция ф(х) — тах(о,1 - 2|х| / к), Н я —, С - неиз-V II/ п
вестные коэффициенты.
С учетом (10) из (1) следует соотношение
ц £ Сj ц (* - у j) « j=1
n 1
! X £ С j j K (*, у )ц( у - у j )dy + f (*), 0 < * < 1. j=1 о
x
а
ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.
TECHNICAL SCIENCE. 2021. No 1
Из (11) на основе интерполяционного метода находим СЛАУ для вычисления неизвестных коэффициентов выражения (10):
n i ^ ZCjф (у, - У j=i
j )=
(12)
n
^ zc Jk (y, y) ф (y - yj) dy + f (yi), j=1 0
i = 1,., n.
u(x
1,5 1,2 0,9 0,6 0,3
0
u(x)
1,5 1,2
0,9 0,6 0,3
^ОД-—" 0,6 0,8 1,0 .r 0
0,2 0,4 0,6 0,8 1,0
u(x)
1,5 1,2
0,9 0,6 0,3
0
u(x)
1,5 1,2
0,9 0,6 0,3
0,2 0,4 0,6 0,8 1,0
* 0
б
J
0,2 0,4 0,6 0,8 1,0
(u, > =
1
2l +1
l ( jh Z Ln \ у, +-T; ul,.,un
j=-i
l
Соответствующий эффект иллюстрируется на рис. 5 при l = 10.
u(x)
0,8 0,6 0,4 0,2
0
u(x) 1,0 0,8 0,6 0,4 0,2
0,2 0,4 0,6 0,8 1,0 0
'0,4 0,6 0,8 1,0
Рис. 5. Решение ИУ (9) при n = 100: сплошная линия -метод центральных прямоугольников (со сглаживанием по двум точкам), кружки - метод квази-Монте-Карло для последовательности Хальтона (слева со сглаживанием) / Fig. 5. Solution of the integral equations (9) for n = 100: solid line - the method of Central rectangles (with two point smoothing), circles-the quasi-Monte-Carlo method for the Halton sequence (on the left with smoothing)
Адаптация. В [4] предложен метод вычислений с адаптивной плотностью, который основан на соотношении
Р1
(у' )'
cu,■
(yj, У)
n K2
j=j p(yj)
i = 1,
n; c = const,
Рис. 4. Приближенное решение ИУ (9): а - сплошная линия -метод центральных прямоугольников (п = 100), кружки - метод Монте-Карло для 10 псевдослучайных точек (слева s — ^/ОД , справа 8 = 0,1); б -при 8 = 0,1: сплошная линия - метод центральных прямоугольников (n = 100), кружки - метод Монте-Карло для 100 псевдослучайных точек (справа последовательность точек модифицирована) / Fig. 4. Approximate solution of the integral equations (9): а - solid line - the method of Central rectangles (n = 10), circles - the Monte-Carlo method for 10 pseudo-random points (left s — , right s = 0,1); б - for e = 0,1: solid line - the method of Central rectangles (n = 100), circles - the Monte-Carlo method for 100 pseudo-random points (on the right, the sequence of points is modified)
Построенный метод может быть все же недостаточным по точности при решении ИУ с сингулярным ядром даже при большом количестве узлов интегрирования (рис. 5, справа). Уточнение может быть достигнуто различными способами, например, на основе соотношения (11) с применением метода Галеркина. Однако более экономичным является применение сглаживания, например, по схеме
где р у) - начальная плотность распределения случайных точек; р1 (у1) - адаптивная («среднеквадра-тическая») плотность распределения; с - нормировочный множитель; и- - сеточное решение на первом этапе вычислений.
Могут быть рассмотрены также другие варианты адаптивной плотности, более подходящие (по причине сингулярности ядра) для рассмотренного примера:
Р2
Рз
(У > (у' )'
C\u,
Z
K
(yj, У)
n -1 j=1,j*i
Z -
j=1,j *i
n-1
p(yj)
(yj, У)
---i
p(yj)
i = 1,.,n; c = const;
,i = 1,.,n; c = const.
Отмеченные усредненные адаптивные плотности преимущественно дают эффект лишь для ИУ с гладким ограниченным ядром. Для сингулярных ИУ построение эффективной адаптивной плотности является еще не полностью решенной задачей.
Примеры. Решение пространственной краевой задачи (КЗ) Дирихле для уравнения Лапласа. Подобные КЗ описывают стационарные и квазистационарные температурное, элек-
X
X
а
х
n
C
ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION.
TECHNICAL SCIENCE. 2021. No 1
трическое и магнитное поля, часто решаются в различных приложениях и имеют вид
д2и д2и д2и _ / ч
—г + —г + —т — 0, (х,у,^)еГ;
дх2 ду2 дг2 ^ } (12)
: = Us , (X, y, z )e S,
где Ыв - значение искомой функции на границе области. Решение КЗ (12) отыскивается в виде потенциала двойного слоя (ПДС):
/ , Л
u (M ) = -J w (P)
д
S
диг
1
V 4nrMP
dSP, M eV,
плотность которого отыскивается из интегрального уравнения
Us(M ) = 1 w(M )-J w( P )Л
í
1
Л
Таблица 7 / Table 7
Решение для шара / The results of calculation for ball field
n U1 51, % U2 52, %
100 0,577 44,92 -0,080 1794,7
200 0,447 12,18 -0,050 1157,0
300 0,377 5,22 0,007 48,00
400 0,355 10,71 0,013 170,1
500 0,362 8,95 0,001 75,30
600 0,424 6,53 0,008 58,79
700 0,447 12,20 0,003 32,31
800 0,423 6,23 0,001 81,73
900 0,411 3,35 0,009 101,1
1000 0,394 0,90 0,005 1,34
dSP, M e S.
S wnP V 4nrMP )
Интегральное уравнение решается на основе метода, описанного выше в подразделе «Метод выделения особенности».
Область - шар. На поверхности шара единичного радиуса задано распределение гармонической функции:
us (О, ф) = 0,15sin(20)cos(ф) + +1,5(1 - cos (20)) sin (2ф),
где 0 - сферический, а ф - полярный угол, О e (0,2п).
Ядро ПДС:
К (гм ,0м ,Фм, rP, 0P ,ф p ) =
= nx (XP - XM ) + ny (yp - Ум ) + nz (ZP - zM ) ¡ 2 2 2 \3/2
4п((xp - хм) +(Ур - Ум) +(zp - zM) ) xR2sin (0P),
xP = Rsin (0P) cos (фр), yP = Rsin (0P) sin (фр), zP = Rcos (0 P),
Хм = rsin (0) cos ( ф ), Ум = rsin (0) sin ( ф ), ZM = rcos (0),
nx = sin (0P)cos (фр ), ny = sin (0P)sin (фр ),
nz = cos (0 p),
где r - сферический радиус, re(0, R).
В табл. 7 даны результаты вычислений решения в точках с координатами ri = 0,9, Г2 = 0,1, 0i,2 = л/8, ф1,2 = л/6, при различных числах узлов интегрирования n, определяемых последовательностью Хальтона.
Область - тор. На поверхности тора с единичным средним радиусом продольного сечения задано распределение гармонической функции:
us (с,ф) = sin (с)cos то)-cos (с),
где (с, х, ф) - тороидальные координаты, связанные с декартовыми координатами соотношениями:
a • shx x — —-cos ф,
У:
chx - cos a
a • shx chx - cos a a • sina
sin ф,
z = -
chx - cos a
-п < о < п, 0 < t < да, 0 < ф < 2п.
Размеры тора определяются параметрами а -1, х0 = 1. Ядро ПДС:
К (rM, тм ,фм,о р ,тр ,фр ) =
_ [nx (xP - Хм) + ny (yp - Ум) + nz (ZP - zм)]
/ 2 2 2\3/2
((Хр - Хм ) +(Ур - Ум ) +(zP - zM ) )
X p02sh ( х0 )
Po =
ch (т0 )- cos (с) '
nx — -p0-1dxP / dx,
ny —-Po-1dyp / dx,
nz — -po~ldzp / dx.
В табл. 8 даны результаты вычислений решения в точках с координатами Х1 = 2, Х2 = 3; С1,2 = л/4, ф1,2 = л/6, при различных числах узлов интегрирования п.
a
ISSN 1560-3644 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION. TECHNICAL SCIENCE. 2021. No 1
Таблица 8 / Table 8
Решеие для тора / The results of the calculation for the toroidal field
n U1 61, % U2 62, %
100 0,287 34,41 0,123 50,09
200 0,274 28,26 0,112 36,62
300 0,230 7,93 0,083 1,44
400 0,232 8,58 0,091 10,59
500 0,237 11,11 0,094 14,77
600 0,234 9,38 0,088 7,08
700 0,230 7,78 0,087 5,58
800 0,221 3,39 0,085 3,22
900 0,218 2,25 0,083 1,00
1000 0,224 5,04 0,089 8,24
Заключение
В статье исследован статистический метод решения интегральных уравнений Фредгольма большой размерности как с гладким, так и с сингулярным ядром. Рассматриваемые подходы позволяют расширить круг задач теории интегральных уравнений, решаемых методами Монте-Карло и квази-Монте-Карло, поскольку отсутствуют ограничения на величину нормы интегрального оператора. Рассмотрена серия примеров, демонстрирующих степень эффективности исследуемого метода
Выражаем благодарность профессору С.А. Некрасову за помощь в подготовке статьи.
Литература
1. Ермаков С.М. Метод Монте-Карло в вычислительной математике (вводный курс). СПб.: Изд-во Бином. 2011. 192 с.
2. Иванов В.М., Кульчицкий О.Ю., Кореневский М.Л. Комбинированный метод решения интегральных уравнений
// Дифференциальные уравнения и процессы управления. № 1. 1998. С. 1 - 40.
3. Сипин А.С. Бессеточные методы Монте-Карло решения краевых задач для уравнений в частных производных: автореф. дис. ... д-ра физ.-мат. наук. СПб., 2016. 32 с.
4. Иванов В.М., Кульчицкий О.Ю. Метод численного решения интегральных уравнений на случайной сетке // Диф. уравнения. 1990. Т. 26, № 2. С. 333 - 341.
5. Берковский Н.А. Модернизация полустатистического метода численного решения интегральных уравнений: автореф. дис. ... канд. физ.-мат. наук. СПб., 2006. 15 с.
6. Кореневский М.Л. Разработка адаптивно-статистических методов вычисления определенных интегралов: дис. ... канд. физ.-мат. наук. СПб., 2000. 161 с.
7. Скорикова О.В. Сильная равномерная распределенность системы функций Ван дер Корпута-Хеммерсли // Изв. Тульского гос. ун-та. Естественные науки. 2013. Вып. 3. С. 91 - 102.
8. Антонов А.А., Ермаков С.М. Эмпирическая оценка погрешности интегрирования методом квази-Монте-Карло // Вестн. Санкт-Петербургского ун-та. Математика. Механика. Астрономия. 2014. Т.1, вып. 1. С. 3 - 11.
9. Антонов А.А. Алгоритм численного интегрирования
методом квази Монте-Карло с апостериорной оценкой погрешности // Вестн. Санкт-Петербургского ун-та. Математика. Механика. Астрономия. 2015. Т.2, вып. 1. С. 3 - 11.
10. Корн Г., Корн Т. Справочник по математике (для научных работников и инженеров). М.: Наука, 1973. 832 с.
11. Стоянцев В. Т. Решение задачи Коши для параболического уравнения методом квази-Монте-Карло // Журн. вычисл. математики и математической физики. 1973. Т.13, № 5. С. 1153 - 1160.
12. Стоянцев В.Т. Решение задачи Дирихле методом квази-Монте-Карло // УМН. 1975. Т.30, № 1(181). С. 263 - 264.
13. Некрасов С.А., Ткачев А.Н. Теория вероятностей и ее приложения: учеб. пособие. Новочеркасск: ЮРГТУ, 2007. 148 с.
14. Некрасов С.А. Методы ускоренного статистического моделирования и их применение в электротехнических задачах / Изв. вузов. Электромеханика. 2008. № 5. С. 13 - 19.
15. Некрасов С.А. Решение интегральных уравнений методом Монте-Карло / NovaInfo.Ru (эл. журн.). 2016. № 53.
16. Некрасов С.А. Решение п-мерного уравнения Шредин-гера методом интегральных уравнений на псевдослучайной сетке // NovaInfo.Ru (эл. журн.). 2016. № 55.
References
1. Ermakov S.M. The Monte-Carlo Method in computational mathematics (introductory course). Saint Petersburg: Binom Publishing House. 2011. 192 p.
2. Ivanov V.M., Kulchitsky O.Yu., Korenevsky M.L. Combined method for solving integral equations / Differential equations and control processes. No. 1. 1998. Pp. 1 - 40.
3. Sipin A.S. Grid-Free Monte-Carlo methods for solving boundary value problems for partial differential equations. The author's abstract Diss. for the degree of doctor of physics.-math. sciences'. Saint Petersburg, 2016. 32 p.
4. Ivanov V.M., Kulchitsky O.Yu. Method of numerical solution of integral equations on a random grid / Differents. Equations. Vol. 26, No. 2. 1990. Pp. 333 - 341.
5. Berkovsky N.A. Modernization of the semi-statistical method for numerical solution of integral equations: the author's abstract diss. for the degree of candidate of physical sciences. - math. sciences'. Saint Petersburg. 2006. 15 p.
6. Korenevsky M.L. Development of adaptive statistical methods for calculating certain integrals: diss. for the degree of candidate of physical sciences. - math. sciences'. Saint Petersburg, 2000. 161 p.
ISSN 0321-2653 IZVESTIYA VUZOV. SEVERO-KAVKAZSKIYREGION. TECHNICAL SCIENCE. 2021. No 1
7. Skorikova O.V. Strong uniform distribution of the system of van der Korput-Hemmersley functions // Proceedings of the Tula state University. Natural science. 2013. Issue 3. Pp. 91 - 102.
8. Antonov A.A., Ermakov S.M. Empirical estimation of integration error by the quasi-Monte-Carlo method // Vestnik of the St. Petersburg University: Mathematics. 2014. Issue 1. Vol. 1. Pp. 3 - 11.
9. Antonov A.A. Algorithm for numerical integration by the quasi-Monte-Carlo method with a posteriori error estimation / Vestnik of the St. Petersburg University: Mathematics. Issue 1. Vol. 2. 2015. Pp. 3 - 11.
10. Korn G., Korn T. Handbook of mathematics of engineers. Moscow: Nauka. 1973. 832 p.
11. Stoyantsev V.T. Solution of the Cauchy problem for a parabolic equation by the quasi-Monte-Carlo // Computational Mathematics and Mathematical Physics. Vol. 13. No. 5. 1973. Pp. 1153 - 1160.
12. Stoyantsev V.T. Solution of the Dirichlet problem by the quasi-Monte-Carlo method / UMN. Vol. 30. No. 1 (181). 1975. Pp. 263 - 264.
13. Nekrasov S.A., Tkachev A.N. Probability Theory and its applications: training manual. Novocherkassk: YURSTU. 2007. 148 p.
14. Nekrasov S.A. Methods of accelerated statistical modeling and their application in electrical problems / Scientific and Technical Journal. Russian Electromechanics. 2008. No 5. Pp. 13 - 19.
15. Nekrasov S.A. Solution of integral equations by the Monte Carlo method // NovaInfo.Ru (E-journal). 2016. No. 53.
16. Nekrasov S.A. Solution of the n-dimensional Schrodinger equation by the method of integral equations on a pseudo-random grid // NovaInfo.Ru (E-journal). 2016. No. 55.
Поступила в редакцию /Received 09 декабря 2020 г. /December 09, 2020