Журнал вычислительной математики и математической физики, 2021, T. 61, № 6, стр. 926-935

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

О. В. Лешонков^1, *, Е. А. Соболева^2, **, А. Д. Чернышов^2, ***

¹ АО НИИЭТ
394033 Воронеж, ул. Старых большевиков, 5, Россия

² ВГУИТ
394036 Воронеж, пр-т Революции, 19, Россия

^* E-mail: chernyshovad@mail.ru
^** E-mail: sobol5661@yandex.ru
^*** E-mail: leekripper@yandex.ru

Поступила в редакцию 15.12.2018
После доработки 16.11.2020
Принята к публикации 16.12.2020

DOI: 10.31857/S0044466921060065

Полный текст (PDF)

Аннотация

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

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

1. ВВЕДЕНИЕ

Известны различные методы вычисления сложных интегралов в классическом смысле с переменным верхним пределом, достаточно подробно изложенные в [1]. В [2] приведены квадратурные формулы наивысшей тригонометрической степени точности, когда используются значения функции в смещенных точках по отношению к равномерно нанесенным узлам данного отрезка. В [3], [4] рассмотрены интегралы от функций с некоторыми особенностями. Для вычисления интегралов часто применяют различные разложения. Оценка погрешности разложений в L₂ приведена в [5]. В [6] используются быстрые синус-разложения, когда функцию представляют суммой граничной функции четного порядка M_2p специальной конструкции и рядом Фурье, исследуется сходимость ряда, дается оценка погрешности. Скорость сходимости такого ряда значительно повышается при увеличении порядка 2p граничной функции.

При рассмотрении более простой проблемы – вычисления интегралов с постоянными пределами, для прикладных целей чаще используется метод Симпсона, или пакет Maple. Известен также метод Ромберга, как техника, позволяющая повышать точность вычислений интегралов путем составления определенной линейной комбинации их значений, полученных на разных сетках интегрирования [7]. В методе Клиншоу–Куртиса предлагается вводить замену переменной и затем использовать классические ряды Фурье [8].

Ниже предлагается новый метод быстрых синус-разложений (МБСР), основанный на быстрых разложениях [6], который эффективен по своей простоте и особенно высокой точности и удобен для практического применения. Это будет видно из численных экспериментов. Метод применим для решения многомерных краевых задач, позволяет решать различные сложные краевые и нелинейные интегродифференциальные задачи [9], [10]. В методе используется поточечный способ вычисления коэффициентов Фурье.

В связи с использованием в данной работе метода коллокаций следует отметить работу [11], в которой решается однa из сложных задач для уравнений Навье–Стокса. Здесь используется комплексный метод, состоящий из следующих действий: метод коллокаций+наименьших квадратов+итерационный+вычисление поправок по методу Крылова для организации итераций. При этом используются разложения по мономам с некоторыми степенями. Каждый из перечисленных шагов улучшает конечный результат, который получается с высокой точностью. Но сравнить его с предлагаемым методом быстрых разложений не представляется возможным, так как неясно, каким образом метод из [11] можно применить для нахождения неявно заданной функции при любых $x \in \left[ {0,a} \right]$, а также и для вычисления определенного интеграла с переменным верхним пределом в любой точке заданного отрезка.

1. ВЫЧИСЛЕНИЕ ОПРЕДЕЛЕННЫХ ИНТЕГРАЛОВ С ПЕРЕМЕННЫМ ВЕРХНИМ ПРЕДЕЛОМ ОТ ФУНКЦИЙ, ДЛЯ КОТОРЫХ ПЕРВООБРАЗНАЯ НЕ ВЫРАЖАЕТСЯ ЧЕРЕЗ ЭЛЕМЕНТАРНЫЕ ИЛИ СПЕЦФУНКЦИИ

Функцию $f\left( x \right)$, $x \in \left[ {0,a} \right]$, будем называть сложной, если на данном отрезке первообразная в точном виде не выражается конечным образом через элементарные или спецфункции.

Пусть задана сложная функция одной переменной

(1.1)

$y = f\left( x \right) \in L_{2}^{{(2p + 1)}},\quad 0 \leqslant x \leqslant а,$

где $p$ – некоторое заданное целое число. Здесь рассматриваются только непрерывные гладкие функции вместе со своими производными до порядка $\left( {2p + 1} \right)$ включительно. Разрывы не допускаются в связи с тем, что в работе используются ряды Фурье, которые в случае разрывных функций медленно сходятся и не допускают почленное дифференцирование, что затрудняет их применение в исследованиях прикладного характера. Запишем интеграл:

(1.2)

$\int\limits_0^x {f\left( t \right)dt} ,\quad x \in \left[ {0,a} \right].$

Полагаем, что функция $f\left( x \right)$ сложная и интеграл (1.2) невозможно представить точно в явном конечном виде через элементарные или спецфункции. Подобные интегралы являются особенно проблемными. Для решения таких задач на отрезке $x \in \left[ {0,a} \right]$ представим $f\left( x \right)$ быстрым разложением [6] с оператором $C{{h}_{{2p}}}$:

(1.3)

$\begin{gathered} f\left( x \right) = {{M}_{2}}\left( x \right) + \sum\limits_{m = 1}^N {{{f}_{{2,m}}}\sin m\pi \frac{x}{a}} ,\quad {{M}_{2}}\left( x \right) = f\left( 0 \right)\left( {1 - \frac{x}{a}} \right) + f\left( a \right)\frac{x}{a} + \\ + \;f{\kern 1pt} {\text{''}}\left( 0 \right)\left( {\frac{{{{x}^{2}}}}{2} - \frac{{{{x}^{3}}}}{{6a}} - \frac{{ax}}{3}} \right) + f{\kern 1pt} {\text{''}}\left( a \right)\left( {\frac{{{{x}^{3}}}}{{6a}} - \frac{{ax}}{6}} \right), \\ \end{gathered} $

где для определенности и простоты изложения вначале положим $p = 1$, ${{f}_{{2,m}}}$ – коэффициенты Фурье для разности $f\left( x \right) - {{M}_{2}}\left( x \right)$, $N$ – количество учитываемых членов в ряде Фурье, ${{M}_{2}}\left( x \right)$ – граничная функция, конструкция которой обеспечивает быструю сходимость ряда Фурье и потому $N$ можно брать небольшим для обеспечения высокой точности вычисления интеграла.

Функция ${{M}_{{2p}}}\left( x \right)$ с четным индексом внизу называется граничной, так как ее коэффициенты определяются через значения $f\left( x \right)$ и ее производных четного порядка от 0 до $2p$ включительно на концах отрезка.

Скорость сходимости ряда Фурье в (1.3) существенно возрастает при увеличении порядка $2p$ используемой граничной функции ${{M}_{{2p}}}\left( x \right)$. Быстрое разложение (1.3) обладает замечательными свойствами: оно быстро сходится и его можно почленно трижды дифференцировать, ряд Фурье при этом остается сходящимся на всем отрезке $\left[ {0,a} \right]$. Сравнение метода быстрых синус-разложений с классическим синус-разложением Фурье на примере с функцией ${{\left( {1 + x} \right)}^{{10}}}$ на отрезке $\left[ {0,1} \right]$ показывает, что при использовании МБСР для $N = 3$ погрешность будет порядка ${{10}^{{ - 3}}}$, тогда как в случае рядов Фурье такая точность достигается только при $N > 2000$. Так как $f\left( x \right)$ – сложная функция, то непосредственно применить МБСР не удается из-за невозможности вычисления коэффициентов ${{f}_{{2,m}}}$ – их интегральными выражениями Фурье. Поэтому ниже предлагается поточечный способ вычисления ${{f}_{{2,m}}}$.

Для этого на отрезок $\left[ {0,a} \right]$ нанесем равномерную сетку с $N$ внутренними точками и быстрое разложение функции $f\left( x \right)$ из (1.3) заменим на другое

(1.4)

$f\left( x \right) = {{M}_{2}}\left( x \right) + \sum\limits_{m = 1}^N {f_{{2,m}}^{*}\sin m\pi \frac{x}{a}} ,$

где выражение ${{M}_{2}}\left( x \right)$ остается прежним, но коэффициенты $f_{{2,m}}^{*}$ вычисляются иным способом – поточечным из алгебраической системы, когда равенство (1.4) записывается в узлах равномерной сетки $\left\{ {{{x}_{k}}} \right\}$:

(1.5)

$\sum\limits_{m = 1}^N {f_{{2,m}}^{*}\sin m\pi \frac{{{{x}_{k}}}}{a}} = {{\varphi }_{2}}\left( {{{x}_{k}}} \right),\quad {{\varphi }_{2}}\left( {{{x}_{k}}} \right) = f\left( {{{x}_{k}}} \right) - {{M}_{2}}\left( {{{x}_{k}}} \right),\quad {{x}_{k}} = \frac{{ka}}{{N + 1}},\quad k = 1{\kern 1pt} \div {\kern 1pt} N.$

Ее решение удается представить в явном конечном виде. Для этого левую и правую части равенства (1.5) умножим на $\sin n\pi {{{{x}_{k}}} \mathord{\left/ {\vphantom {{{{x}_{k}}} a}} \right. \kern-0em} a}$ и просуммируем по $k$:

(1.6)

$\sum\limits_{k = 1}^N {\left( {f\left( {{{x}_{k}}} \right) - {{M}_{2}}\left( {{{x}_{k}}} \right)} \right)\sin n\pi \frac{{{{x}_{k}}}}{a}} = \sum\limits_{k = 1}^N {\sum\limits_{m = 1}^N {f_{{2,m}}^{*}\sin m\pi \frac{{{{x}_{k}}}}{a}\sin n\pi \frac{{{{x}_{k}}}}{a}} } .$

В [12] доказано, что множество $\left\{ {\sin n\pi {{{{x}_{k}}} \mathord{\left/ {\vphantom {{{{x}_{k}}} a}} \right. \kern-0em} a}} \right\}$, $\left( {n,k} \right) = 1{\kern 1pt} \div {\kern 1pt} N$ сохраняет свойство ортогональности на бесконечном отрезке, имеющее вид

(1.7)

$\begin{gathered} \sum\limits_{k = 1}^N {\sin m\pi \frac{{{{x}_{k}}}}{a}\sin n\pi \frac{{{{x}_{k}}}}{a}} = 0\quad {\text{при}}\quad m \ne n,\quad \left( {m,n} \right) = 1{\kern 1pt} \div {\kern 1pt} N, \\ \sum\limits_{k = 1}^N {\sin m\pi \frac{{{{x}_{k}}}}{a}\sin n\pi \frac{{{{x}_{k}}}}{a}} = \sum\limits_{k = 1}^N {{{{\sin }}^{2}}n\pi \frac{{{{x}_{k}}}}{a}} \quad {\text{при}}\quad m = n. \\ \end{gathered} $

Можно показать, что свойство (1.7) остается верным и для конечного отрезка $x \in \left[ {0,a} \right]$. Кроме (1.7), отметим еще одно важное предельное свойство коэффициентов $f_{{2,m}}^{*}$ и ${{f}_{{2,m}}}$ при использовании быстрых разложений и равномерном разбиении отрезка $\left[ {0,a} \right]$ на мелкие части с $N$ внутренними точками:

(1.8)

$\mathop {\lim }\limits_{N \to \infty } {{f}_{{2,m}}} = \mathop {\lim }\limits_{N \to \infty } f_{{2,m}}^{*}\quad \forall m = 1{\kern 1pt} \div {\kern 1pt} N.$

Теперь из (1.6) при помощи свойства (1.7) найдем коэффициенты

(1.9)

$f_{{2,m}}^{*} = {{\sum\limits_{k = 1}^N {\left( {f\left( {{{x}_{k}}} \right) - {{M}_{2}}\left( {{{x}_{k}}} \right)} \right)\sin m\pi \frac{{{{x}_{k}}}}{a}} } \mathord{\left/ {\vphantom {{\sum\limits_{k = 1}^N {\left( {f\left( {{{x}_{k}}} \right) - {{M}_{2}}\left( {{{x}_{k}}} \right)} \right)\sin m\pi \frac{{{{x}_{k}}}}{a}} } {\sum\limits_{k = 1}^N {{{{\sin }}^{2}}m\pi \frac{{{{x}_{k}}}}{a}} }}} \right. \kern-0em} {\sum\limits_{k = 1}^N {{{{\sin }}^{2}}m\pi \frac{{{{x}_{k}}}}{a}} }},\quad m = 1{\kern 1pt} \div {\kern 1pt} N.$

Свойство (1.8) позволяет предположить: если для рассматриваемых гладких функций разложение (1.3) с рядом Фурье в его правой части быстро сходится, то и разложение (1.4), где $f_{{2,m}}^{*}$ определяются из поточечной системы (1.9), тоже быстро сходится. Это предположение подтверждается многочисленными численными экспериментами, которые приводятся ниже. Доказательство такого свойства приводится в [13].

Представляя $f\left( x \right)$ поточечным равенством (1.4), где $f_{{2,m}}^{*}$ находятся по формуле (1.9), с высокой точностью получим искомое выражение для определенного интеграла (1.2) с переменным верхним пределом:

(1.10)

Правая часть (1.10) допускает возможность дифференцирования до второго порядка включительно. С ростом порядка производной относительная погрешность растет (см. табл. 1).

Таблица 1

$\Delta f\left( x \right)$	$\Delta f{\kern 1pt} {\text{'}}\left( x \right)$	$\Delta f{\kern 1pt} {\text{''}}\left( x \right)$
$\begin{gathered} \Delta = 2.38 \times {{10}^{{ - 3}}} \hfill \\ x = 1 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 7 \times {{10}^{{ - 5}}} \hfill \\ x = 0.96 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 4.98 \times {{10}^{{ - 4}}} \hfill \\ x = 1 \hfill \\ \end{gathered} $

В табл. 1 приведены максимальные относительные погрешности при использовании быстрого разложения с оператором $C{{h}_{2}}$ второго порядка при $N = 10$ для функции $f\left( x \right) = {{\left( {1 + x} \right)}^{{10}}}$.

Точность вычисления интеграла (1.10) значительно повышается при использовании в быстром разложении (1.3) вместо ${{M}_{2}}\left( x \right)$ граничной функции ${{M}_{4}}\left( x \right)$ более высокого (четвертого) порядка

(1.11)

$\begin{gathered} f\left( x \right) = {{M}_{4}}\left( x \right) + \sum\limits_{m = 1}^N {f_{{4,m}}^{*}\sin m\pi \frac{x}{a}} ,\quad {{M}_{4}}\left( x \right) = f\left( 0 \right)\left( {1 - \frac{x}{a}} \right) + f\left( a \right)\frac{x}{a} + \\ + \;f{\kern 1pt} {\text{''}}\left( 0 \right)\left( {\frac{{{{x}^{2}}}}{2} - \frac{{{{x}^{3}}}}{{6a}} - \frac{{ax}}{3}} \right) + f{\kern 1pt} {\text{''}}\left( a \right)\left( {\frac{{{{x}^{3}}}}{{6a}} - \frac{{ax}}{6}} \right) + \\ + \;{{f}^{{\left( 4 \right)}}}\left( 0 \right)\left( {\frac{{{{x}^{4}}}}{{24}} - \frac{{{{x}^{5}}}}{{120a}} - \frac{{a{{x}^{3}}}}{{18}} + \frac{{{{a}^{3}}x}}{{45}}} \right) + {{f}^{{\left( 4 \right)}}}\left( a \right)\left( {\frac{{{{x}^{5}}}}{{120a}} - \frac{{a{{x}^{3}}}}{{36}} + \frac{{7{{a}^{3}}x}}{{360}}} \right). \\ \end{gathered} $

Введем обозначения

(1.12)

${{\varphi }_{4}}\left( {{{x}_{k}}} \right) = f\left( {{{x}_{k}}} \right) - {{M}_{4}}\left( {{{x}_{k}}} \right),\quad {{\varphi }_{4}}\left( {{{x}_{k}}} \right) = \sum\limits_{m = 1}^N {f_{{4,m}}^{*}\sin m\pi \frac{{{{x}_{k}}}}{a}} ,\quad {{x}_{k}} = \frac{{ka}}{{N + 1}},\quad k = 1{\kern 1pt} \div {\kern 1pt} N.$

Используя обозначение ${{\varphi }_{4}}\left( {{{x}_{k}}} \right)$, данное в первом равенстве (1.12), коэффициенты $f_{{4,m}}^{*}$, подобно (1.9), можно записать в виде

(1.13)

$f_{{4,m}}^{*} = {{\sum\limits_{k = 1}^N {{{\varphi }_{4}}\left( {{{x}_{k}}} \right)\sin m\pi \frac{{{{x}_{k}}}}{a}} } \mathord{\left/ {\vphantom {{\sum\limits_{k = 1}^N {{{\varphi }_{4}}\left( {{{x}_{k}}} \right)\sin m\pi \frac{{{{x}_{k}}}}{a}} } {\sum\limits_{k = 1}^N {{{{\sin }}^{2}}m\pi \frac{{{{x}_{k}}}}{a}} }}} \right. \kern-0em} {\sum\limits_{k = 1}^N {{{{\sin }}^{2}}m\pi \frac{{{{x}_{k}}}}{a}} }},\quad m = 1{\kern 1pt} \div {\kern 1pt} N.$

Подставляя $f\left( x \right)$ и ${{M}_{4}}\left( x \right)$ из (1.11) в (1.2), после интегрирования будем иметь формулу для вычисления определенного интеграла с переменным верхним пределом, когда используется БСР с оператором $C{{h}_{4}}$:

(1.14)

Квадратура (1.14) допускает вычисление производных до четвертого порядка включительно, но при этом растет относительная погрешность (см. табл. 2), составленной для $f\left( x \right) = {{\left( {1 + x} \right)}^{{10}}}$ при $N = 10$.

Таблица 2

$\Delta f\left( x \right)$	$\Delta f{\kern 1pt} {\text{'}}\left( x \right)$	$f{\kern 1pt} {\text{''}}\left( x \right)$	$\Delta f{\kern 1pt} {\text{'''}}\left( x \right)$	$\Delta {{f}^{{\left( 4 \right)}}}\left( x \right)$
$\begin{gathered} \Delta = 1.1 \times {{10}^{{ - 5}}} \hfill \\ x = 1 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 3.38 \times {{10}^{{ - 7}}} \hfill \\ x = 0.96 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 2.1 \times {{10}^{{ - 6}}} \hfill \\ x = 1 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 1.85 \times {{10}^{{ - 5}}} \hfill \\ x = 0.96 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 1.8 \times {{10}^{{ - 4}}} \hfill \\ x = 1 \hfill \\ \end{gathered} $

Скорость сходимости ряда для квадратуры еще более повышается при использовании в (1.11) вместо ${{M}_{4}}\left( x \right)$ граничной функции ${{M}_{6}}\left( x \right)$, т.е.

(1.15)

$\begin{gathered} f\left( x \right) = {{M}_{6}}\left( x \right) + \sum\limits_{m = 1}^N {{{f}_{{6,m}}}\sin m\pi \frac{x}{a}} , \\ {{M}_{6}}\left( x \right) = f\left( 0 \right)\left( {1 - \frac{x}{a}} \right) + f\left( a \right)\frac{x}{a} + f{\kern 1pt} {\text{''}}\left( 0 \right)\left( {\frac{{{{x}^{2}}}}{2} - \frac{{{{x}^{3}}}}{{6a}} - \frac{{ax}}{3}} \right) + \\ + \;f{\kern 1pt} {\text{''}}\left( a \right)\left( {\frac{{{{x}^{3}}}}{{6a}} - \frac{{ax}}{6}} \right) + {{f}^{{\left( 4 \right)}}}\left( 0 \right)\left( {\frac{{{{x}^{4}}}}{{24}} - \frac{{{{x}^{5}}}}{{120a}} - \frac{{a{{x}^{3}}}}{{18}} + \frac{{{{a}^{3}}x}}{{45}}} \right) + \\ + \;{{f}^{{\left( 4 \right)}}}\left( a \right)\left( {\frac{{{{x}^{5}}}}{{120a}} - \frac{{a{{x}^{3}}}}{{36}} + \frac{{7{{a}^{3}}x}}{{360}}} \right) + {{f}^{{\left( 6 \right)}}}\left( 0 \right)\left( {\frac{{{{x}^{6}}}}{{6{\kern 1pt} !}} - \frac{{{{x}^{7}}}}{{7{\kern 1pt} !a}} - \frac{{a{{x}^{5}}}}{{4 \times 90}} + \frac{{{{a}^{3}}{{x}^{3}}}}{{3 \times 90}} - \frac{{4{{a}^{5}}x}}{{6! \times 3}}} \right) + \\ + \;{{f}^{{\left( 6 \right)}}}\left( a \right)\left( {\frac{{{{x}^{7}}}}{{7{\kern 1pt} !a}} - \frac{{a{{x}^{5}}}}{{6{\kern 1pt} !}} + \frac{{7{{a}^{3}}{{x}^{3}}}}{{36 \times 60}} - \frac{{31{{a}^{5}}x}}{{7{\kern 1pt} !\, \times 3}}} \right). \\ \end{gathered} $

Вспомогательную функцию ${{\varphi }_{6}}\left( {{{x}_{k}}} \right)$, подобно (1.12), определим равенством

(1.16)

${{\varphi }_{6}}\left( {{{x}_{k}}} \right) = f\left( {{{x}_{k}}} \right) - {{M}_{6}}\left( {{{x}_{k}}} \right),\quad {{\varphi }_{6}}\left( {{{x}_{k}}} \right) = \sum\limits_{m = 1}^N {f_{{6,m}}^{*}\sin m\pi \frac{{{{x}_{k}}}}{a}} ,\quad {{x}_{k}} = \frac{{ka}}{{N + 1}},\quad k = 1{\kern 1pt} \div {\kern 1pt} N.$

Коэффициенты $f_{{6,m}}^{*}$, по аналогии с (1.13), найдем из выражения

(1.17)

$f_{{6,m}}^{*} = {{\sum\limits_{k = 1}^N {{{\varphi }_{6}}\left( {{{x}_{k}}} \right)\sin m\pi \frac{{{{x}_{k}}}}{a}} } \mathord{\left/ {\vphantom {{\sum\limits_{k = 1}^N {{{\varphi }_{6}}\left( {{{x}_{k}}} \right)\sin m\pi \frac{{{{x}_{k}}}}{a}} } {\sum\limits_{k = 1}^N {{{{\sin }}^{2}}m\pi \frac{{{{x}_{k}}}}{a}} }}} \right. \kern-0em} {\sum\limits_{k = 1}^N {{{{\sin }}^{2}}m\pi \frac{{{{x}_{k}}}}{a}} }},\quad m = 1{\kern 1pt} \div {\kern 1pt} N.$

При использовании разложения (1.15) интеграл (1.2) с переменным верхним пределом представляется формулой

(1.18)

$\begin{gathered} \int\limits_0^x {f\left( t \right)dt} = f\left( 0 \right)\left( {x - \frac{{{{x}^{2}}}}{{2a}}} \right) + f\left( a \right)\frac{{{{x}^{2}}}}{{2a}} + f{\kern 1pt} {\text{''}}\left( 0 \right)\left( {\frac{{{{x}^{3}}}}{6} - \frac{{{{x}^{4}}}}{{24a}} - \frac{{a{{x}^{2}}}}{6}} \right) + f{\kern 1pt} {\text{''}}\left( a \right)\left( {\frac{{{{x}^{4}}}}{{24a}} - \frac{{a{{x}^{2}}}}{{12}}} \right) + \\ + \;{{f}^{{\left( 4 \right)}}}\left( 0 \right)\left( {\frac{{{{x}^{5}}}}{{120}} - \frac{{{{x}^{6}}}}{{720a}} - \frac{{a{{x}^{4}}}}{{72}} + \frac{{{{a}^{3}}{{x}^{2}}}}{{90}}} \right) + {{f}^{{\left( 6 \right)}}}\left( 0 \right)\left( {\frac{{{{x}^{7}}}}{{7 \times 6{\kern 1pt} !}} - \frac{{{{x}^{8}}}}{{8 \times 7{\kern 1pt} !a}} - \frac{{a{{x}^{6}}}}{{24 \times 90}} + \frac{{{{a}^{3}}{{x}^{4}}}}{{12 \times 90}} - \frac{{4{{a}^{5}}{{x}^{2}}}}{{6 \times 6{\kern 1pt} !}}} \right) + \\ + \;{{f}^{{\left( 6 \right)}}}\left( a \right)\left( {\frac{{{{x}^{8}}}}{{8 \times 7{\kern 1pt} !a}} - \frac{{a{{x}^{6}}}}{{6 \times 6{\kern 1pt} !}} + \frac{{7{{a}^{3}}{{x}^{4}}}}{{36 \times 240}} - \frac{{31{{a}^{5}}{{x}^{2}}}}{{6 \times 7{\kern 1pt} !}}} \right) + \\ + \;{{f}^{{\left( 4 \right)}}}\left( a \right)\left( {\frac{{{{x}^{6}}}}{{720a}} - \frac{{a{{x}^{4}}}}{{144}} + \frac{{7{{a}^{3}}{{x}^{2}}}}{{720}}} \right) + \sum\limits_{m = 1}^N {\frac{a}{{m\pi }}f_{{6,m}}^{*}\left( {1 - \cos m\pi \frac{x}{a}} \right)} . \\ \end{gathered} $

Выражение квадратуры (1.18) можно дифференцировать 6 раз, при этом ряды Фурье будут оставаться сходящимися. С ростом порядка производной относительная погрешность возрастает (см. табл. 3, составленную для $f\left( x \right) = {{\left( {1 + x} \right)}^{{10}}}$ при $N = 10$).

Таблица 3

$\Delta f{\text{''}}\left( x \right)$	$\Delta f{\text{'''}}\left( x \right)$	$\Delta {{f}^{{\left( 4 \right)}}}\left( x \right)$	$\Delta {{f}^{{\left( 5 \right)}}}\left( x \right)$	$\Delta {{f}^{{\left( 6 \right)}}}\left( x \right)$
$\begin{gathered} \Delta = 4.15 \times {{10}^{{ - 9}}} \hfill \\ x = 1 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 3.71 \times {{10}^{{ - 8}}} \hfill \\ x = 0.96 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 3.12 \times {{10}^{{ - 7}}} \hfill \\ x = 1 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 3.55 \times {{10}^{{ - 6}}} \hfill \\ x = 0.96 \hfill \\ \end{gathered} $	$\begin{gathered} \Delta = 4.93 \times {{10}^{{ - 5}}} \hfill \\ x = 1 \hfill \\ \end{gathered} $

Из табл. 3 видно, что производная шестого порядка при $N = 10$ вычисляется слишком грубо, поэтому для ее вычисления необходимо брать $N > 10$.

При одинаковом значении $N$ квадратура (1.14) более точная по сравнению с (1.10), а квадратура (1.18) еще точнее по сравнению с (1.14). То есть с ростом порядка $2p$ оператора $C{{h}_{{2p}}}$, как и с увеличением $N$ – числа учитываемых членов в рядах Фурье, погрешность быстро уменьшается. По этой причине предложенный здесь МБСР целесообразно применять для проведения высокоточных расчетов, например, космического характера.

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

Ниже в таблицах приведено сравнение данного в работе метода быстрых синус-разложений с методом Симпсона, когда для рассматриваемой подынтегральной функции $f\left( x \right)$ существует первообразная и оба предела постоянные $\left[ {0,\;1} \right]$. Это дает возможность сравнить данные два метода с точным значением интеграла. В первом примере выбрана гладкая функция $\sin 0,3\pi x$, во втором примере быстро осциллирующая $\sin 5.3\pi x$, в третьем функция, имеющая участок резкого возрастания ${{\left( {1 + x} \right)}^{{10}}}$. Сравнения сделаны для случаев, когда в быстрых разложениях используются три граничные функции четных порядков ${{M}_{2}}$, ${{M}_{4}}$, ${{M}_{6}}$. При этом количество внутренних точек бралось одинаковым и равным $N = 10$. Для метода Симпсона величина $N$ подбиралась так, чтобы абсолютная погрешность $\Delta $ была одинаковой с методом БСР. Результаты приведены в табл. 4.

Таблица 4

Функция	Метод		Метод		Метод
Функция	Симпсон	МБСР Сh₂	Симпсон	МБСР Сh₄	Симпсон	МБСР Сh₆
sin(0.3πx)	Δ = 4.66 × 10^–9		Δ = 2.44 × 10^–12		Δ = 1.42 × 10^–15
N	28	10	170	10	194	10
sin(5.3πx)	Δ = 5.45 × 10^–4		Δ = 9.1 × 10^–5		Δ = 1.677 × 10^–5
N	20	10	30	10	42	10
(1 + x)¹⁰	Δ = 2.381 × 10^–3		Δ = 1.1 × 10^–5		Δ = 2.548 × 10^–8
N	24	10	86	10	378	10

Отсюда можно сделать следующие выводы:

1) погрешность при использовании МБСР с оператором $C{{h}_{6}}$ на несколько порядков меньше по сравнению с методом Симпсона (более 100 раз и еще больше при использовании операторов высокого порядка);

2) с ростом порядка оператора $C{{h}_{{2p}}}$ при неизменном $N$ погрешность быстро уменьшается;

3) с увеличением количества членов $N$ в ряде Фурье (1.15) погрешность также быстро уменьшается.

Алгоритм нахождения определенного интеграла от сложной функции с переменным верхним пределом

1. Считаем заданными гладкую функцию $f\left( x \right)$ и размер отрезка $a$.

2. Учитывая данные в таблицах или другие данные подобного характера, выберем значение порядка $2p$ граничной функции и число учитываемых членов $N$ в ряде Фурье.

3. Внутри отрезка $\left[ {0,a} \right]$ равномерно нанесем $N$ внутренних точек ${{x}_{k}} = {{ak} \mathord{\left/ {\vphantom {{ak} {\left( {N + 1} \right)}}} \right. \kern-0em} {\left( {N + 1} \right)}}$, $k = 1{\kern 1pt} \div {\kern 1pt} N$.

4. Вычислим значения функции ${{\varphi }_{{2p}}}\left( {{{x}_{k}}} \right)$ во внутренних расчетных точках ${{x}_{k}} = {{ak} \mathord{\left/ {\vphantom {{ak} {\left( {N + 1} \right)}}} \right. \kern-0em} {\left( {N + 1} \right)}}$, $k = 1{\kern 1pt} \div {\kern 1pt} N$. Вид функции ${{\varphi }_{{2p}}}\left( {{{x}_{k}}} \right)$ представлен в (1.5) при$p = 1$, в (1.12) при$p = 2$, в (1.16) при $p = 3$.

5. Коэффициенты Фурье $f_{{2p,m}}^{*}$ найдем по одной из формул (1.9), (1.13) или (1.17), в зависимости от значения порядка $2p$.

6. Определенный интеграл (1.2) вычислим по соответственной формуле (1.10), (1.14) или (1.18), имеющей явный и удобный для вычислений вид.

Численное сравнение метода БСР с методами Maple позволяют сделать следующие выводы.

1. Современные методы Maple вычисляют определенный интеграл с постоянными пределами от сложных функций с точностью до 300 знаков и более после запятой, тогда как метод БСР – с точностью до 35 и более знаков.

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

3. В случае интеграла с переменным верхним пределом методы Maple позволяют его вычислить только от элементарных или спецфункций, заложенных в его базу данных. Если же подынтегральная функция сложная, например $f\left( x \right) = {{x}^{2}}\sin \left( {3 + \ln \left( {1 + {{x}^{2}}} \right)} \right)$, то Maple отказывается от вычисления, тогда как метод БСР легко справляется с проблемой. Лишь бы $f\left( x \right)$ была достаточно гладкой и допускала ее представление быстрым разложением.

4. Метод быстрых синус-разложений имеет еще одно важное преимущество. Если функция задана неявно, или требуется вычислить интеграл от сложной функции с переменным верхним пределом, то Maple, если сможет в соответствии со своей базой данных, выдаст ответ в сложном виде через спецфункции. Тогда как метод БСР во всех подобных случаях выдает ответ через элементарные функции с высокой точностью. Это создает существенные удобства в дальнейших исследованиях, где используются подобные интегралы.

2. ПРЕДСТАВЛЕНИЕ НЕЯВНО ЗАДАННОЙ ФУНКЦИИ В ЯВНОМ “АНАЛИТИЧЕСКОМ” ВИДЕ И ВЫЧИСЛЕНИЕ ОПРЕДЕЛЕННОГО ИНТЕГРАЛА С ПЕРЕМЕННЫМ ВЕРХНИМ ПРЕДЕЛОМ

Решение поставленной здесь проблемы представляет определенный интерес при рассмотрении различных прикладных вопросов, а также при необходимости вычисления определенных интегралов с переменным верхним пределом или при проведении других исследований. На первый взгляд постановка самой задачи о вычислении подобного интеграла от неявно заданной функции является противоречивой. Ведь функция в явной форме отсутствует и потому как же можно вычислить подобный определенный интеграл, да еще и с переменным верхним пределом. В классической литературе и научных статьях данная проблема вообще не обсуждается. Вполне возможно, что здесь она рассматривается впервые. В этой связи вначале решим вспомогательную задачу о приближенном представлении неявно заданной функции в явном “аналитическом” виде при помощи метода быстрых разложений [5]. Это позволит вычислять $f\left( x \right)$ и ее производные до $\left( {2p + 1} \right)$-го порядка включительно в любой внутренней точке и на концах отрезка $\left[ {0,a} \right]$.

Обычно под аналитической функцией понимают функцию, которая непрерывна вместе со всеми своими производными любого конечного порядка в области определения. Ниже мы будем называть функцию аналитической, если она непрерывна вместе со всеми своими производными до некоторого заданного конечного m-го порядка включительно на данном отрезке $\left[ {0,a} \right]$ и ее можно вычислить не только в узловых точках, но и при любом $x \in \left[ {0,a} \right]$. Подобное определение необходимо в связи с применением теории быстрых разложений, так как здесь используются ряды Фурье, допускающие почленное дифференцирование некоторое ограниченное заданное число раз. Поэтому над словом «аналитическая» в дальнейшем кавычки будем опускать.

Пусть некоторая функция

(2.1)

$y = f\left( x \right) \in L_{2}^{{\left( {2p + 1} \right)}}\quad \forall x \in \left[ {0,a} \right],\quad S{\kern 1pt} {\text{*}}\left( f \right)$

задана на отрезке $x \in \left[ {0,a} \right]$ неявной формой:

(2.2)

$F\left( {x,y} \right) = 0,\quad F\left( {x,y} \right) \in {{C}^{{\left( {2p + 1} \right)}}}\left( {{{\Omega }_{\square }}} \right),\quad \left( {x,y} \right) \in {{\Omega }_{\square }},\quad 0 \leqslant x \leqslant a,\quad 0 \leqslant y \leqslant b.$

Здесь $2p + 1$ – целое число. Под $S{\kern 1pt} {\text{*}}\left( f \right)$ в (2.1) будем понимать некоторое дополнительное условие, позволяющее выделить единственное из множества возможных решений уравнения $F\left( {{{x}_{k}},{{y}_{k}}} \right) = 0$ относительно ${{y}_{k}}$ при любом заданном $0 \leqslant {{x}_{k}} \leqslant a$. Предположим, что неявная функция (2.1), определенная зависимостью (2.2), существует, единственна и удовлетворяет требованиям, указанным в (2.1). Размеры прямоугольника ${{\Omega }_{\square }} = a \times b$ в (2.2) считаем конечными и известными. Полагая $x = 0$ в (2.2) и затем $x = a$, определим значения функции на концах отрезка $\left[ {0,a} \right]$

(2.3)

${{y}_{0}} = f\left( 0 \right)\quad {\text{и}}\quad {{y}_{a}} = f\left( a \right).$

Для существования зависимости (2.1) в приближенном явном виде при использовании метода БСР необходимо наложить условия на частные производные

(2.4)

${{F}_{y}}\left( {x,y} \right) \ne 0,\quad \,\left| {{{F}_{{\underbrace {x \ldots x}_i,\underbrace {y \ldots y}_j,}}}\left( {x,y} \right)} \right| < \infty \quad \forall \left( {i + j} \right) = 0{\kern 1pt} \div {\kern 1pt} 2p + 1,\quad \left( {x,y} \right) \in {{\Omega }_{\square }}.$

Здесь $\left( {i + j} \right)$ – суммарный порядок частной производной. Для возможности представления $f\left( x \right)$ быстрым разложением вычислим производную от равенства (2.2) как от сложной функции по переменной $x$, учитывая возможную зависимость $y = f\left( x \right)$, записанную в (2.1):

(2.5)

${{F}_{x}}\left( {x,y} \right) + {{F}_{y}}\left( {x,y} \right)y{\text{'}} = 0.$

Из (2.5) получаем известное выражение для первой производной, что позволяет вычислить ее значения на концах отрезка, используемые при применении метода БСР:

(2.6)

$\begin{gathered} y{\text{'}}\left( x \right) = - {{{{F}_{x}}\left( {x,y} \right)} \mathord{\left/ {\vphantom {{{{F}_{x}}\left( {x,y} \right)} {{{F}_{y}}\left( {x,y} \right)}}} \right. \kern-0em} {{{F}_{y}}\left( {x,y} \right)}}, \\ y{\text{'}}\left( 0 \right) = - {{{{F}_{x}}\left( {0,{{y}_{0}}} \right)} \mathord{\left/ {\vphantom {{{{F}_{x}}\left( {0,{{y}_{0}}} \right)} {{{F}_{y}}\left( {0,{{y}_{0}}} \right)}}} \right. \kern-0em} {{{F}_{y}}\left( {0,{{y}_{0}}} \right)}},\quad y{\text{'}}\left( a \right) = - {{{{F}_{x}}\left( {a,{{y}_{a}}} \right)} \mathord{\left/ {\vphantom {{{{F}_{x}}\left( {a,{{y}_{a}}} \right)} {{{F}_{y}}\left( {a,{{y}_{a}}} \right)}}} \right. \kern-0em} {{{F}_{y}}\left( {a,{{y}_{a}}} \right)}}. \\ \end{gathered} $

После повторного дифференцирования (2.5), как сложную функцию по $x$, найдем $y{\text{''}}\left( x \right)$, необходимую для построения граничной функции ${{M}_{{2p}}}\left( x \right)$:

(2.7)

${{F}_{{xx}}}\left( {x,y} \right) + 2{{F}_{{xy}}}\left( {x,y} \right)y{\text{'}} + {{F}_{{y\,y}}}\left( {x,y} \right){{y}^{{'2}}} + {{F}_{y}}\left( {x,y} \right)y{\text{''}} = 0.$

Из (2.7) вычислим значения второй производной на концах отрезка

(2.8)

$\begin{gathered} y{\text{''}}\left( x \right) = - {{\left( {{{F}_{{xx}}}\left( {x,y} \right) + 2{{F}_{{xy}}}\left( {x,y} \right)y{\text{'}} + {{F}_{{y\,y}}}\left( {x,y} \right){{y}^{{'2}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{F}_{{xx}}}\left( {x,y} \right) + 2{{F}_{{xy}}}\left( {x,y} \right)y{\text{'}} + {{F}_{{y\,y}}}\left( {x,y} \right){{y}^{{'2}}}} \right)} {{{F}_{y}}\left( {x,y} \right)}}} \right. \kern-0em} {{{F}_{y}}\left( {x,y} \right)}}, \\ y{\text{''}}\left( 0 \right) = - {{\left( {{{F}_{{xx}}}\left( {0,{{y}_{0}}} \right) + 2{{F}_{{xy}}}\left( {0,{{y}_{0}}} \right)y{\text{'}} + {{F}_{{yy}}}\left( {0,{{y}_{0}}} \right){{y}^{{'2}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{F}_{{xx}}}\left( {0,{{y}_{0}}} \right) + 2{{F}_{{xy}}}\left( {0,{{y}_{0}}} \right)y{\text{'}} + {{F}_{{yy}}}\left( {0,{{y}_{0}}} \right){{y}^{{'2}}}} \right)} {{{F}_{y}}\left( {0,{{y}_{0}}} \right)}}} \right. \kern-0em} {{{F}_{y}}\left( {0,{{y}_{0}}} \right)}}, \\ y{\text{''}}\left( a \right) = - {{\left( {{{F}_{{xx}}}\left( {a,{{y}_{a}}} \right) + 2{{F}_{{xy}}}\left( {a,{{y}_{a}}} \right)y{\text{'}} + {{F}_{{yy}}}\left( {a,{{y}_{a}}} \right){{y}^{{'2}}}} \right)} \mathord{\left/ {\vphantom {{\left( {{{F}_{{xx}}}\left( {a,{{y}_{a}}} \right) + 2{{F}_{{xy}}}\left( {a,{{y}_{a}}} \right)y{\text{'}} + {{F}_{{yy}}}\left( {a,{{y}_{a}}} \right){{y}^{{'2}}}} \right)} {{{F}_{y}}\left( {a,{{y}_{a}}} \right)}}} \right. \kern-0em} {{{F}_{y}}\left( {a,{{y}_{a}}} \right)}}. \\ \end{gathered} $

Подобным образом находятся все четные производные ${{y}^{{\left( {2p} \right)}}}\left( 0 \right)$, ${{y}^{{\left( {2p} \right)}}}\left( a \right)$, используемые в конструкции граничной функции ${{M}_{{2p}}}\left( x \right)$.

Для построения ${{M}_{{2p}}}\left( x \right)$ произвольного четного порядка приведем достаточно удобную рекуррентную формулу. Для этого запишем ${{M}_{{2p}}}\left( x \right)$ суммой

(2.9)

${{M}_{{2p}}}\left( x \right) = \sum\limits_{m = 0}^p {\left( {{{f}^{{\left( {2m} \right)}}}\left( 0 \right){{P}_{{2m}}}\left( x \right) + {{f}^{{\left( {2m} \right)}}}\left( a \right){{Q}_{{2m}}}\left( x \right)} \right)} ,\quad x \in \left[ {0,a} \right].$

Здесь ${{f}^{{\left( {2m} \right)}}}\left( 0 \right)$, ${{f}^{{\left( {2m} \right)}}}\left( a \right)$ – значения четных производных $2m$-го порядка функции $f\left( x \right)$ на концах отрезка $\left[ {0,a} \right]$, ${{P}_{{2m}}}\left( x \right)$, ${{Q}_{{2m}}}\left( x \right)$ – полиномы, для построения которых найдем рекуррентную формулу следующим образом. Вначале запишем граничную функцию ${{M}_{0}}\left( x \right)$ нулевого порядка

(2.10)

${{M}_{0}}\left( x \right) = f\left( 0 \right)\left( {1 - \frac{x}{a}} \right) + f\left( a \right)\frac{x}{a},\quad {{P}_{0}}\left( x \right) = \left( {1 - \frac{x}{a}} \right),\quad {{Q}_{0}}\left( x \right) = \frac{x}{a}.$

В (2.10) даны полиномы ${{P}_{0}}\left( x \right)$ и ${{Q}_{0}}\left( x \right)$ нулевого порядка. Рекуррентные формулы для ${{P}_{{2m}}}\left( x \right)$, ${{Q}_{{2m}}}\left( x \right)$ при $m = 1{\kern 1pt} \div {\kern 1pt} p$ записываются следующими интегральными выражениями: если известны полиномы ${{P}_{{2m - 2}}}\left( x \right)$, ${{Q}_{{2m - 2}}}\left( x \right)$, то полиномы порядка ${{P}_{{2m}}}\left( x \right)$, ${{Q}_{{2m}}}\left( x \right)$ вычислим через нижеприведенные неопределенные интегралы (2.11) с учетом их значений на концах отрезка $\left[ {0,a} \right]$:

(2.11)

$\begin{gathered} {{P}_{{2m}}}\left( x \right) = \int {\left( {\int {{{P}_{{2m - 2}}}\left( x \right)dx} } \right)dx} ,\quad {{P}_{{2m}}}\left( 0 \right) = {{P}_{{2m}}}\left( a \right) = 0, \\ {{Q}_{{2m}}}\left( x \right) = \int {\left( {\int {{{Q}_{{2m - 2}}}\left( x \right)dx} } \right)} dx,\quad {{Q}_{{2m}}}\left( 0 \right) = {{Q}_{{2m}}}\left( a \right) = 0,\quad m = 1{\kern 1pt} \div {\kern 1pt} p. \\ \end{gathered} $

Теперь можно представить неявно заданную функцию $f\left( x \right)$ с оператором $C{{h}_{{2p}}}$ порядка $2p$ быстрым разложением:

(2.12)

$f\left( x \right) = {{M}_{{2p}}}\left( x \right) + \sum\limits_{m = 1}^N {{{f}_{{2p,m}}}\sin m\pi \frac{x}{a}} ,$

где ${{f}_{{2p,m}}}$ – коэффициенты ряда Фурье для разности $f\left( x \right) - {{M}_{{2p}}}\left( x \right)$, $N$ – число учитываемых членов ряда, ${{M}_{{2p}}}\left( x \right)$ – граничная функция. Ее конкретный вид при $p = 1$ приводится в (1.3), в (1.11) при $p = 2$, в (1.15) при $p = 3$.

Разложение (2.12) обеспечивает быструю сходимость ряда Фурье. Хотя $f\left( x \right)$ и задана неявно, но значения всех четных производных ${{f}^{{\left( {2m} \right)}}}\left( 0 \right)$, ${{f}^{{\left( {2m} \right)}}}\left( a \right)$, $m = 0{\kern 1pt} \div {\kern 1pt} p$, для определения ${{M}_{{2p}}}\left( x \right)$ можно считать найденными. При $p = 1$ данные производные представлены формулами (2.8). Для завершения построения ${{M}_{{2p}}}\left( x \right)$ остается найти коэффициенты Фурье ${{f}_{{2,m}}}$, но их интегральные выражения неизвестны, так как зависимость $f\left( x \right)$ явно не дана. Поэтому в дальнейшем используем поточечный способ.

Для решения задачи на отрезок $\left[ {0,a} \right]$ равномерно нанесем $N$ внутренних точек ${{x}_{k}}$ :

(2.13)

${{x}_{k}} = {{ak} \mathord{\left/ {\vphantom {{ak} {\left( {N + 1} \right)}}} \right. \kern-0em} {\left( {N + 1} \right)}},\quad k = 1{\kern 1pt} \div {\kern 1pt} N.$

Будем полагать, что функция $F\left( {x,y} \right)$ такова, что при выполнении условий $S{\kern 1pt} {\text{*}}\left( f \right)$ в (2.1) неизвестные ${{y}_{k}}$ найдем из алгебраических уравнений

(2.14)

$F\left( {{{x}_{k}},{{y}_{k}}} \right) = 0,\quad {{x}_{k}} = {{ak} \mathord{\left/ {\vphantom {{ak} {\left( {N + 1} \right)}}} \right. \kern-0em} {\left( {N + 1} \right)}},\quad k = 1{\kern 1pt} \div {\kern 1pt} N.$

Используя формулы тригонометрической интерполяции [11], выражение (2.12) заменим на следующее:

(2.15)

$f\left( x \right) = {{M}_{2}}\left( x \right) + \sum\limits_{m = 1}^N {f_{{2,m}}^{*}sinm\pi \frac{x}{a}} ,$

где коэффициенты $f_{{2,m}}^{*}$ определяются из линейной системы, полученной из (2.15) при $x = {{x}_{k}}$:

(2.16)

$f\left( {{{x}_{k}}} \right) = {{y}_{k}} = {{M}_{2}}\left( {{{x}_{k}}} \right) + \sum\limits_{m = 1}^N {f_{{2,m}}^{*}sinm\pi \frac{{{{x}_{k}}}}{a}} ,\quad k = 1{\kern 1pt} \div {\kern 1pt} N.$

Здесь величины ${{y}_{k}}$ считаем уже найденными из уравнений (2.14). Подобным же образом можно построить ${{M}_{{2p}}}\left( {{{x}_{k}}} \right)$ и вычислить коэффициенты $f_{{2p,m}}^{*}$. Затем по формуле, аналогичной формулам (1.10), (1.14), (1.18), в зависимости от величины порядка $2p$, вычисляется интеграл с переменным верхним пределом от неявно заданной функции.

В качестве примера рассмотрим следующую задачу.

Задача. Найти определенный интеграл (1.2) с переменным верхним пределом от следующей неявно заданной функции:

(2.17)

${{x}^{2}} + {{y}^{2}} = 2\sin \left( {xy + 0.9} \right) + 4,\quad x \in \left[ {0,1} \right],\quad y \geqslant 0.$

Полагая $x = {{x}_{k}} = {k \mathord{\left/ {\vphantom {k {11}}} \right. \kern-0em} {11}}$, $k = 0{\kern 1pt} \div {\kern 1pt} 11$, вычислим $\left\{ {{{y}_{k}}} \right\}$ в (2.17).

Возьмем $2p = 6$, так как оператор $C{{h}_{6}}$ шестого порядка дает очень высокую точность. Тогда из (2.17) надо вычислить производные ${{y}^{{\left( n \right)}}}\left( 0 \right)$, ${{y}^{{\left( n \right)}}}\left( 1 \right)$, $n = 1 \div 6$ на концах отрезка $\left[ {0,\;1} \right]$.

Из (1.15) вычислим значения ${{M}_{6}}\left( {{{x}_{k}}} \right)$ в точках $x = {{x}_{k}}$ и затем по формуле ${{\varphi }_{6}}\left( {{{x}_{k}}} \right) = f\left( {{{x}_{k}}} \right) - {{M}_{6}}\left( {{{x}_{k}}} \right)$ найдем ${{\varphi }_{6}}\left( {{{x}_{k}}} \right)$.

Теперь из (1.17) можно определить коэффициенты $f_{{6,m}}^{*}$. После подстановки в (1.15) значений четных производных ${{y}^{{\left( {2n} \right)}}}\left( 0 \right)$, ${{y}^{{\left( {2n} \right)}}}\left( 1 \right)$, $n = 0,\;1,\;2,\;3$, и коэффициентов $f_{{6,m}}^{*}$ будем иметь в (2.17) неявно заданную функцию приближенно в явном виде. Оценка погрешности $\Delta $ вычислялась по формуле

(2.18)

$\Delta = \sup \left| {y\left( x \right) - {{M}_{6}}\left( x \right) - \sum\limits_{m = 1}^{10} {f_{{6,m}}^{*}sinm\pi x} } \right|,\quad x \in \left[ {0,1} \right] \Rightarrow \Delta \left( {x = 0.04} \right) = 4.547 \times {{10}^{{ - 8}}},$

где значения $y\left( x \right)$ вычислялись из выражения (2.17).

Из (1.18) получаем искомую зависимость интеграла (1.2) от переменной $x \in \left[ {0,\;1} \right]$ сразу на всем отрезке в явном виде с высокой точностью.

Приведенные численные примеры показывают очень высокую точность и быструю сходимость метода. Кроме предложенного здесь метода БСР, точно также может быть развит метод быстрых косинус-разложений (МБКР). Оба метода можно применять для решения, например, уравнений Вольтерра I и II родов, для вычисления многомерных интегралов с переменными пределами и других сложных систем интегродифференциальных уравнений.

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

Крылов В.И. Приближенное вычисление интегралов. М.: Наука, 1967. С. 377–432.
Габдулхаев Б.Г. Квадратурные формулы наивысшей тригонометрической степени точности и их приложения // Известия ВУЗов. Математика. 2007. № 7 (542). С. 28–41.
Боголюбский А.И., Скороходов С.Л., Христофоров Д.В. Быстрое вычисление эллиптических интегралов и их обобщений // Ж. вычисл. матем. и матем. физ. 2005. Т. 45. № 11. С. 1938–1953.
Гласко А.В. Повышение устойчивости метода вычисления коэффициентов абсолютно сходящегося ряда, приближающего функциональный интеграл // Ж. вычисл. матем. и матем. физ. 2000. Т. 40. № 1. С. 30–34.
Абилов В.А., Абилова Ф.В., Керимов М.К. О некоторых оценках для интегрального преобразования Фурье-Бесселя в пространстве ${{L}_{2}}\left( R \right)$ // Ж. вычисл. матем. и матем. физ. 2009. Т. 49. № 7. С. 1158–1166.
Чернышов А.Д. Метод быстрых разложений для решения нелинейных дифференциальных уравнений // Ж. вычисл. матем. и матем. физ. 2014. Т. 54. № 1. С. 13–24.
Жидков Е.П., Лобанов Ю.Ю., Рушай В.Д. Применение метода Ромберга для повышения точности вычисления кратных интегралов // Ж. вычисл. матем. и матем. физ. 2009. Т. 49. № 2. С. 232–240.
Clenshaw C.W., Curtis A.R. A method for numerical integration in an automatic computer // Numer. Math. 1960. Bd. 2. S. 197–205.
Чернышов А.Д., Попов В.М., Шахов А.С., Горяйнов В.В. Повышенная точность решения задачи о контактном термосопротивлении между сжатыми шарами методом быстрых разложений // Тепловые процессы в технике. 2014. Т. 6. № 4. С. 179–191.
Чернышов А.Д., Горяйнов В.В. Применение метода быстрых разложений для расчета траекторий космического корабля // Изв. вузов. Авиационная техника. 2015. № 2. С. 41–47.
Isaev V.I., Shapeev V.P. High-Accuracy of the Collocations and Least Squares Method for the Numerical Solution of the Navies-Stokes Equations // Comp. Math. and Math. Phys. 2010. V. 50. № 10. P. 1670–1681.
Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987. С. 166.
Чернышов А.Д. Решение нелинейного уравнения теплопроводности для криволинейной области методом быстрых разложений // Инженерно-физ. ж. Минск. 2018. Т. 91. № 2. С. 456.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал