Журнал вычислительной математики и математической физики, 2022, T. 62, № 4, стр. 580-596

ВВЕДЕНИЕ

Пусть ${\mathbf{C}}[a;b]$ – пространство функций, непрерывных на отрезке $[a;b]$. Классическим подходом к аппроксимации функций из ${\mathbf{C}}[a;b]$ является использование многочленов. Этот подход опирается на известную теорему Вейерштрасса о том, что всякую такую функцию можно с любой заданной точностью в метрике пространства ${\mathbf{C}}[a;b]$ аппроксимировать многочленом. Получить наилучшее представление непрерывной функции многочленом позволяет теория П.Л. Чебышева. Уместно упомянуть здесь теоремы типа С.М. Никольского, которые позволяют при фиксированном $s > 0$ указать класс непрерывных функций, каждая из которых для всякого $n \in \mathbb{N}$ допускает равномерное на $[a;b]$ приближение с точностью порядка ${{n}^{{ - s}}}$ полиномами степени $n$ (подробнее см. [1, гл. VI]). В случае, когда аппроксимируемая функция является достаточно гладкой всюду, за исключением одной или нескольких точек, более эффективно использование рациональных многочленов вида ${{\mathcal{P}}_{n}}(x){\text{/}}{{\mathcal{Q}}_{n}}(x)$, где ${{\mathcal{P}}_{n}}(x)$, ${{\mathcal{Q}}_{n}}(x)$ – алгебраические многочлены степени $n$ (см. [1, гл. VI]). (О других способах аппроксимации см. в [2, гл. 2, ${{\S}}$ 2.7, с. 66], [3, гл. 10].)

В теории обработки сигналов, изображений и т.д. используются также подходы, основанные, в частности, на теории вейвлетов (всплесков). Эти подходы оказались весьма эффективными в смысле возможности сжатия и преобразования информации (см. [4], [5]). Конкретное семейство вейвлетов порождается сдвигами и сжатиями из так называемого материнского вейвлета – некоторой функции, удовлетворяющей специальным условиям. Типичным примером является материнский вейвлет “мексиканская шляпа” (подробнее см. [4, гл. 3]). Сдвиги и сжатия квадратичной экспоненты (т.е. функции Гаусса) – это тоже квадратичные экспоненты. В [6] было показано, что материнский вейвлет “мексиканская шляпа” на любом конечном отрезке можно сколь угодно точно в равномерной метрике приблизить линейной комбинацией в точности двух квадратичных экспонент. Поэтому на фиксированном отрезке аппроксимация квадратичными экспонентами сопоставима по эффективности с аппроксимацией вейвлетом “мексиканская шляпа”. В [7] доказано, что любой многочлен степени $n$ на любом конечном отрезке можно сколь угодно точно в равномерной метрике приблизить линейной комбинацией не более $(n + 1)$ квадратичной экспоненты. Отметим, что аппроксимация квадратичными экспонентами строится неоднозначно и допускает бесконечное количество вариантов и, вообще, имеет ряд общих черт с представлением фреймами (см. [5, ${{\S}}$ 1.8]). Что касается указанной неоднозначности, то для инженерных приложений она бывает полезна (см. [5, ${{\S}}$ 1.8]).

В представленном выше обзоре речь шла об аппроксимации без дополнительных условий. Однако бывают ситуации (например, при решении некоторых задач математической статистики, см. [8]), когда требуется строить монотонную аппроксимацию для монотонной непрерывной функции. В качестве содержательного примера такой ситуации укажем следующую. Прежде всего, приведем формулировку модификации теоремы А.Н. Колмогорова в стиле работы [9, разд. “Введение”, “Построение внутренних функций”], опирающейся в свою очередь на [10].

Лемма 1. При произвольно заданных числах: натуральном $n \geqslant 2$, алгебраическом $\lambda > 0$ степени $\deg \lambda \geqslant n$, натуральном $\gamma \geqslant 2(n + 1)$ и $\sigma = 1{\text{/}}(\gamma - 1)$, существует действительная непрерывная монотонно возрастающая функция одного переменного $\Psi (y)$, $y \geqslant 0$, такая, что каждая определенная на $n$-мерном единичном кубе ${{E}^{n}}$ непрерывная действительная функция $f({{x}_{1}}, \ldots ,{{x}_{n}})$ представима в виде

где семейство действительных непрерывных функций одного переменного $\{ {{F}_{q}}(y),q = \overline {0,2n} \} $ определяется по функции $f$.

Способ из [9] построения функции $\Psi $ (в виде бесконечного ряда) сложный и многоэтапный. Поэтому его практическое использование затруднительно. При этом функции $\Psi $, ${{F}_{q}}$ определяются неоднозначно.

В [11] было предложено аппроксимировать непрерывные функции многих переменных на основе леммы 1 и аппроксимации функций $\Psi (y)$, ${{F}_{q}}(y)$ линейными комбинациями квадратичных экспонент (там же см. теоретическое обоснование этого подхода, описание его преимуществ и в случае размерности $2$ результаты численных экспериментов). Как видно из этих экспериментов, для различных функций $f(x)$ аппроксимация функции, играющей роль $\Psi (y)$ в выражении ${{G}_{q}}(x)$ из леммы 1, может выглядеть по-разному, быть не только возрастающей, но и убывающей (что соответствует умножению аргумента функций ${{F}_{q}}(y)$ на $( - 1)$; возможен сдвиг или растяжение этой функции); условие монотонности может оказаться нарушенным на некоторых участках. Указанное нарушение монотонности, вообще говоря, не удивительно, поскольку аппроксимируемая функция $\Psi (y)$ входит в аргумент внешних функций ${{F}_{q}}(y)$, которые одновременно с ней подвергаются аппроксимации. Может оказаться и так, что поведение функции $\Psi $ на некоторых участках не оказывает влияния на суперпозиции внешних и внутренних функций. С другой стороны, нарушение монотонности, по крайней мере, по общим соображениям может приводить к замедлению или затруднению процесса сходимости при оптимизации квадратичной невязки, либо минимизации функционала в задаче бесконечномерной оптимизации, при конечномерной аппроксимации искомой функции, по параметрам аппроксимации. Естественно возникает вопрос о разработке эффективного способа монотонной аппроксимации функции $\Psi (y)$. Желание получить ответ на этот вопрос как раз и послужило исходной мотивацией написания данной статьи.

В настоящей работе для непрерывных монотонных функций, заданных на конечном отрезке $[ - b;b]$, строится монотонная аппроксимация с любой заранее заданной точностью в метрике пространства ${\mathbf{C}}[ - b;b]$ с помощью сдвигов и сжатий функции (интеграла) Лапласа. Этот способ построения аппроксимации основан, в общих чертах, на следующих соображениях.

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

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

3. Интеграл от квадрата линейной комбинации квадратичных экспонент можно представить в виде линейной комбинации (некоторого специального шаблона) сдвигов и сжатий интеграла Лапласа. Отсюда следует вывод, что гладкую монотонную функцию можно эффективно аппроксимировать специальной линейной комбинацией сдвигов и сжатий интеграла Лапласа, обладающей свойством монотонности.

4. Для непрерывной монотонной функции всегда можно построить бесконечно дифференцируемую монотонную аппроксимацию (в виде многочлена). Таким образом, тот же вывод можно сделать и в отношении непрерывных монотонных функций. Для краткости построенную таким образом аппроксимацию будем обозначать $Q(x)$.

По сравнению с применением для аппроксимации непосредственно квадратичных экспонент, здесь можно усмотреть лишь некоторое неудобство, связанное с тем, что интеграл Лапласа не выражается в элементарных функциях. Но с практической точки зрения, при использовании любого математического пакета, квадратичная экспонента (так же, как и интеграл Лапласа) вычисляется все равно приближенно. Кроме того, существуют математические пакеты (в частности, MATLAB), в которых реализована функция ошибки (error function): $\operatorname{erf} (x) = 2{\text{/}}\sqrt \pi \int_0^x {{{e}^{{ - {{t}^{2}}}}}dt} $. Ясно, что интеграл Лапласа ${{\widehat \Phi }_{0}}(x) = 1{\text{/}}\sqrt {2\pi } \int_0^x {{{e}^{{ - {{t}^{2}}/2}}}dt} = 1{\text{/}}2\operatorname{erf} (x{\text{/}}\sqrt 2 )$. Аналогично функции Гаусса, интеграл Лапласа обладает свойством локализации: достаточно знать его значения на $[0;5]$, поскольку это нечетная функция, а при $x > 5$ значения ее практически не отличаются от $1{\text{/}}2$; на $[0;5]$ из математической статистики известны таблицы значений, которые можно использовать в готовом виде, не тратя времени на вычисления (в точках между узлами сетки можно использовать линейную интерполяцию). Наконец, в разд. 3 мы показываем, что интеграл Лапласа можно с высокой степенью точности аппроксимировать суммой линейной функции и линейной комбинации двух-трех квадратичных экспонент, и указываем конкретные значения ее коэффициентов. Представляется, что этот последний результат интересен и сам по себе, поскольку интеграл Лапласа играет важную роль при решении многих прикладных задач, связанных с теорией вероятностей, математической статистикой и вычислительной математикой. При повышении количества квадратичных экспонент в линейной комбинации точность такой аппроксимации (обозначим ее ${{\widehat \Phi }_{{0,\varepsilon }}}(x)$), разумеется, возрастает. При этом точность аппроксимации оценивается по разности вторых производных функций ${{\widehat \Phi }_{0}}(x)$ и ${{\widehat \Phi }_{{0,\varepsilon }}}(x)$ в пространстве ${{L}_{2}}$ и, соответственно, первых производных в пространстве ${\mathbf{C}}$. Более того, мы показываем, что при достаточной малости отклонения $\left| {\widehat \Phi _{0}^{'}(x) - \widehat \Phi _{{0,\varepsilon }}^{'}(x)} \right|$ в метрике ${\mathbf{C}}$ имеет место устойчивость свойства монотонности аппроксимации $Q(x)$ при замене в ней интеграла Лапласа на ${{\widehat \Phi }_{{0,\varepsilon }}}(x)$. Кроме того, мы представляем результаты численных экспериментов по аппроксимации изучаемым способом различных монотонных функций: четырех гладких и одной непрерывной кусочно линейной. Во всех указанных случаях эффективность изучаемого способа аппроксимации оказалась достаточно высокой, а устойчивость свойства монотонности аппроксимации подтвердилась экспериментально.

Сделаем краткий обзор по работам, посвященным аппроксимации монотонных функций. По-видимому, впервые эта проблема была поставлена в [12], [13]. В [14] был получен результат об аппроксимации непрерывной возрастающей функции возрастающими многочленами. В дальнейшем этот результат уточнялся, в основном, на случай монотонных функций различной степени гладкости (см., например, [15]–[17]). Вопрос о достижимости точности аппроксимации непрерывной неубывающей функции $f(x)$ неубывающим многочленом ${{\mathcal{P}}_{n}}(x)$ на уровне порядка ${{\omega }_{k}}(f,{{n}^{{ - 1}}})$, $k = 2,3,4$, где ${{\omega }_{k}}$ – модуль непрерывности порядка $k$, исследовался в [18]. Отметим, что процедура построения аппроксимирующих многочленов была достаточно сложной и многоэтапной. Аппроксимация достаточно гладких $s$-монотонных функций изучалась в [19]. Функция $f:I = ( - 1;1) \to \mathbb{R}$ называется $s$-монотонной для целого $s \geqslant 0$, если $\Delta _{\tau }^{s}f(t) \geqslant 0$ для всех $t \in I$, $\tau > 0$. Здесь $\Delta _{\tau }^{s}f(t) = \sum\nolimits_{k = 0}^s {{{{( - 1)}}^{{s - k}}}C_{s}^{k}f(t + k\tau )} $ при $(t,t + s\tau ) \subset I$ (иначе равная нулю) – $s$-я конечная разность с шагом $\tau $. В [20] изучалась проблема ${{L}_{1}}$-аппроксимации монотонных функций $d$ переменных по $n$ вычислениям функции методом Монте-Карло.

Упомянутые работы [14]–[20] носили, между тем, чисто теоретический характер; результаты численных экспериментов в них представлены не были (строго говоря, в [20] есть один пример с численными расчетами, но, как там сказано, он имеет лишь теоретическое значение).

Подчеркнем, что автора проблема аппроксимации функций интересует с точки зрения приближения бесконечномерных задач оптимизации конечномерными. В такой ситуации значения искомой функции в контрольных точках нам неизвестны, поскольку заранее неизвестна сама функция (и в этом состоит принципиальное отличие от обычной задачи аппроксимации). В этом случае неизвестные параметры аппроксимации определяются путем минимизации значения целевого функционала, вычисляемого на конечномерном множестве аппроксимаций $\Omega $, если они являются допустимыми по смыслу задачи либо на суперпозициях, содержащих функции класса $\Omega $ и принимающих допустимые значения. Указанное принципиальное отличие следует учитывать при сравнении различных подходов к аппроксимации (подробнее об этом уже было сказано в [11]): здесь главное – это достижение приемлемой точности аппроксимации произвольной функции уже при небольшом количестве параметров. Как показано в [7], [11] (там же см. результаты численных экспериментов), достаточно эффективными (в указанном смысле) оказываются аппроксимации с помощью квадратичных экспонент. Однако вопрос сохранения монотонности аппроксимации там не исследовался.

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

Для теории вероятностей и математической статистики важное значение имеет так называемый интеграл Лапласа

Функцию ${{\Phi }_{0}}(x)$ мы, соответственно, называем функцией Лапласа. Эта функция будет играть существенную роль в наших построениях.

Далее, называя непрерывную функцию возрастающей на отрезке, мы имеем в виду, что она допускает продолжение до непрерывной строго возрастающей функции на некотором объемлющем интервале. Называя гладкую функцию возрастающей на отрезке, мы имеем в виду, что ее производная строго положительна на всем отрезке. Следующее утверждение реализует соображения 1–3 построения монотонной аппроксимации гладких монотонных функций, описанные во Введении.

Теорема 1. Пусть $b > 0$, $\Psi \in {{{\mathbf{C}}}^{1}}[ - b;b]$ – неубывающая функция. Тогда для любого $\varepsilon > 0$ найдутся константа $C$ и непрерывная функция $r(t)$, $t \in [ - b;b]$, а также параметры $\nu \in \mathbb{N}$, $(\alpha ,\beta ,\gamma ) \in {{\mathbb{R}}^{{3\nu }}}$ такие, что

Всякая функция видa (1) не убывает на $\mathbb{R}$. В случае возрастающей функции $\Psi (x)$ функцию $Q(x)$ тоже можно считать возрастающей на $[ - b;b]$ (при некотором наборе параметров).

Следующее утверждение реализует соображение 4 построения монотонной аппроксимации непрерывных монотонных функций из тех, которые были описаны во Введении.

Теорема 2. Пусть $b > 0$, $\Psi \in {\mathbf{C}}[ - b;b]$ – возрастающая функция. Тогда для любого $\varepsilon > 0$ найдутся константа $C$ и непрерывная функция $R(x)$, а также параметры $\nu \in \mathbb{N}$, $(\alpha ,\beta ,\gamma ) \in {{\mathbb{R}}^{{3\nu }}}$ такие, что

где функция $Q(x)$ имеет вид $(1)$ с параметрами $(2)$ и является возрастающей на $[ - b;b]$, ${\text{|}}R(x){\kern 1pt} {\text{|}} \leqslant \varepsilon $ $\forall x \in [ - b;b]$. Всякая функция вида $(1)$ не убывает на $\mathbb{R}$.

Доказательства теорем 1, 2 представлены в следующем разделе. Аппроксимация интеграла Лапласа с помощью квадратичных экспонент обсуждается в разд. 3. Как видно из теоремы 1, выражение (1) можно заменить следующим:

В этом случае $C = \Psi (0)$. Выражение (3) более удобно, если предполагается использовать производные по параметрам. Примем обозначения Нетрудно заметить, что ${{\beta }_{{ii}}} = {{\beta }_{i}}$, ${{\gamma }_{{ii}}} = {{\gamma }_{i}}{\text{/}}2$, ${{\sigma }_{{ii}}} = 0$, откуда

Следующее утверждение необходимо при использовании методов первого порядка для минимизации квадратичной невязки между аппроксимацией $Q(x)$ вида (3) и аппроксимируемой функцией.

Теорема 3. Для функции $(3)$ имеем

Доказательство теоремы 3 представлено в следующем разделе.

Теоремы 1–3 являются теоретической основой алгоритма построения монотонной аппроксимации непрерывных монотонных функций и состоящего в минимизации квадратичной невязки (по заданному набору контрольных точек) между шаблоном вида (1) или, что эквивалентно, вида (3) и аппроксимируемой функцией. При его практической реализации на компьютере возникает вопрос о конкретном способе вычисления значений интеграла Лапласа в формулах (1) или (3). Как уже было отмечено во Введении, интеграл Лапласа не выражается в элементарных функциях. Поэтому при расчетах на компьютере приходится использовать его аппроксимацию. При этом возможно нарушение монотонности полученной компьютерной реализации $P(x)$ аппроксимации $Q(x)$. Из общих соображений ясно, что это нарушение будет достаточно незначительным (и в этом смысле приемлемым) при условии, что аппроксимация интеграла Лапласа достаточно точна. Говоря о приемлемости, можно иметь в виду, в частности, что монотонность аппроксимации требуется нам для того, чтобы облегчить процесс минимизации невязки и повысить шансы на получение более точной аппроксимации непрерывной монотонной функции. Соответственно, если отклонения от монотонности аппроксимации достаточно невелики, то это существенным образом отсекает посторонние “рыскания” процесса минимизации невязки, как бы направляя его “в нужное русло”, и в практическом плане это уже достаточно ценное достижение. Вопросы построения высокоточной аппроксимации интеграла Лапласа обсуждаются далее в разд. 3. Естественно, что при построении аппроксимации интеграла Лапласа необходимо иметь такие способы, которые позволили бы нам обеспечить сохранение свойства монотонности для реализации $P(x)$ аппроксимации $Q(x)$.

Следующее утверждение естественно понимать как теорему об устойчивой монотонности аппроксимирующей функции.

Теорема 4. Пусть (при заданном наборе параметров) функция $Q(x)$, определяемая формулой $(1)$, является возрастающей на $[ - b;b]$. Тогда при всех достаточно малых $\varepsilon > 0$ функция ${{P}_{\varepsilon }}(x)$, определяемая аналогичной формулой при замене интеграла Лапласа ${{\Phi }_{0}}(x)$ его гладкой аппроксимацией ${{\Phi }_{{0,\varepsilon }}}(x)$, удовлетворяющей условиям

остается возрастающей на $[ - b;b]$.

Доказательство теоремы 4 представлено в следующем разделе.

Описываемый в разд. 3 способ IV аппроксимации интеграла Лапласа (основанный на аппроксимации $\Phi _{0}^{{''}}(x)$ в пространстве ${{L}_{2}}[0;5]$ линейной комбинацией вторых производных квадратичных экспонент) позволяет обеспечить выполнение условий теоремы 4. На самом деле при расчетах на компьютере вместо числовой оси мы всегда имеем дело с ее дискретизацией в виде конечного набора точек. Да и кроме того, нет смысла повышать точность расчетов выше некоторой неустранимой (или диктуемой какими-то иными причинами) погрешности вычислений. В этом смысле всегда присутствует некий уровень неразличимости точек на числовой оси. В соответствии с этим будем говорить, что функция $P(x)$ приближенно возрастает на $[a;b]$ на уровне неразличимости $\delta > 0$, если для любого разбиения $a = {{x}_{0}} < {{x}_{1}} < \ldots < {{x}_{n}} = b$ при условии $\Delta {{x}_{i}} = {{x}_{i}} - {{x}_{{i - 1}}} \geqslant \delta $ имеем $P({{x}_{i}}) > P({{x}_{{i - 1}}})$ для всех $i = \overline {1,n} $.

Теорема 5. Пусть (при заданном наборе параметров) функция $Q(x)$, определяемая формулой (1), является возрастающей на $[ - b;b]$. Тогда при заданном уровне неразличимости $\delta $ и всех достаточно малых $\varepsilon > 0$ функция ${{P}_{\varepsilon }}(x)$, определяемая аналогичной формулой при замене интеграла Лапласа ${{\Phi }_{0}}(x)$ его гладкой аппроксимацией ${{\Phi }_{{0,\varepsilon }}}(x)$, удовлетворяющей условиям

будет приближенно возрастающей на $[ - b;b]$.

Доказательство теоремы 5 представлено в следующем разделе.

Описываемый в разд. 3 способ III аппроксимации интеграла Лапласа (основанный на аппроксимации $\Phi _{0}^{'}(x)$ в пространстве ${{L}_{2}}[0;5]$ линейной комбинацией первых производных квадратичных экспонент) позволяет обеспечить выполнение условий теоремы 5.

2. ДОКАЗАТЕЛЬСТВО РЕЗУЛЬТАТОВ ОБ АППРОКСИМАЦИИ С ПОМОЩЬЮ ИНТЕГРАЛА ЛАПЛАСА

Лемма 2. Имеем ${{\left( {\sum\nolimits_{i = 1}^n {{{a}_{i}}} } \right)}^{2}} = \sum\nolimits_{i = 1}^n {a_{i}^{2}} + 2\sum\nolimits_{i < j}^{} {{{a}_{i}}{{a}_{j}}} $ для любых ${{a}_{i}} \in \mathbb{R}$, $i = \overline {1,n} $.

Доказательство легко проводится математической индукцией.

Доказательство теоремы 1. Очевидно, что

В свою очередь, ${{\psi }_{0}}(t) = {{\psi }^{2}}(t)$, где $\psi (t) = \sqrt {{{\psi }_{0}}(t)} $ – непрерывная функция. Согласно [7], функция $\psi (t)$ с любой степенью точности в метрике ${\mathbf{C}}[ - b;b]$ представима линейной комбинацией квадратичных экспонент где ${{\varepsilon }_{1}} > 0$ – произвольно фиксированное число. Положим $r(t) = 2\{ \psi (t) - \sigma (t)\} \sigma (t) + {{\sigma }^{2}}(t)$. Тогда $\int_0^x {{{\psi }^{2}}(t)dt = K(x)} + R(x)$.

Далее будем считать, что ${{\varepsilon }_{1}}$ выбрано так, чтобы

Пользуясь леммой 2, рассмотрим Преобразуем сумму (см. (2)). Положим для краткости ${{\ell }_{{ij}}} = \sqrt 2 {{\gamma }_{{ij}}}$. Заметим, что Таким образом, Стало быть, Делая замену $t - \beta \ell = \tau {\text{/}}\sqrt 2 $, $\gamma = \ell {\text{/}}\sqrt 2 $, получаем Следовательно, $K(x) = c + Q(x)$, где $Q(x)$ определяется формулой (1), Таким образом, $\Psi (x) = C + Q(x) + R(x)$ $\forall x \in [ - b;b]$; $C = c + \Psi (0)$.

По построению $Q{\kern 1pt} {\kern 1pt} '{\kern 1pt} (x) = K{\kern 1pt} {\kern 1pt} '{\kern 1pt} (x) \geqslant 0$, $x \in \mathbb{R}$. Следовательно, функция $Q(x)$ не убывает. Предположим теперь, что функция $\Psi (x)$ возрастает, т.е. $\Psi {\kern 1pt} {\kern 1pt} '{\kern 1pt} (x) > 0$, $x \in [ - b;b]$. По теореме Вейерштрасса существует минимум $\omega = \mathop {\min }\limits_{x \in [ - b;b]} \Psi {\kern 1pt} '{\kern 1pt} (x)$, и он достигается в некоторой точке ${{x}_{0}} \in [ - b;b]$. Таким образом, $\omega = \Psi {\kern 1pt} '{\kern 1pt} ({{x}_{0}}) > 0$. Имеем

если ${\text{|}}r(x){\kern 1pt} {\text{|}} < \omega $, т.е. достаточно, чтобы $u({{\varepsilon }_{1}}) < \omega $. Теорема 1 доказана.

Пусть $\omega (f,\delta ) = \sup \{ {\text{|}}f(x) - f(y){\kern 1pt} {\text{|}}:{\text{|}}x - y{\kern 1pt} {\text{|}} \leqslant \delta ,\;x,y \in X\} $ – модуль непрерывности. Следующее утверждение доказано в [14, theorem 2].

Лемма 3. Для любой возрастающей функции $f \in {\mathbf{C}}(X)$, $X = [ - 1;1]$, и для всякого $n \in \mathbb{N}$ существует возрастающий многочлен ${{P}_{n}}(x)$ степени $n$ такой, что ${{\left\| {f - {{P}_{n}}} \right\|}_{{{\mathbf{C}}(X)}}} \leqslant c\omega (f,{{n}^{{ - 1}}})$.

Доказательство теоремы 2. В силу теоремы Кантора для компакта $X \subset \mathbb{R}$ и функции $f \in {\mathbf{C}}(X)$ получаем, что $\omega (f,\delta ) \to 0$ при $\delta \to + 0$. Таким образом, согласно лемме 3, непрерывную возрастающую функцию можно с любой степенью точности приблизить возрастающим многочленом (т.е. гладкой возрастающей функцией). Отсюда и непосредственно из теоремы 1 вытекает теорема 2. Теорема 2 доказана.

Доказательство теоремы 3. Непосредственно из (3) получаем

т.е. (4). Аналогично, Рассмотрим интеграл Отсюда очевидным образом получаем (5). Вычислим теперь Заметим, что Беря по частям при $u = (t - {{\beta }_{{ij}}})$, $v = - \gamma _{{ij}}^{2}{{b}_{{ij}}}(t)$, получаем Далее, Таким образом, Отсюда сразу следует (6). Теорема 3 доказана.

Доказательство теоремы 4. Поскольку функция $Q(x)$ гладкая и возрастающая, согласно теореме Вейерштрасса, найдется число $\omega > 0$ такое, что $Q{\kern 1pt} {\kern 1pt} '{\kern 1pt} (x) \geqslant \omega $ на $[ - b;b]$. В соответствии с формулой (1) для всех $x \in [ - b;b]$ имеем

Таким образом, если $2\kappa \varepsilon \leqslant \omega $, то Теорема 4 доказана.

Доказательство теоремы 5. Поскольку функция $Q(x)$ гладкая и возрастающая, согласно теореме Вейерштрасса, найдется число $\omega > 0$ такое, что $Q{\kern 1pt} '{\kern 1pt} (x) \geqslant \omega $ на $[ - b;b]$. Выберем произвольное разбиение отрезка $ - b = {{x}_{0}} < {{x}_{1}} < \ldots < {{x}_{n}} = b$, $\Delta {{x}_{i}} \geqslant \delta $, $i = \overline {1,n} $. Для произвольных $\beta \in \mathbb{R}$, $\gamma > 0$ рассмотрим точки ${{\tilde {x}}_{i}} = ({{x}_{i}} - \beta ){\text{/}}\gamma $, $i = \overline {0,n} $. Имеем

Очевидно, что $\left| r \right| \leqslant 2{{\left\| {\Phi _{0}^{'} - \Phi _{{0,\varepsilon }}^{'}} \right\|}_{{{{L}_{1}}[0;5]}}} \leqslant 2\varepsilon $. В соответствии с (1) где ${\text{|}}R{\text{|}} \leqslant \kappa {\text{|}}r{\text{|}} \leqslant 2\kappa \varepsilon $. Таким образом, при условии $2\kappa \varepsilon < \omega \delta $ имеем ${{P}_{\varepsilon }}({{x}_{i}}) > {{P}_{\varepsilon }}({{x}_{{i - 1}}})$, $i = \overline {1,n} $. Теорема 5 доказана.

3. АППРОКСИМАЦИЯ ИНТЕГРАЛА ЛАПЛАСА С ПОМОЩЬЮ КВАДРАТИЧНЫХ ЭКСПОНЕНТ

Рассмотрим четыре подхода к аппроксимации интеграла Лапласа с помощью квадратичных экспонент. Первые два из этих четырех подходов применимы, вообще говоря, к аппроксимации произвольной непрерывной функции; третий – к аппроксимации произвольной непрерывной функции с первой производной из класса ${{L}_{2}}$; четвертый – к аппроксимации произвольной функции со второй производной из класса ${{L}_{2}}$. Эти подходы и сами по себе представляют определенный интерес. Третий подход соответствует теореме 5, четвертый – теореме 4.

Применительно к интегралу Лапласа, прежде всего, заметим, что он обладает рядом замечательных свойств. Для нас важно, главным образом, то, что он, во-первых, является нечетной функцией, и, во-вторых, значения, принимаемые им при $x > 5$, практически неотличимы от числа $1{\text{/}}2$. Поэтому аппроксимацию интеграла Лапласа достаточно строить лишь на отрезке $[0;5]$. Именно этим отрезком мы и ограничимся. Подходы к аппроксимации мы ранжируем по степени возрастания сложности. Что касается иных характеристик, сделаем следующие замечания.

Наша цель состоит в том, чтобы при возможно меньшем количестве $\mu $ квадратичных экспонент получить высокоточную аппроксимацию интеграла Лапласа (вполне достаточную для практических приложений). Упомянутые четыре подхода – это четыре возможности добиться поставленной цели. Первый – наиболее простой и быстрый в смысле реализации. Однако он носит эвристический характер, поскольку сводится к минимизации квадратичной невязки по фиксированному набору контрольных точек. Таким образом, малость невязки по норме ${{L}_{p}}[a;b]$ или ${\mathbf{C}}[a;b]$ строго теоретически не гарантируется. Второй, третий и четвертый имеют теоретически обоснованную оценку точности.

Подход I. При заданном $\mu \in \mathbb{N}$ и наборе параметров $(\xi ,\eta ,\zeta ) \in {{\mathbb{R}}^{{3\mu }}}$ в качестве аппроксимирующих функций (аппроксимаций) будем использовать функции вида

где $\varepsilon > 0$ – некоторое достаточно малое число (нужное для того только, чтобы избежать нуля в знаменателе; в численных расчетах мы везде брали $\varepsilon = {{10}^{{ - 10}}}$). Пусть $a = {{z}_{1}} < \ldots < {{z}_{m}} = b$ – некоторое разбиение отрезка $[a;b] = [0;5]$ с помощью достаточно большого количества контрольных точек. Мы брали $m = 1000$. Подбор оптимальных параметров аппроксимации производится путем минимизации квадратичной невязки за счет варьирования всей тройки $(\xi ,\eta ,\zeta ) \in {{\mathbb{R}}^{{3\mu }}}$. Задачу минимизации мы решали различными методами: Левенберга–Марквардта, Гаусса–Ньютона, Хука–Дживса, DFP, BFGS, BFS. Для одномерного поиска использовались метод квадратичной интерполяции, Брента и More–Thuente (см. [21]). Наилучшие результаты показал метод Левенберга–Марквардта в сочетании с методом More–Thuente.

Для случая $\mu = 2$ были получены следующие результаты:

$\xi = (0.50102, - 0.74041)$, $\eta = (0.16707, - 1.0234)$, $\zeta = (91.92,1.6321)$, $W = 7.548{\text{e}}{\kern 1pt} - {\kern 1pt} 5$ (запись ${\text{e}}{\kern 1pt} - {\kern 1pt} 5$ означает умножение на ${{10}^{{ - 5}}}$), максимальное отклонение в контрольных точках аппроксимации от аппроксимируемой функции ${{\delta }_{{{\text{max}}}}} = 0.0013038$, максимальное отклонение в контрольных точках производной аппроксимации от производной аппроксимируемой функции $\delta _{{{\text{max}}}}^{'} = 0.014958$. Графики аппроксимации и аппроксимируемой функции не приводим, поскольку различия в них можно увидеть только при очень большом увеличении (на малом участке между точками 3 и 4). Таким образом, достаточно хорошее качество аппроксимации удается получить уже этим простейшим способом при $\mu = 2$. Для сравнения укажем, что при использовании классического способа (взять частичную сумму ряда, полученного почленным интегрированием разложения подынтегральной функции в ряд Маклорена) приходится учитывать порядка сорока членов разложения для достижения сопоставимой точности.

Приведем результаты для случая $\mu = 3$:

$\xi \, = \,(0.89868, - 0.80742, - 0.012592)$, $\eta \, = \,( - 110.5396, - 1.1287,2.6466)$, $\zeta \, = \,(151.7562,1.7226,2.094)$, $W = 1.0964{\text{e}}{\kern 1pt} - {\kern 1pt} 5$, ${{\delta }_{{{\text{max}}}}} = 0.00052134$, $\delta _{{{\text{max}}}}^{'} = 0.0072675$. Графики не приводим, поскольку различия в них можно увидеть только при очень большом увеличении.

Подход II. Этот подход основан на следующем общем утверждении об аппроксимации.

Теорема 6. Пусть $E$ – линейное пространство со скалярным произведением $\langle .,.\rangle $, ${{x}_{i}} \in E$, $i = \overline {1,m} $, ${{\zeta }_{{ij}}} = \langle {{x}_{i}},{{x}_{j}}\rangle $, $i,j = \overline {1,m} $, $\Gamma = ({{\zeta }_{{ij}}}{{)}_{{i,j = \overline {1,m} }}}$ – матрица Грама указанной системы векторов, ${{\Gamma }_{1}} = ({{\zeta }_{{ij}}}{{)}_{{i,j = \overline {1,m - 1} }}}$, $\Delta = {\text{det}}\Gamma $, ${{\Delta }_{1}} = {\text{det}}{{\Gamma }_{1}}$. Предположим, что ${{\Delta }_{1}} \ne 0$ (отсюда по свойствам определителя Грама вытекает, что ${{\Delta }_{1}} > 0$). Тогда найдутся вектор $\widetilde \xi \in {{\mathbb{R}}^{{m - 1}}}$ и соответствующая ему линейная комбинация $\tilde {x} = \sum\nolimits_{i = 1}^{m - 1} {{{{\widetilde \xi }}_{i}}{{x}_{i}}} $ такие, что для разности $h = \tilde {x} - {{x}_{m}}$ имеем ${{\left\| h \right\|}^{2}} = \Delta {\text{/}}{{\Delta }_{1}}$. Более того, вектор $\tilde {x}$ минимизирует расстояние от вектора ${{x}_{m}}$ до линейного подпространства $L$, натянутого на векторы ${{x}_{i}}$, $i = \overline {1,m - 1} $, причем вектор $\widetilde \xi $ определяется однозначно как решение системы

Доказательство. Согласно [22, ${{\S}}$ XVIII.1, пример 9, с. 601], существует единственный вектор $\widetilde \xi $ (определяемый как решение системы (7)) такой, что

По построению вектор $h = \tilde {x} - {{x}_{m}}$ ортогонален пространству $L$, причем $\tilde {x} \in L$. Согласно второму следствию из [23, гл. 12, теорема 98.2, с. 349–350], справедливо равенство где $\left| {\Gamma ({{x}_{1}}, \ldots ,{{x}_{{m - 1}}},h)} \right|$ – определитель Грама системы ${{x}_{1}}, \ldots ,{{x}_{{m - 1}}},h$, и соответственно, указанное равенство переписывается как Известно (см. [23, гл. 12, свойства 2, 3, с. 348]), что определитель Грама не изменяется от прибавления к любому вектору системы линейной комбинации остальных векторов, и при умножении любого вектора системы на число $\alpha $ умножается на ${{\alpha }^{2}}$. С учетом принадлежности $\tilde {x} \in L$ отсюда ясно, что $\left| {\Gamma ({{x}_{1}}, \ldots ,{{x}_{{m - 1}}},h)} \right| = \Delta $. Теорема доказана.

Замечание 1. В частности, при $\Delta = 0$ получаем в качестве следствия классический факт: из равенства нулю определителя Грама вытекает линейная зависимость системы векторов. Вообще, как известно, равенство нулю определителя Грама необходимо и достаточно для линейной зависимости системы векторов.

Замечание 2. Для аппроксимации ${{x}_{m}} \approx \sum\nolimits_{i = 1}^{m - 1} {{{{\widetilde \xi }}_{i}}{{x}_{i}}} $ точность оценивается как $\sqrt {\Delta {\text{/}}{{\Delta }_{1}}} $.

Обозначим для краткости ${{\varphi }_{j}}(x) = {{\varphi }_{j}}[{{\eta }_{j}},{{\zeta }_{j}}](x)$, $j = \overline {1,\mu } $. Рассмотрим систему функций $\{ {{\varphi }_{j}},j = \overline {1,\mu } ;{{\widehat \Phi }_{0}}\} $ как систему из $(\mu + 1)$ векторов в пространстве ${{L}_{2}}[a;b]$. Положим $\Gamma $ – матрица Грама этой системы, ${{\Gamma }_{1}}$ – матрица Грама системы $\{ {{\varphi }_{j}},j = \overline {1,\mu } \} $, $\Delta = \det \Gamma $, ${{\Delta }_{1}} = \det {{\Gamma }_{1}}$. Непосредственно из теоремы 6 вытекает

Теорема 7. Пусть ${{\Delta }_{1}} \ne 0$ (следовательно, ${{\Delta }_{1}} > 0$). Тогда существует вектор $\xi \in {{\mathbb{R}}^{\mu }}$ такой, что для аппроксимации ${{\Phi }_{\mu }}(x) = \sum\nolimits_{j = 1}^\mu {{{\xi }_{j}}{{\varphi }_{j}}(x)} $ имеем ${{\left\| {{{\Phi }_{\mu }} - {{{\widehat \Phi }}_{0}}} \right\|}_{{{{L}_{2}}[a;b]}}} = \sqrt {\Delta {\text{/}}{{\Delta }_{1}}} $. При этом в качестве $\xi $ можно взять решение системы

где

Из [24, гл. 4, упражнение 15, с. 92] видно, что в теореме 7 вектор $\xi $ как раз и обеспечивает минимум расстояния аппроксимируемой функции ${{\widehat \Phi }_{0}}$ до подпространства, натянутого на систему $\{ {{\varphi }_{j}},j = \overline {1,\mu } \} $.

Исходя из теоремы 7, при заданном $\mu \in \mathbb{N}$ возникает идея подобрать параметры $\eta ,\zeta \in {{\mathbb{R}}^{\mu }}$ путем минимизации отношения $R(\eta ,\zeta ) = \Delta (\eta ,\zeta ){\text{/}}{{\Delta }_{1}}(\eta ,\zeta )$, скажем, методом Хука–Дживса (нулевого порядка), а затем взять в качестве вектора коэффициентов $\xi \in {{\mathbb{R}}^{\mu }}$ решение системы (8). Впрочем, здесь следует соблюдать осторожность: при неудачном выборе начального приближения минимизируемая величина может повести себя как неопределенность типа $\left( {\frac{0}{0}} \right)$. С другой стороны, эту трудность можно обойти, если ввести штраф за слишком малые значения ${{\Delta }_{1}}$ (но, как показывает опыт, в таком случае происходит сильное замедление процесса численной минимизации). Отметим, наконец, что отношение $R(\eta ,\zeta )$ нельзя вычислить точно (аналитически). Поэтому приходится его вычислять приближенно с помощью численного интегрирования. Тем не менее близость к оптимуму полученных в итоге параметров $\xi ,\eta ,\zeta $ можно оценить с помощью квадратичной невязки $W[\xi ,\eta ,\zeta ]$, а также оценкой максимального отклонения аппроксимации от аппроксимируемой функции в контрольных точках.

Описанный способ позволил нам получить следующие результаты.

Для $\mu = 2$: $\xi = ( - 0.74267,1.6952)$, $\eta = ( - 1.0268, - 2798.7635)$, $\zeta = (1.6348,2536.7604)$, $W = 6.6945{\text{e}}{\kern 1pt} - {\kern 1pt} 6$, ${{\delta }_{{{\text{max}}}}} = 0.0012945$, $\delta _{{{\text{max}}}}^{'} = 0.014745$. Графики не приводим, поскольку различия в них можно увидеть только при очень большом увеличении.

Для $\mu \, = \,3$: $\xi \, = \,( - 0.32242, - 0.74678,0.82234)$, $\eta \, = \,(4.2016, - 1.05,4.1338)$, $\zeta = (8.2595,1.65,13.411)$, $W = 7.4844{\text{e}}{\kern 1pt} - {\kern 1pt} 5$, ${{\delta }_{{{\text{max}}}}} = 0.00079101$, $\delta _{{{\text{max}}}}^{'} = 0.011011$. Графики не приводим, поскольку различия в них можно увидеть только при очень большом увеличении.

Подход III аналогичен подходу II, с тем лишь отличием, что вместо аппроксимации самого интеграла Лапласа ${{\widehat \Phi }_{0}}(x)$ точно тем же способом строится аппроксимация (в пространстве ${{L}_{2}}[a;b]$) по функциям $\varphi _{i}^{'}(x)$ его производной $\widehat \Phi _{0}^{'}(x) = 1{\text{/}}\sqrt {2\pi } \exp \left[ { - \frac{{{{x}^{2}}}}{2}} \right]$. После этого в качестве аппроксимации функции ${{\widehat \Phi }_{0}}(x)$ принимается ${{Q}_{\mu }}[\xi ,\eta ,\zeta ](x) = C + {{\Phi }_{\mu }}[\xi ,\eta ,\zeta ](x)$, где константа $C$ определяется из условия ${{Q}_{\mu }}(a) = {{\widehat \Phi }_{0}}(a)$.

Итак, рассмотрим систему функций

как систему из $(\mu + 1)$ векторов в пространстве ${{L}_{2}}[a;b]$. Положим $\Gamma {\kern 1pt} '$ – матрица Грама этой системы, $\Gamma _{1}^{'}$ – матрица Грама системы $\Delta {\kern 1pt} ' = \det \Gamma {\kern 1pt} '$, $\Delta _{1}^{'} = \det \Gamma _{1}^{'}$. Непосредственно из теоремы 6 вытекает

Теорема 8. Пусть $\Delta _{1}^{'} \ne 0$ (следовательно, $\Delta _{1}^{'} > 0$). Тогда существует вектор $\xi \in {{\mathbb{R}}^{\mu }}$ такой, что для аппроксимации ${{Q}_{\mu }}(x) = C + {{\Phi }_{\mu }}(x)$, $C = {{\widehat \Phi }_{0}}(a) - {{\Phi }_{\mu }}(a)$, имеем ${{\left\| {{{Q}_{\mu }} - {{{\widehat \Phi }}_{0}}} \right\|}_{{{\mathbf{C}}[a;b]}}} \leqslant \sqrt {(b - a)\Delta {\kern 1pt} '{\text{/}}\Delta _{1}^{'}} $. При этом в качестве $\xi $ можно взять решение системы

где

Доказательство. Используя неравенство Гёльдера, получаем

Дальнейшее очевидно. Теорема доказана.

Исходя из теоремы 8, при заданном $\mu \in \mathbb{N}$ возникает идея подобрать параметры $\eta ,\zeta \in {{\mathbb{R}}^{\mu }}$ путем минимизации отношения $R{\kern 1pt} '(\eta ,\zeta ) = \Delta {\kern 1pt} '(\eta ,\zeta ){\text{/}}\Delta _{1}^{'}(\eta ,\zeta )$, скажем, методом Хука–Дживса (нулевого порядка), а затем взять в качестве вектора коэффициентов $\xi \in {{\mathbb{R}}^{\mu }}$ решение системы (9). Отношение $R{\kern 1pt} '(\eta ,\zeta )$ опять же нельзя вычислить точно (аналитически). Поэтому приходится его вычислять приближенно с помощью численного интегрирования. Тем не менее близость к оптимуму полученных в итоге параметров $\xi ,\eta ,\zeta $ можно оценить с помощью квадратичной невязки $W[\xi ,\eta ,\zeta ]$, а также оценкой максимального отклонения аппроксимации от аппроксимируемой функции в контрольных точках.

Описанный способ позволил нам получить следующие результаты.

Для $\mu = 2$: $C = 0.50001$, $\xi = ( - 0.30633, - 0.7523)$, $\eta = ( - 0.43757, - 1.9128)$, $\zeta = (1.453,1.7258)$, квадратичное отклонение аппроксимации от аппроксимируемой функции в контрольных точках $\widetilde W = 6.3837{\text{e}}{\kern 1pt} - {\kern 1pt} 8$, ${{\delta }_{{{\text{max}}}}} = 1.4717{\text{e}}{\kern 1pt} - {\kern 1pt} 5$, $\delta _{{{\text{max}}}}^{'} = 7.9605{\text{e}}{\kern 1pt} - {\kern 1pt} 5$. Рисунков не приводим, поскольку визуально график аппроксимации полностью совпадает с графиком аппроксимируемой функции. Для практических приложений такой точности уже вполне достаточно.

Для $\mu = 3$: $C = 0.5$, $\xi = ( - 0.025654, - 0.27286, - 0.76869)$, $\eta = (0.47302, - 0.40261, - 1.824)$, $\zeta = (1.3121,1.3683,1.6522)$, $\widetilde W = 3.5108{\text{e}}{\kern 1pt} - {\kern 1pt} 9$, ${{\delta }_{{{\text{max}}}}} = 3.2107{\text{e}}{\kern 1pt} - {\kern 1pt} 6$, $\delta _{{{\text{max}}}}^{'} = 4.8288{\text{e}}{\kern 1pt} - {\kern 1pt} 5$. Рисунков не приводим, поскольку визуально график аппроксимации полностью совпадает с графиком аппроксимируемой функции.

Подход IV совершенно аналогичен подходу III, с тем лишь отличием, что вместо аппроксимации производной интеграла Лапласа $\widehat \Phi _{0}^{'}(x)$ точно тем же способом строится аппроксимация (в пространстве ${{L}_{2}}[a;b]$) по функциям $\varphi _{i}^{{''}}(x)$ его второй производной $\widehat \Phi _{0}^{{''}}(x)$. После этого в качестве аппроксимации ${{\widehat \Phi }_{0}}(x)$ принимается ${{Q}_{\mu }}[\xi ,\eta ,\zeta ](x) = {{C}_{1}} + {{C}_{2}}x + {{\Phi }_{\mu }}[\xi ,\eta ,\zeta ](x)$, где ${{C}_{1}},{{C}_{2}}$ определяются из условий ${{Q}_{\mu }}(0) = {{\widehat \Phi }_{0}}(0)$, $Q_{\mu }^{'}(0) = \widehat \Phi _{0}^{'}(0)$.

Описанный способ позволил нам получить следующие результаты.

Для $\mu = 2$: ${{C}_{1}} = 0.50023$, ${{C}_{2}} = - 4.8295{\text{e}}{\kern 1pt} - {\kern 1pt} 005$, $\xi = ( - 0.50169, - 0.60577)$, $\eta = ( - 0.64487, - 2.0532)$, $\zeta = (1.5027,1.456)$, по функции ($\widetilde W$ – квадратичное отклонение, ${{\delta }_{{{\text{max}}}}}$ – максимальное отклонение) $\widetilde W = 6.125{\text{e}}{\kern 1pt} - {\kern 1pt} 7$, ${{\delta }_{{{\text{max}}}}} = 4.9302{\text{e}}{\kern 1pt} - {\kern 1pt} 5$; по первой производной $\widetilde W{\kern 1pt} ' = 1.4818{\text{e}}{\kern 1pt} - {\kern 1pt} 6$, $\delta _{{{\text{max}}}}^{'} = 7.5314{\text{e}}{\kern 1pt} - {\kern 1pt} 5$; по второй производной $\widetilde W{\kern 1pt} '{\kern 1pt} ' = 9.9166{\text{e}}{\kern 1pt} - {\kern 1pt} 6$, $\delta _{{{\text{max}}}}^{{''}} = 0.00056108$. Рисунков не приводим, поскольку визуально график аппроксимации полностью совпадает с графиком аппроксимируемой функции. Как видим, разница с подходом III не очень большая, из чего можно заключить, что и подход III будет достаточно работоспособен. То, что по функции отклонение получилось даже несколько больше, чем для подхода III, не должно нас смущать, поскольку отклонение вычислялось не от “идеального” интеграла Лапласа, а от его приближения, полученного численным интегрированием. Поэтому отклонение по производной здесь более информативно.

Для $\mu \, = \,3$ результаты следующие: ${{C}_{1}}\, = \,0.49988$, ${{C}_{2}}\, = \,2.5236{\text{e}}{\kern 1pt} - {\kern 1pt} 5$, $\xi \, = \,(10.9656, - 0.28145, - 11.7683)$, $\eta = ( - 1.9261, - 0.46093, - 1.9349)$, $\zeta = (1.9636,1.4548,1.9586)$, по функции $\widetilde W = 7.5652{\text{e}}{\kern 1pt} - {\kern 1pt} 9$, ${{\delta }_{{{\text{max}}}}} = 9.2179{\text{e}}{\kern 1pt} - {\kern 1pt} 6$; по первой производной $\widetilde W{\kern 1pt} ' = 8.9872{\text{e}}{\kern 1pt} - {\kern 1pt} 8$, $\delta _{{{\text{max}}}}^{'} = 2.202{\text{e}}{\kern 1pt} - {\kern 1pt} 5$; по второй производной $\widetilde W{\kern 1pt} '{\kern 1pt} ' = 1.222{\text{e}}{\kern 1pt} - {\kern 1pt} 6$, $\delta _{{{\text{max}}}}^{{''}} = 0.00021175$. Рисунков не приводим, поскольку визуально график аппроксимации полностью совпадает с графиком аппроксимируемой функции.

В заключение этого раздела приведем также следующие данные. Нам было интересно сравнить время счета при использовании нашего способа аппроксимации интеграла Лапласа с помощью линейной функции и двух квадратичных экспонент (подход IV), с одной стороны, и при использовании функции ошибки erf, реализованной в системе MATLAB (седьмой версии), с другой стороны. Оценивалось совокупное время, необходимое для вычисления интеграла Лапласа (точнее, его аппроксимации) на массиве узловых точек сетки отрезка $[0;5]$ с шагом ${{10}^{{ - 6}}}$ (т.е. на массиве из порядка пяти миллионов точек). Разумеется, время счета определенным образом зависит от других процессов, функционирующих в системе (а также от типа и конфигурации компьютера и операционной системы). Поэтому ради чистоты эксперимента мы повторили его 20 раз. Во всех 20 из них вычисления по нашему способу оказались несколько быстрее. Приведем типичный результат. Для нашего способа время счета $0.317164$ с; при использовании функции erf – $0.384069$ с. Максимум отклонений в узлах сетки: $4.4487 \times {{10}^{{ - 5}}}$. Графики визуально полностью совпадают. При этом следует понимать, что MATLAB – это коммерческий продукт и алгоритм, реализованный в программе erf – это коммерческая тайна. Можно предположить, что там реализован табличный способ в сочетании с линейной интерполяцией, но автору об этом ничего не известно. Остается лишь отметить, что подобный способ требует хранения в памяти компьютера соответствующего массива данных, а при собственной реализации такой массив еще надо создать.

4. РЕЗУЛЬТАТЫ ЧИСЛЕННЫХ ЭКСПЕРИМЕНТОВ

Для неубывающих функций $f(t)$, $t \in [a;b]$, различного вида мы строили аппроксимацию по формуле (3), $C = f(a)$, при значениях $\nu \geqslant 2$ – до достижения достаточно высокой точности. Параметры $(\alpha ,\beta ,\gamma ) \in {{\mathbb{R}}^{{3\nu }}}$ определялись, исходя из минимизации квадратичной невязки по тысяче контрольных точек. Вместо интеграла Лапласа использовалась его аппроксимация в рамках подхода IV при $\mu = 2$. При использовании подхода III и встроенной функции MATLAB ${\text{erf}}(t)$ для вычисления интеграла Лапласа отличия практически неразличимы. Например, для теста № 1 при $\nu = 2$ для подхода III получаем $\delta = 0.00085953$, ${{\delta }_{{\max }}} = 0.0026705$; $\alpha = (0.36985,0.62644)$, $\beta = (1.6809{\text{e}}{\kern 1pt} - {\kern 1pt} 5,\;8.4071{\text{e}}{\kern 1pt} - {\kern 1pt} 11)$, $\gamma = (0.97662,3.437)$. При использовании функции ${\text{erf}}(t)$ $\delta = 0.0008385$, ${{\delta }_{{\max }}} = 0.0023415$; $\alpha = (0.37384,0.62218)$, $\beta = ( - 9.5407{\text{e}}{\kern 1pt} - {\kern 1pt} 12,\; - {\kern 1pt} 3.7628{\text{e}}{\kern 1pt} - {\kern 1pt} 12)$, $\gamma = (0.98342,3.4653)$. Здесь (и далее) используются следующие обозначения: $\delta $ – значение квадратичной невязки, ${{\delta }_{{\max }}}$ – модуль максимального отклонения аппроксимации от аппроксимируемой функции по набору контрольных точек. Графики по сравнению с подходом IV визуально никак не различаются (см. далее соответствующиe результаты для подхода IV). Кроме того, проводилась проверка полученной аппроксимации на монотонность. При использовании функции ${\text{erf}}(t)$ аппроксимация во всех тестах оказалась неубывающей (разумеется, на конечной сетке – из $1000$ контрольных точек). При использовании подхода III наблюдалось незначительное отклонение от монотонности в одной-двух контрольных точках: тест № 3, $\nu = 2$ – порядка ${{10}^{{ - 3}}}$; тест № 4, $\nu = 2$, $\nu = 3$ – порядка ${{10}^{{ - 3}}}$; тест № 5, $\nu = 5$ – порядка ${{10}^{{ - 4}}}$. Отметим, что нарушение положительности (отрицательности) на величину порядка ${{10}^{{ - 5}}}$ и меньше может возникать вследствие накопленной погрешности вычислений (в том числе и достаточно простых) даже там, где этого теоретически не должно быть. Автор сталкивался с подобным явлением, например, при программной реализации симплекс-метода. Итак, приведем полученные результаты для подхода IV. Здесь $\omega = \mathop {\min }\limits_{i = \overline {1,n} } \{ Q({{x}_{i}}) - Q({{x}_{{i - 1}}})\} $ – сеточный показатель монотонности.

Тест № 1. Функция $f(t) = \operatorname{arctg} t$, $t \in [ - 3;3]$. При $\nu = 2$: $\delta = 0.00085827$, ${{\delta }_{{\max }}} = 0.0027058$, $\omega = 0.00051423$; $\alpha = (0.36984,0.62644)$, $\beta = (2.4837{\text{e}}{\kern 1pt} - {\kern 1pt} 5,\;2.0265{\text{e}}{\kern 1pt} - {\kern 1pt} 10)$, $\gamma = (0.97661,3.437)$. При $\nu = 3$: $\delta = 3.1132{\text{e}}{\kern 1pt} - {\kern 1pt} 5$, ${{\delta }_{{\max }}} = 0.00050259$, $\omega = 0.00062426$; $\alpha = (0.20861,0.46246,0.328)$, $\beta = (3.2261{\text{e}}{\kern 1pt} - {\kern 1pt} 5,\;1.2145{\text{e}}{\kern 1pt} - {\kern 1pt} 5,\;1.7191{\text{e}}{\kern 1pt} - {\kern 1pt} 11)$, $\gamma = (0.79389,1.7984,9.0161)$. Графики не приводим, поскольку визуально они полностью совпадают.

Тест № 2. Функция $f(t) = \sin t$, $t \in [ - \pi {\text{/}}2;\pi {\text{/}}2]$. При $\nu = 2$: $\delta = 0.015425$, ${{\delta }_{{\max }}} = 0.013986$, $\omega = 0.00036968$, $\alpha = (5.4264, - 4.4197)$, $\beta = ( - 0.023435, - 0.14461)$, $\gamma = (3.9465,9.0997)$. Графики см. на фиг. 1: аппроксимация изображена пунктирной линией, аппроксимируемая функция – сплошной линией. При $\nu = 3$: $\delta = 0.0010741$, ${{\delta }_{{\max }}} = 0.0049779$, $\omega = 0.00025331$; $\alpha = ( - 0.38012, - 0.941, - 0.38081)$, $\beta = (1.0384,0.00052416, - 1.038)$, $\gamma = ( - 0.6554,1.0154,0.65592)$. Графики не приводим, поскольку визуально они полностью совпадают.

Тест № 3. Функция $f(t) = \cos t$, $t \in [ - \pi ;0]$. При $\nu = 2$: $\delta = 0.027529$, ${{\delta }_{{\max }}} = 0.015041$, $\omega = - 0.0011892$: нарушение монотонности лишь в одной контрольной точке с индексом $905$; $\alpha = ( - 4.5679,5.4914)$, $\beta = ( - 0.29867, - 1.3424)$, $\gamma = (8.3551,3.8693)$. При использовании для аппроксимации интеграла Лапласа подхода IV при $\mu = 3$ нарушение монотонности устраняется: $\delta = 0.028168$, ${{\delta }_{{\max }}} = 0.015552$, $\omega = 8.1668{\text{e}}{\kern 1pt} - {\kern 1pt} 5$. Графики см. на фиг. 2. При $\nu = 3$: $\delta = 0.0020035$, ${{\delta }_{{\max }}} = 0.0085406$, $\omega = 4.0173{\text{e}}{\kern 1pt} - {\kern 1pt} 5$; $\alpha = (0.28284, - 1.6716,0.95752)$, $\beta = ( - 2.4766,0.57694, - 1.3151)$, $\gamma = (0.76281,0.44424,1.4416)$. Графики не приводим, поскольку визуально они полностью совпадают.

Тест № 4. Функция $f(t) = \exp t$, $t \in [ - 3;3]$. При $\nu = 2$: $\delta = 0.010551$, ${{\delta }_{{\max }}} = 0.015452$, $\omega = 2.173{\text{e}}{\kern 1pt} - {\kern 1pt} 4$; $\alpha = (13.0709,2.7665)$, $\beta = (7.5862,5.8611)$, $\gamma = (3.5002,5.5059)$. При $\nu = 3$: $\delta = 0.0051879$, ${{\delta }_{{\max }}} = 0.013392$, $\omega = 0.00030832$; $\alpha = (6.5988,4.0538,0.13998)$, β = $\beta = (6.1664,5.9167,0.72595)$, $\gamma = (2.6424,4.7098,8.6292)$. Графики не приводим, поскольку визуально они полностью совпадают.

Тест № 5. Непрерывная кусочно-линейная функция, см. фиг. 3–4, при $t \in [0;3]$. При $\nu = 3$: $\delta = 1.376$, ${{\delta }_{{\max }}} = 0.12811$, $\omega = 0.00012769$; $\alpha = (13.8495,3.6601, - 0.22289)$, β = $\beta = (11.0969,5.8107,2.2459)$, $\gamma = (3.8892,3.9442,9.1205)$. Графики см. на фиг. 3. При $\nu = 4$: $\delta = 0.21604$, ${{\delta }_{{\max }}} = 0.068932$, $\omega = 0.00011881$; $\alpha = ( - 6.7454,5.2712,2.4614, - 1.9459)$, β = $\beta = (1.0425,0.89296,1.3388,2.7021)$, $\gamma = (0.61874,0.49289,0.3392,0.66608)$. Графики см. на фиг. 4. При $\nu = 5$: $\delta = 0.049677$, ${{\delta }_{{\max }}} = 0.030583$, $\omega = - 0.00024438$; – нарушение монотонности лишь в одной контрольной точке с индексом $256$; $\alpha = ( - 0.33219, - 7.8355,6.5769, - 0.49263, - 0.58114)$, $\beta = (0.22402,2.3705,2.337,2.4544,0.9735)$, $\gamma = (0.54247,0.70341,0.58191,0.13697, - 0.27281)$. При использовании для аппроксимации интеграла Лапласа подхода IV при $\mu = 3$ нарушение монотонности существенно уменьшается (остается в одной контрольной точке с индексом $204$): $\delta = 0.049884$, ${{\delta }_{{\max }}} = 0.030384$, $\omega = - 8.6895{\text{e}}{\kern 1pt} - {\kern 1pt} 5$. При использовании для аппроксимации интеграла Лапласа функции ${\text{erf}}$ нарушение монотонности полностью устраняется: $\delta = 0.049594$, ${{\delta }_{{\max }}} = 0.030323$, $\omega = 0.00023876$. Графики не приводим, поскольку визуально они практически совпадают. При $\nu = 6$: $\delta = 0.027565$, ${{\delta }_{{\max }}} = 0.031594$, $\omega = 2.3785{\text{e}}{\kern 1pt} - {\kern 1pt} 4$; α = $\alpha = - (0.46309,0.43636,1.8492,1.1983,0.33155, - 0.37226)$, β = (0.92851, 2.462, 2.8463, 1.7451, 0.21825, $1.4669)$, $\gamma = (0.22713,0.12356, - 0.56682,0.62113,0.53356, - 0.086174)$.

Графики не приводим, поскольку визуально они полностью совпадают.

Сделаем некоторые комментарии. Функции тестов № 1–4 гладкие, поэтому подпадают под ситуацию, описанную в теореме 1. Функция теста № 5 негладкая, но непрерывная и строго возрастающая, следовательно, подпадает под ситуацию, описанную в теореме 2. В численных экспериментах этого раздела для минимизации невязки мы использовали численные методы первого порядка (в основном, метод Левенберга–Марквардта и метод DFP – в зависимости от того, какой из них показывал лучшие результаты, с одномерным поиском по методу More–Thuente). Поэтому требовалось вычисление производных невязки по параметрам аппроксимации, и тут использовались формулы из теоремы 3.