Квадратурная формула для сингулярных интегралов с ядром Гильберта

Фундаментальные основы теории сингулярных интегральных уравнений и их многочисленные приложения подробно изложены в классических трудах Н. И. Мусхелишвили [1, 2], И. Н. Векуа [3], Ф. Д. Гахова [4] и других авторов. Задача получения конструктивных решений, доведенных до числового результата, неоднократно определялась ведущими математиками как одно из приоритетных направлений исследований. В частности, важность "доведения до числа" подчеркивали Н. И. Мусхелишвили, М. А. Лаврентьев и И. Н. Векуа. Пионерские результаты в данной области были представлены М. А. Лаврентьевым в работе 1932 года [5]. В предисловии к этой работе академик С. А. Чаплыгин указал на высокую значимость полученных результатов для практики и выразил надежду на дальнейшее развитие исследований в этом направлении. Значение работы М. А. Лаврентьева получило высокую оценку и в трудах Н. И. Мусхелишвили, который во всех изданиях своей монографии [2] неизменно отмечал, что разработка приближенных методов решения сингулярных интегральных уравнений является одной из наиболее актуальных задач современной теории. Существенный прогресс в развитии вычислительных подходов, вызванный потребностями прикладных задач, связан с исследованиями С. М. Белоцерковского (1954 г.) в области аэродинамики. Эти работы заложили основы метода дискретных особенностей - эффективного численного метода решения сингулярных интегральных уравнений, получившего в настоящее время широкое распространение и признание. Отдельного упоминания заслуживает также значительный вклад академика С. Л. Соболева в развитие теории квадратурных и кубатурных формул [6], составивший теоретическую базу для последующего обоснования многих численных алгоритмов. Сингулярные интегральные уравнения являются универсальным языком описания широкого круга физических явлений, характеризующихся наличием особенностей, и служат основой для построения высокоточных численных моделей. Аналитические решения таких задач, как правило, представимы в виде сингулярных интегралов с известными ядрами.

В классической теории сингулярных интегральных уравнений часто используются неконструктивные методы: предельные переходы, принцип исключенного третьего, аппроксимации без оценки погрешности. В рамках конструктивной математики и теории алгоритмов [7] требуется, чтобы все построения были эффективно вычислимы, а сходимость - конструктивно обоснована.

В данной работе предлагается квадратурная формула, удовлетворяющая принципам конструктивного анализа: узлы сетки вычислимы; веса формулы выражаются через элементарные функции; отсутствует деление на величины, которые могут обратиться в нуль в недетерминированных точках; интеграл понимается в смысле главного значения, но с явной вычислимой оценкой.

Отметим, что в работе [8] получена следующая квадратурная формула:

I'

где

^л P(t,%)

, sin(t — %)

I t^l ¹ -(—y ' cos'.-cos^l,   (1)

2и + 1    sin( tk — %) [                      cos % J як tk = о—Тк = = 0,1, -,2п. 2п + 1

Следует отметить, что соотношение (1) для некоторых значений аргумента % приводит к потере точности из-за деления на cos % в знаменателе. Вторым недостатком квадратурной формулы (1) является особенность, возникающая при % = я/2, где cos % = 0 и cos(2n + 1)% также равен нулю. Указанные особенности актуализируют получение новых квадратурных формул, позволяющих исключить эти ситуации.

Следует отметить также квадратурную формулу [9] интерполяционного типа, не содержащую деления на cos %:

т-1

Г^п p(t,%)       я v^-1            t_k — %

I ^Пр— 0^dt “m ^Ip(t ‘ ' ^X)ctg— ' ^m = ²ⁿ ' ⁰                      к=0

Данная формула является приближенной, о точности которой автор в своей работе ничего не конкретизирует, кроме того, что в пределе она стремится к нулю. Это информация не применима для практического применения и не позволяет определить значение и в квадратурной формуле для заданной точности в расчетах. В формуле веса содержат котангенс, что при tk ~ % требует отдельной аналитической обработки особенности.

Конструктивная постановка задачи и метод решения

Рассмотрим некоторую непрерывную по своим аргументам функцию F(t, %) на квадрате [0, я] х [0, я], за исключением t = %, где она имеет сингулярные особенности вида sin(t — %). Пусть F(t, %) является я — периодической по каждой переменной, т. е.

F(я + t, %) = F(t, %), F(t, я + %) = F(t, %). (2)

Требуется вычислить (в смысле вычислимой аппроксимации) значение интеграла:

/(%) = I F(t, %)dt.

В силу периодичности функции F(t, %) имеем следующее представление: p ^(t, %) F ⁽ t,%) = —71 ---v sin(t — %)

где p(t, %) — вычислимая функция на всем квадрате [0, я] х [0, я] функция, удовлетворя- ющая условиям:

р(2я + t,%) =—p(t,%), p(t, 2я + %) =—p(t, %).

Поэтому интеграл (3) можно представить в виде:

^п p(t, %)

/(%)=!^-dt.

0 sin(t — %)

Замечание. Ядро sln(t — %)

называется ядром Гильберта для периодических функций на отрезке [0, п]. Оно имеет особенность первого порядка при t = х. Интеграл понимается как главное значение по Коши, но в конструктивной формулировке мы заменяем его конечной суммой с оценкой погрешности.

При фиксированном значении х функцию p(t, х) будем рассматривать как функцию по переменной t на интервале [0, т] . В силу периодичности ядра интегрального оператора удобно использовать тригонометрические полиномы.

Обозначим ^(t) = p(t, х). Введем N = 2n + 1 равноотстоящих узлов на периоде: пк

tk = ;-----т, к = 0,1,..., 2п. ^k 2п +1,      ,

Построим тригонометрический интерполяционный полином ^_n(t), интерполирующий функцию ^(t) в узлах t k . Для нечетного числа узлов на периоде [0, п] ядро интерполяционного процесса может быть представлено в виде

2п

Фп(Г) = ^ФШ ■ ik(t), к=0

где l k (t) — фундаментальные полиномы Лагранжа, соответствующие тригонометрической интерполяции. При равномерном распределении узлов с шагом h = ~~ указанные базисные функции принимают вид:

I (t) =------k     2n + 1

sin ((п + 2) ⁽ t - t k )) s ⁱⁿ (2 ⁽ t — t k ))

,

что представляет собой классическое ядро Дирихле. Однако для последующего интегрирования с ядром

1/sin(t — х)

более удобной оказывается эквивалентная форма записи, использующая произведение синусов. Применяя тригонометрические преобразования, интерполяционную формулу можно привести к следующему виду:

2п

2т у

P ^(t,x) ~ 2nЧ^^P^(t k ^,X) ’ к=0

sin((n + 1)(t_k — t))sin(n(t_k — t)) sin(2(t k — t))           .

Подставим аппроксимацию (8) в исходный интеграл (5):

Я         2n

Ч (Ч-llP^x)- 0          k=0

sin((n + 1)(t_k — t))sin(n(t_k — t)) sin(2(t k — t))

)

dt sin(t — x)'

Меняем местами знаки суммы и интеграла и получаем:

2 ^{п                Я} Si ⁿ⁽⁽ⁿ + ^1)(t k ^{— t))} Si ⁿ⁽ⁿ⁽ t k ^{— t))}

2n + 12j^{P k,X} J_o     sin(2(t_k — t)) • sin(t — х)        ^.

Из [10] известно, что тригонометрические полиномы обладают свойством "воспроизведения" для сингулярных интегралов с ядром Гильберта.

Лемма (Основное тождество) [10]. Для любых t k ^ х справедливо: ^resin((n + 1)(t k — t))sin(n(t k — t))

I ------—-----7----;----;---dt sin(2(tk — t)) • sin(t — х)

sin((n + 1)(tk — x))sin(n(tk — х))

sin(2(tk — %))

.

При t * = х правая часть понимается как предел:

sln((n + 1)u)sln(nu)

lim------7ТТ7----- u^o      si:(2u)

= 0,

где и = t_k — х. В этом случае левая часть интеграла становится регулярной, и тожде-

ство сохраняется.

Обоснование предела:

sln((n + 1)и) ~ (n + 1)и, sln(nu) ~ пи, sln(2u) ~ 2и(и ^ 0),

откуда

((п + 1)и)(пи) п(п + 1)

------------=-------и ^ 0.

2и

Подставляя результат леммы (10) в выражение (9), получаем искомую квадратурную

формулу:

[ ^                    2п

®(t, х)          2п У'¹         sln((n + 1)(t — х))sln(n(t — х))

sSd^^Ь^-²¹— -х))      (1 ¹ )

*=0

где пк

t * = 2—+Г, к = 0,1, ■", 2п, п — нечетное положительное число.

Численные эксперименты

Для проверки эффективности предложенной квадратурной формулы (11) проведем

ее тестирование на примерах, как допускающих точное аналитическое вычисление, так и не имеющих аналитического решения в элементарных функциях.

Задача 1. Для проверки эффективности предложенной квадратурной формулы (11)

проведем ее тестирование на примере, допускающем точное аналитическое вычисление.

Рассмотрим интеграл (5) с функцией

^(С,х) = sln(t — х).

Подставляя (12) в (5), получаем я sln( t — х)

1(х) = I             dt = I dt = и.

0 sin(t — х)

Таким образом, точное значение интеграла не зависит от х и равно п.

Вычислим данный интеграл с использованием квадратурной формулы (11). По

условию

^(t * ,х) = sln(t * — х).

Подстановка последнего в (11) дает

/(х)

2п

2п у 2n__^ _S ^ln ₍ ^t * ^-X ₎

*=0

sln((n + 1)(t * — х)) sln(n(t * — х))

*

—

.

Обозначим

2п

^ * =^ sln(t * — х)

sln((n + 1)(t * — х)) sln(n(t * — х))

.

*=0

*

—

Для демонстрации работоспособности формулы выберем наименьшее нечетное значение п = 1.

Тогда 2п + 1 = 3, и узлы сетки имеют вид t * = —,к = 0,1, 2.

При п = 1 выражение (15) упрощается:

5 * = sln(t * — х)

sln(2(t * — х))sln(t * — х) sln(2(t * — х))

sln²(t * — х).

Следовательно, квадратурная сумма (15) принимает вид

S = sin2(x) + sin2 (я—’) + sin2 (^—x).

Применяя формулу понижения степени для квадрата синуса из (16), получаем:

S = - — - [cos 2x + cos (2Я —2x) + cos (4^—2’)] •

Сумма косинусов в квадратных скобках последнего выражения равна нулю, что проверяется непосредственно.

Таким образом,

S = - .

Подставляя найденное значение в формулу (14), получаем

2я 3 "’'з    ^я

Следовательно, результат, полученный по квадратурной формуле (11), совпадает с точным значением интеграла (13).

Задача 2. Рассмотрим интеграл (5) с функцией

^(t,x) = sin(t + x).

Требуется вычислить интеграл resin(t + x) I(x) = I              dt.

0 sin(t — x)

Данный интеграл понимается в смысле главного значения по Коши, так как подынтегральная функция имеет особенность при t = x при x £ (0, я).

Аналитическое вычисление. В формуле (18) в числителе подынтегральной функции на основании тригонометрических формул осуществляем преобразование:

sin(t + x) = sin((t — x) + 2x) = sin(t — x)cos 2x + cos(t — x)sin 2x.

Тогда подынтегральная функция принимает вид sin(t + x)                   cos(t — x)

—------ = cos 2x + sin 2x--;----- = cos 2x + sin 2x • ctg(t — x).

sin(t — x)                   sin(t — x)

Следовательно,

7Г                                      7Г

I(x) = I cos 2x dt + sin 2x I ctg(t — x) dt = яcos 2x + sin 2x • /(x),

где

/(x) = I ctg(t — x) dt.

Вычислим J(x). Выполним замену переменной u = t — x ^ Пределы интегрирования: при t = 0: u = —x; при t = я: u = я — Тогда

du = dt.

x.

n-x

/(x) = I ctg u du.

-x

При x £ (0, я) точка u = 0 лежит внутри интервала (—x, я — x),

потому что —x < 0 <

я — x. Последний интеграл является несобственным интегралом 2 рода, для которого по определению:

-         f-~x

J(x) = limlj ctg u du + J ctgudu) =

= lim [In | sinu |] £^0 ⁺

--

:X + lim [In | sinu |] £ ^x

= lim (In | sin(—e) | -ln | sin(—x) |) + lim (In | sin(л — x) | -ln | sine |).

£^0+                                  £^0+

Вычисление последних пределов приводит к результату/(x) = 0. Следовательно, исход ный интеграл (18), имеет значение:

/(x) = лcos2x.                                  (21)

Вычислим значение интеграла (18) с помощью квадратурной формулы (11) для функции

^(t,x) = sin(t + x).

nk

Введем обозначение функции ^(tk, x) = sin(tk + x). При n = 1 и узлах tk = —, к = 0,1,2, квадратурная формула (11) принимает вид

2л

^(x) « — > sin(t k + x)

У sin( t_k + x)sin(t_k — x).

sin(2(tk — x)) sin(tk — x) 2л sin(2(tk —x))         3

На основании тригонометрических формул из последнего следует:

sin(tk + x)sin(tk — x) = - (cos 2x — cos 2tk).

Следовательно,

^; (x) ~ ^т’чУ (cos2x 32

— cos 2 tk) =

k=0 2n

С учетом t k = 0, П , — из последнего имеем

^{л 2}

— I 3cos 2x — у cos 2tk .

X           k=0        /

У cos 2t k = 1 —  —  = 0.

k=0

л

I(x) ~  - 3 cos 2x = л cos 2x.

Что совпадает с точным аналитическим значением (21).

Результаты численного интегрирования задачи 2. Проверяем квадратурную формулу (11) для данной задачи. Вычисления проведены с помощью программы, реализованной на JavaScript для функции

^(t, x) = sin(t + x)

при различных значениях аргумента x . В качестве параметра дискретизации выбрано n = 5, что соответствует количеству узлов сетки равное 2n + 1 = 11 . В таблице 1 представлены значения интеграла, полученные по аналитическому выражению (15) и по квадратурной формуле (11), а также абсолютная погрешность:

△ 5 =1 ^(x) — 1 5

Таблица 1. Численные расчеты интеграла по аналитическому выражению и квадратурной формуле

X	/(X)	* 5 (х)	△ 5
0	3.1415926536	3.1415926536	0.0000000000
л/4	0.0000000000	0.0000000000	0.0000000000
1.5	-3.1101531544	-3.1101531544	0.0000000000
2	-2.0534819974	-2.0534819974	0.0000000000
л	3.1415926536	3.1415926536	0.0000000000

Как видно из табл. 1, при n = 5 квадратурная формула (11) дает значение интеграла, совпадающее с аналитическим результатом.

где

Задача 3. Рассмотрим следующий интеграл:

ад = f

ⁿ^₁(t — a)^₄(t + а)

I d 1 (t — x)d₄(t + x) ’

d 1 (t)  = 2q^1/4 ^(-1)ⁿqⁿ⁽ⁿ⁺¹⁾sm(2n + 1)t,

^ 4 (0

n=0

= 1 + 2 У (—1)nqn2 cos 2nt — n=1

тета-функции Якоби [11] с нормой q = e^mT , Im т > 0, | т |= — l^, а и x - заданные па-

раметры.

Замечание. Подынтегральная функция имеет простые полюсы при t = x (в нуле ^ 1 ) и при t = —x (с учетом периодичности). Интеграл понимается в смысле главного значения по Коши.

Получим аналитическое выражение для I(x) в терминах логарифмических производных тета-функций, пригодное для дальнейших численных расчетов.

На основании тождество Вейерштрасса–Якоби [12]:

^ 1 (t — a)^ 4 (t + а) ₌ B 1 (t — x)      ^ 4 (t + x)

B 1 (t — x)d₄(t + x)      ^(t — x)      d₄(t + x) '

из двоякопериодичности левой части по t (с периодами п и пт} следует, что В = —А.

Определим коэффициенты А и С. Для нахождения А умножим обе части тождества

(23) на B 1 (t — x) и перейдем к пределу t ^ x.

Слева:

B 1 (t — a)B₄(t + a) ^ 1 (x — a)B₄(x + а)

lim---^---:---=----—г-----.

t^x      d₄(t + x)             ^₄(2x)

Справа:

-d ' ₄(t + x) lim[A • ^ 1 (t — x) + В •           $ 1 (t — x) + C^ 1 (t — x)].

t^x                   V₄(t + x)

Так как d 1 (t — x) ~ ^ 1 (0)(t — x) при t ^ x, второе и третье слагаемые стремятся к нулю.

Остается А • ^ 1 (0).

Приравнивая:

Следовательно,

А • ^ 1 (0)

B 1 (x — a~)B₄(x + a) ^₄(2x)

^ 1 (x — a)^₄(x + a)

^А =   ^ 1 (0)^₄(2x)

Коэффициент С находим из (23) при t = a:

С = А-

‘ ^1             ^ 4

— ⁽x — a' ⁾ + — ⁽ x + a' ⁾

Интегрируя (23) от 0 до пи используя свойства логарифмических производных тета-функций, получаем:

I(x) = А • f ^r(t — x}dt + В • f 0 ^ 1                    0

= А1п

В 1 (п — x) ^ 1 (—x)

+ В1п

” $4 ,

— (t + x}dt + Сп =

I ^ 4

d₄(n + x} ^ 4 ⁽ x)

+ Cn.

С учетом (24), (25) и (26) окончательное аналитическое выражение интеграла (22) принимает вид:

01(х - а)<Ш + a) pi              J л%)  ;  ~—[ё;(%-а)+54(%+а)]'           (27)

Численное тестирование квадратурной формулы (11) проводилось для двух значений нормы q при фиксированном а = ;/2.

Таблица 2. Значения интеграла (22) при а = ;/2 , q = 0.1 и и = 11

X	/ (x)	/ 11 (x)	△ 11
0.0785	0.661863650185544	0.661863650268924	-0.8338 х 10-10
0.2356	1.749342391815649	1.749342392037261	-0.2216 х 10-9
0.3927	2.369877414935797	2.369877415240109	-0.3043 х 10-9
0.5498	2.627698464144353	2.627698464489365	-0.3450 х 10-9
0.7069	2.705974181006287	2.705974181373787	-0.3675 х 10-9

Как видно из табл. 2, при п = 11 абсолютная погрешность не превышает 3,7 х 10^-10.

Таблица 3. Значения интеграла (20) при а = ;/2 , q = 0,5 и п = 25

X	/ (X)	hsC*)	△ 2 s
0.0785	78.118507750736	78.118507788496	-0.3776 х 10-7
0.2356	79.645634438020	79.645634475880	-0.3786 х 10-7
0.3927	43.047040002891	43.047040022731	-0.1984 х 10-7
0.5498	22.258565385914	22.258565395741	-0.9827 х 10-8
0.7069	12.178767398389	12.178767403511	-0.5122 х 10-8

Для случая q = 0.5 требуется более мелкая дискретизация; при п = 25 абсолютная погрешность составляет величину порядка 10-7 ^ 10-8, что подтверждает сходимость квадратурной формулы.

Задача 4. Рассмотрим функцию

^(t, %) = esint,  % = 1.

Тогда

re _e slnt

^;(%) = I ^—77-- dt^dt'

0 sm(t - 1)

В заданной области при попытке воспользоваться программными комплексами "Математика" результата в элементарных функциях для интеграла (28) не получаем.

Для интеграла (28) численные значения, полученные по квадратурной формуле (11) при п = 5 и п = 11 , составили 1₅ » 4.695574969512963 и 1 11 « 4.693786801341791 . На основе экстраполяции Ричардсона с отношением шагов г = h 1 /h₂ = 11/5 = 2.2 определен фактический порядок точности метода. Анализ показал, что при предположении о четвертом порядке (р = 4) экстраполированное значение / _экст р ~ 4.69378 практически совпадает с 1 11 , тогда как при р = 2 наблюдается заметное расхождение (/ _экстр « 4.69118), что свидетельствует о фактическом порядке точности р = 4. Оценка погрешности составила £ 11 ~ 7 х 10^-6 для п = 11 и £₅ ~ 1.6 х 10^-4 для п = 5.

Таким образом, использование экстраполяции Ричардсона позволяет уточнить численный результат до ~ 10-6 и подтверждает высокую эффективность примененной квадратурной формулы.

Конструктивные особенности и связь с неклассическими логиками

Предложенная квадратурная формула (11) обладает рядом преимуществ по сравнению с известной формулой (1). В ней отсутствует деление на cos %, что исключает потерю точности в окрестности точек х = п/2 + пт. Структура формулы позволяет сохранять устойчивость при любых значениях аргумента х, включая случаи, когда знаменатель формулы (11) sin(2(t_fc — х)) обращается в нуль; при этом соответствующие предельные переходы могут быть выполнены аналитически. Все проведенные операции в задачах 1–4 выполнимы и согласуются с теорией алгоритмов [13].

Проведенные численные эксперименты с тригонометрическими функциями и тета-функциями Якоби демонстрируют высокую точность метода. Для гладких функций уже при малых п достигается точность, близкая к машинной. Для более сложных функций, таких как тета-функции Якоби, требуется увеличение числа узлов, однако полученные результаты подтверждают сходимость и устойчивость квадратурной формулы.

В интуиционистской логике не допускается закон исключенного третьего в бесконечных конструкциях. В данной работе: все аппроксимации конечны; точность устанавливается конструктивно: для любого £ > 0 можно указать п (зависящее от ^ и х), такое что погрешность меньше £. В работе не используется принцип выбора в неконструктивной форме. Формула может быть верифицирована в системах типа CIC (Calculus of Inductive Constructions) или Agda [14-16], где требуется явное предъявление вычислимой функции.

Оценка погрешности

Для гладких функций (дважды непрерывно дифференцируемых по t) тригонометрическая интерполяция обеспечивает погрешность аппроксимации 0(1/п²). Интегрирование с ядром Гильберта, являющимся сингулярным оператором, сохраняет этот порядок.

Для произвольных вычислимых функций полное аналитическое обоснование выходит за рамки данной работы. В рамках конструктивного подхода для каждого конкретного вычисления предлагается контроль погрешности методом Ричардсона, что позволяет получить результат с гарантированной точностью.

Заключение

В работе получена новая квадратурная формула интерполяционного типа для сингулярных интегралов с ядром Гильберта, удовлетворяющая требованиям: конструктивности (все шаги алгоритмически выполнимы); вычислимости (отсутствие неэффективных предельных переходов); совместимости с неклассическими логиками (финитность, явные оценки). Формула построена на основе аппроксимации плотности тригонометрическим полиномом в равноотстоящих узлах и точного интегрирования полученного выражения. Предложенная квадратурная формула свободна от недостатков, присущих известным подходам, связанных с делением на cos х и численной неустойчивостью в особых точках. Проведенное тестирование на аналитических примерах с тригонометрическими функциями и тета-функциями Якоби подтвердило корректность полученного соотношения и его высокую точность. Представленный метод может быть использован при численном решении сингулярных интегральных уравнений, возникающих в задачах аэродинамики, теории упругости, математической физики, также в алгоритмических системах символьных вычислений, формальных доказательствах в интуиционистской математике, где требуется явное предъявление алгоритма и других областях.