УДК 517.58, 519.65

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

f(x) @ Ґ е n=1  fn Jn(jn x),  n і -1/2,  0 Ј x Ј 1, (1)

где {j_n}-положительные корни уравнения J_n(z) = 0, пронумерованные в порядке возрастания, и

fn = 2 Jn+1(jn)2 у х 1 0  Jn(jn t) t f(t) dt,  n = 1,2,... (2)

   Ряд (1) сходится в весовом пространстве L₂([0,1],x) (т.е. при т₀¹x|f(x)|² dx < Ґ).
   2. Предположим, что x₀ = 0 < x₁ < .. < x_l = 1, q і 1, f(x) О C^2q[x_k-1,x_k], k = 2,...,l, f(x) О C^2q[e,x₁] при любом e, 0 < e < x₁ и введем в рассмотрение оператор

Ln · = ж з и - 1 x   d dx  x  d dx  + n2 x2 ц ч ш ·   , x > 0. (3)

   Имеем (см.[2], глава 7)  L_n ( J_n(j x)) = j² J_n(j x), j > 0. В окрестности точки x О (0,1), x П {x_k}, многократным интегрированием по частям получим (L_n⁰-единичный оператор)

у х Jn(j x) x f(x) dx = j-2q у х Lnq(f(x)) x Jn(j x) dx + + x  ж и Jn(j x) q-1 е r=0  j-2r-2 (Lnrf(x))ў- Jn ў(j x) q-1 е r=0  j-2r-1 (Lnrf(x)) ц ш . (4)

   Предположим теперь, что (n x d / dx (L_n^rf(x)) - L_n^rf(x)) ® 0, x ® 0 при r = 0,1,...,(q - 1) и L_n^qf(x)) О L₂(0,x₁) . Отсюда и из (4), учитывая поведение J_n(x) при x ® 0 (см.[1,2]), придем к асимптотическому разложению вида f_n = Y_n^q + R_n ^q по степеням j_n^-1 с точностью до R_n^q = o(n^-2q), n ® Ґ, где

Ynq = 2 Jn+1(jn)2 l е k=1  xk ж и Jn(jnxk) q е r=1  jn-2r A1rk - Jnў(jnxk) q е r=1  jn-2r+1 A0rk ц ш ,                          Rn q = 2 jn2q Jn+1(jn)2 l-1 е k=0  у х xk+1 xk  Lnq(f(t)) t Jn(jnt) dt,  n = 1,2,...                               (5)

и обозначено A⁰_rl = L_n^rf(1), A¹_rl = (L_n^rf(1))ў,  A⁰_rk = L_n^rf(x_k - 0) - L_n^rf(x_k + 0), A¹_rk = (L_n^rf(x_k - 0))ў-L_n^rf(x_k + 0))ў,  k = 1,..., l - 1,  r = 1,...,q;  ( (·)^ў = d/dx(·) ).
   3. Приведем свод известных формул (см. главу 18 в [1] и формулы (55) и (58)-(61)из глав 7.2 и 7.15 в [2]), лежащих в основе дальнейшего исследования:

±nJn(x) + x Jn ў (x) = ±x Jn-±1(x); (6)

Fn(x,X,z) = Ґ е n = 1  Jn(x jn) Jn(X jn) ( jn2 - z2 )   J1 + n(jn)2 = = p Jn(x z)  ( ( Jn(X z)  Yn(z) ) - Jn(z) Yn(X z) ) 4  Jn(z) ,  z П {±jn}, (7)

где 0 Ј x Ј X Ј 1. Отсюда при z ® 0 (см.[2], глава 7)

Fn(x,X,0) = Ґ е n = 1  Jn(x jn) Jn(X jn) jn2   J1 + n(jn)2 = p xn ( -1 + X2 n )  csc(p n) 4 Xn  G(1 - n)  G(1 + n) . (8)

Заметим, что если перейти к пределу n ® 0 и выражение справа примет вид logX /2.
   Функция F_n(x,X,z) будет определена при любом z О C, если при целом положительном m и z = ±j_m принять

Fn(x,X,jm) =     Jn(x jn) Jn(X jn) ( jn2 - jm2 )   J1 + n(jn)2 =     ж з и Fn(x,X,z) - Jn(x jm) Jn(X jm) ( j m2 - z2 )   J1 + n(jm)2 ц ч ш ,   0 Ј x Ј X Ј 1. (9)

   Предел в (9) справа существует и без труда может быть представлен в явном виде (например, применением системы Mathematica, см. [10]).
   Как известно (см., например, [1], главы 12 и 15), j_n = p(n + n/2 - 1/4) + O(n^-1)  , n ® Ґ,  и при x > 0 функция J_m(x) ограничена при любом m і -1/2, а также |J_n+1(j_n) | і const > 0. Отсюда следует равномерная сходимость рядов в формулах (6)-(9) по указанным значениям x, X и z О C.
   Рассмотрим теперь функцию (x О [0,1])

Gn(x,x,z) = м н о Fn(x,x,z) при    0 Ј x < x; Fn(x,x,z) при     x Ј x Ј 1.                                               (10)

   Учитывая асимптотику функций Бесселя на бесконечности, функцию G_n можно дифференцировать по x  (при x № x), а также многократно дифференцировать по w = z².
   4. Предположим теперь, что нам известна функция f_q(x) =. Тогда разложение f_N(x) = f_q(x) + сходится к f(x) в метрикеL₂[(0,1),x] (см.[1]) со скоростью ||f - f_N|| = o(N^-2q+1), N ® Ґ. Тот же порядок равномерной сходимости будет обеспечен и на отрезках вида [e,1], 0 < e < 1.
   Если нам известны только коэффициенты и точки {x_k}, то практическая реализация этой схемы ускоренной сходимости состоит из следующих алгоритмов.
   A. Нахождение скaчков , s = 0,1, k = 1,..., l - 1, r = 1,...,q. Для достаточно большого N и некоторого 0 < d = const < 1 подберем целые значения {n_s} так, чтобы было dN Ј n₁ < n₂ < ј < n_p-1 < n_p = N, p = 2(l + q), и найдем приближенные значения ,-с точностью -= o(N^2(r-q)-1), -= o(N^2(r-q)), N ® Ґ,-как решения линейной системы уравнeний = 0, n = n_r, r = 1, 2 ,..., 2(l + q). Не останавливаясь на вопросах разрешимости этой системы, отметим только, что, как показывают асимптотические формулы для функций Бесселя, она аналогична соответствующей системе в КЭГ-методе.
   B. Нахождение функции f_q(x). Из формул и заключительных замечаний п.3 следует, что суммы вида

Ґ е n=1   Jn(jn x) Jn(jn x) Jn+1(jn)2  jn 2r  ,     Ґ е n=1   jn  Jn(jn x)  Jn ў(jn x) Jn+1(jn)2  jn 2r  ,  r і 1 (11)

явно выражаются через функцию G_n(x,x,z) и ее соответствующие производные при z = 0. Таким образом, функция f_q(x) находится приближенно, в явном виде, с использованием результатов алгоритма А.
   Описанный алгоритм ускорения сходимости (назовем его алгоритмом AB или АВ-методом) является естественным распространением КЭГ-метода наразложения Фурье-Бесселя, причем в роли полиномов Бернулли выступают упомянутые значения функции G_n(x,x,z) и ее производных при z = 0.
   5. Переходя к распространению результатов работы [8] на ряды Фурье-Бесселя, заметим, что в представлении (5) для коэффициентов {Y_n^q},-при каждом фиксированном индексе суммирования k,-фигурируют два асимптотических ряда, соответственно по четным и нечетным степеням j_n^-1. При q і 2 к каждому из них применим аппроксимацию Паде для рядов по степеням w = j_n^-2 некоторого, зависящего от k, порядка [M/m], M + m = (q - 1), m і 1 (см.[9]).
   В результате этой основной процедуры, вместо упомянутых коэффициентов {Y_n^q}, будут фигурировать следующие:

Fnq = l е k=1  xk ж з и jn 2(mk1-Mk1) Jn(jnxk) Pk1(jn2) Jn+1(jn) Qk1(jn2) + jn 2(mk2-Mk2) jn Jn ў(jnxk) Pk2(jn2) Jn+1(jn)2 Qk2(jn2) ц ч ш , n і 1, (12)

где P_k^s и Q_k^s-многочлены от j_n² порядков (M_ks-1) и m_ks соответственно, s = 1,2.
   Здесь мы приходим к алгоритму, основанному на явном виде функции j_q(x) = .

   C. Нахождение функции j_q(x). Разложив рациональные функции P_k^s(x)/Q_k^s(x), s = 1,2, на простые дроби, нетрудно убедиться, что функция j_q(x) является линейной комбинацией значений G_n(x,x,z), ее первых производных по x, x № x, и (r - 1) производных по w = z², при x = x_k, но уже при z² = w_kts, k = 1, ..., l; t = 1, ...m_ks, s = 1,2, где величины {w_kts} зависят от нулей знаменателей примененных аппроксимаций Паде, а r-максимум их кратности. Таким образом, функция j_q(x) также находится приближенно (с учетом результатов алгоритма А), в явном виде.
   На этот раз,-уже по известной функции j_q(x),-найдем разложение g_N(x) = j_q(x) + которое,-при условии нахождения скачков с достаточной точностью,-сходится (см. [9]) к f(x) (как и f_N) со скоростью ||f - g_N|| = o(N^-2q+1), N ® Ґ.
   Алгоритм A+C (назовем его алгоритмом AC или АС-методом) является естественным аналогом алгоритма квазиполиномиального разложения работы [8].
   6. Численную реализацию алгоритмов АВ и АС проиллюстрируем в простейшем (и наиболее важном для приложений) случае, когда f(x)-гладкая при x О (0,1] функция, удовлетворяющая в окрестности нуля указанным выше условиям. В этом случае, пользуясь известными формулами для функцийБесселя, из (5) (учитывая, что J_n(j_n) = 0 и l = 1), получим

Ynq = 2 jn Jn+1(jn)   q е r=0  Ar jn-2 r,  Ar = Ar10,  A0 = 0  n і 1. (13)

   Продифференцировав формулу (7) по X, получим при X®1 (Gў_n = dG_n /dX,см. также (6))

Gўn(1,x,z) = Ґ е n = 1  jn Jn(x jn) ( jn2 - z2 )   J1 + n(jn) = Jn(x z) 2Jn(z) . (14)

   Дифференцируя теперь Gў_n(1,x,z) по z² и устремляя z®0, найдем суммы рядов r = 0,1, ..., в виде явных функций, являющихся в АВ-методе аналогами полиномов Бернулли в КЭГ-методе. Вот как выглядят первые три из них (r = 0,1,2):

м н о xn,     xn  (1 - x2 ) 4  ( 1 + n) ,     xn ( x2 -1)   ( -3 - n+ x2 ( 1 + n) ) 32  ( 1 + n) 2  ( 2 + n) ь э ю .

   Заметим теперь, что, как нетрудно видеть, что для g(x) О C^q[0,1], q і 1 и n і -1/2 функция f(x) = xⁿg(x²) удовлетворяет всем условиям, приведенным в п.2 выше.
   Алгоритмы АВ и АС применялись к такой функции

f(x) = xn ( cos(25 x) + log(1 + x2) ). (15)

   Эксперимент проводился на персональном компьютере применениемсистемы Mathematica ([10]). Предполагалось, что q = 2 p-четное и в алгоритме С к разложению  w = j_n^-2, применялась диагональная аппроксимация Паде порядка [p/p].