Методы интерполяции и сглаживания сплайнами

Постановка задачи и основные положения

Как отмечалось ранее, в вычислительной практике при глобальном способе аппроксимации высокие степени интерполяционных многочленов использовать нецелесообразно. Поэтому часто применяется кусочно-глобальный способ на основе, например, линейной и квадратичной (параболической) интерполяции. Однако производные таких интерполяционных многочленов, получающихся на частичных отрезках, в местах их стыка терпят разрывы. На рис. 4.10,а показано поведение совокупности трех таких параболических многочленов (звеньев). Легко видеть, что производные в точках стыка частичных отрезков (шаблонов [math](x_{i},x_{i+1},x_{i+2}),[/math] [math](x_{i+2}, x_{i+3},x_{i+4}),[/math] [math](x_{i+4},x_{i+5},x_{i+6})[/math]) терпят разрывы. Это является крайне нежелательным, если вся совокупность [math]\textstyle{\bigcup_{i}S_{2,i}(x)}[/math] используется, например, для аппроксимации обводов конструктивных элементов сложных поверхностей. Для изготовления таких поверхностей требуется задание в узлах не только значений [math]S_{k,i}(x_i)[/math], но и первой производной [math]S'_{k,i}(x_i)[/math] или даже второй производной [math]S''_{k,i}(x_i)[/math]. Таким образом, вместе с функцией [math]y_i= f(x_i),~ i=\overline{0,n}[/math], требуется достаточно хорошо аппроксимировать и ее производные. Если восполняемая функция достаточно гладкая, то и ее производные тоже гладкие, и поэтому они должны быть непрерывными во всех внутренних узлах. В связи с этим к аппроксимирующим функциям [math]S_m(x)[/math] ([math]m[/math]-степень многочлена) предъявляются требования непрерывности [math]S_m(x),\, S'_m(x)[/math] и по возможности [math]S''_m(x)[/math] или [math]S'''_m(x)[/math] во всех точках [math]x\in[a,b][/math] (рис. 4.10,б). Это и обусловливает необходимость построения сплайн-функций, обладающих указанными свойствами и имеющих интерполяционный или сглаживающий характер.

Пусть на отрезке [math][a,b][/math] задана сеточная функция [math]y_i= f(x_i)[/math] на сетке [math]\Omega_n[/math]. Требуется восполнить ее функцией [math]S_m(x)[/math], где [math]m[/math] — степень многочлена, кусочно-глобальным способом. Дадим сначала упрощенное определение сплайна, которое затем уточним.

Сплайн-функцией или сплайном называется совокупность [math]S_{m,i}(x)[/math] — алгебраических многочленов степени [math]m[/math] (звеньев), определенных на частичных отрезках [math][x_{i},x_{i+1}],~ i=\overline{0,n-1}[/math], и соединенных вместе по всем частичным отрезкам так, чтобы можно было составить многозвенную функцию [math]\textstyle{S_m(x)= \bigcup\limits_{i=0}^{n-1} S_{m,i}(x)}[/math], определенную и непрерывную на всем отрезке [math][a,b][/math] вместе со всеми своими производными [math]S_{m}^{(p)}(x)[/math] до некоторого их порядка [math]p=1,2,\ldots[/math]. Разность между [math]m[/math] и наибольшим порядком производной, непрерывной на отрезке [math][a,b][/math], определяет дефект сплайна [math]q[/math]. Условиями согласования звеньев [math]S_{m,i}(x)[/math] сплайна с исходной функцией [math]y_i= f(x_i)[/math] на соответствующем частичном отрезке [math][x_{i}, x_{i+1}][/math] называются условия, накладываемые на невязки дифференциального и интегрального типов [math]\delta S_{m,i}^{(p)}(x_i),\, \delta S_{m,i}(I_{i}^{i+1})[/math] и использующиеся для вывода формулы одного звена сплайна на указанном частичном отрезке.

Сплайн, удовлетворяющий условиям (4.3) нулевых функциональных невязок, называется интерполяционным, а удовлетворяющий только интегральным условиям (4.5), — сглаживающим или интегрально-сглаживающим.

Сплайн, построение которого основано только на дифференциальных условиях согласования (4.2) при [math]p=0,\!1[/math] или [math]p=0,\!2[/math], называется дифференциальным (Д-сплайном). Если же наряду с дифференциальными условиями согласования используется интегральное условие вида (4.5), то сплайн называется интегрально-дифференциальным (ИД-сплайном).

Количество условий согласования, необходимых для получения формулы одного звена сплайна, должно соответствовать степени сплайна (число условий на единицу больше [math]m[/math]). После определения формулы [math]S_{m,i}(x)[/math] ее правая часть выражается через известные и неизвестные величины (параметры): [math]f^{(p)}(x_i)[/math] — для дифференциального сплайна или через совокупность [math]f^{(p)}(x_i)[/math] и [math]\textstyle{I_{i}^{i+1}= \int\limits_{x_i}^{x_{i+1}} f(x)\,dx}[/math] —для интегрально-дифференциального. Эти величины будут называться параметрами сплайна. В зависимости от того, заданы те или иные параметры сплайна в постановке задачи или нет, они именуются определенными или неопределенными. Последние вычисляются на основе условий непрерывности (гладкости) сплайна [math]S_m(x)[/math], которые здесь называются условиями стыковки. Для некоторого узла [math]x_i[/math], общего для звена [math]i-1[/math], относящегося к частичному отрезку [math][x_{i-1},x_i][/math], и звена [math]i[/math], определенного на частичном отрезке [math][x_i,x_{i+1}][/math], это условие имеет вид

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

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

При решении задачи сплайн-аппроксимации необходимо предварительно вычислить неопределенные параметры сплайна для всех его звеньев. Этот этап соответствует решению параметрической задачи. Для пояснения ее характера введем понятие порядка параметров сплайн-функции, который обозначим через [math]l[/math]. Этот порядок зависит от порядка производных исходной функции [math]y=f(x)\colon[/math] для параметра [math]I_i^{i+1}[/math] примем [math]l=-1[/math], а для [math]f^{(p)}(x)[/math] положим [math]l=p~ (p=0,1,2,\ldots)[/math].

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

Две совокупности параметров, относящиеся к параметрическим уравнениям и сплайнам разных степеней, называются подобными, если они могут быть получены друг из друга путем соответствующего изменения порядков параметров [math]l[/math] на одинаковую величину.

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

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

Сплайн называется локальным, если все неопределенные параметры, относящиеся к каждому его звену [math]S_{m,t}(x),~ t=i[/math], при [math]x\in [x_i, x_{i+1}],~ i=\overline{0,n-1}[/math], находятся локально, т.е. независимо от параметров, характеризующих все остальные (или почти все остальные) звенья [math]S_{m,t}(x),~ (t=0,1,\ldots,i-1,i+1,\ldots,n-1)[/math].

Сплайн-аппроксимация (как правило, это интерполяция) локальными сплайнами сводится, таким образом, к получению конкретных значений коэффициентов многочленов [math]S_{m,i}(x)[/math] для каждого звена результирующего сплайна [math]S_{m}(x)[/math] и неопределенных параметров путем их вычисления по аппроксимационным формулам. Недостаток локальной интерполяции состоит в том, что таким способом не удается обеспечить минимальный дефект сплайна, т.е. максимально возможную его гладкость в смысле непрерывности производных как можно большего порядка [math](p=0,1,2,\ldots)[/math].

Альтернативой локальному сплайну является глобальный сплайн, для которого неопределенные параметры, относящиеся к каждому его звену [math]S_{m,t}(x),~ y=i[/math] при [math]x\in [x_i, x_{i+1}],~ i=\overline{0,n-1}[/math], находятся совместно с параметрами, характеризующими все остальные звенья [math]S_{m,t}(x),~ t\ne i[/math]. Неопределенные параметры в глобальных сплайнах для всех звеньев вычисляются, как правило, путем решения системы линейных алгебраических уравнений трех-диагонального вида методом прогонки, причем это обеспечивает непрерывность одной из производных, а именно той, которая не включена в условия согласования. Если в условие согласования включена производная другого порядка, то ее непрерывность также гарантируется. Глобальный способ аппроксимации по сравнению с локальным способом обеспечивает минимально возможный дефект сплайна. Поэтому глобальные сплайны широко используются в вычислительной практике.

В заключение данного раздела дадим развернутое определение глобального интегрально-дифференциального сплайна.

Функция [math]\textstyle{S_{m}(x)= \bigcup\limits_{i=0}^{n-1} S_{m,i}(x)}[/math] которая определена на отрезке [math][a,b][/math], принадлежит классу гладкости [math]C_r[a,b][/math] и составлена из совокупности звеньев [math]S_{m,i}(x),~ i=\overline{0,n-1}[/math], определенных на каждом частичном отрезке [math][x_i,x_{i+1}][/math] сетки [math]\Omega_n[/math], называется алгебраическим интегрально-дифференциальным сплайном степени [math]m[/math] и дефекта [math]q~ (0 \leqslant r \leqslant m,~ q=m-r)[/math] с узлами на сетке [math]\Omega_n[/math], если каждое его звено [math]S_{m,i}(x)~ x\in[x_i,x_{i+1}],~ i=\overline{0,n-1}[/math], представляется в виде алгебраического многочлена степени [math]m\colon[/math]

с коэффициентами [math]a_{k,i}[/math], выражаемыми из совокупности [math](m+1)[/math] интегральных и (или) дифференциальных условий согласования:

и из условий стыковки (непрерывности) звеньев [math]S_{m,i}(x)[/math] по производным [math]S_{m,i}^{(p_2)}(x)\colon[/math]

Здесь [math]p_1~(0 \leqslant p_1 \leqslant r)[/math] принимает значения из совокупности нескольких целых чисел, соответствующих порядкам производных, а [math]p_2[/math] — целое число или в общем случае несколько целых чисел, таких, что [math]0 \leqslant p_2 \leqslant r[/math], причем ни одно из этих чисел не равно ни одному значению чисел [math]p_1[/math].

Замечания

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

2. Условия стыковки (4.67) вместе с дифференциальными условиями согласования (4.66) обеспечивают непрерывность [math]S_{m}^{(p)}(x),~ p=0,1,\ldots,r[/math], во всех внутренних узлах [math]x_i[/math] определяющих точки стыковки звеньев, т.е. [math]\bigl(x_i, S_{m}(x_i)\bigr),~ i=\overline{1,n-1}[/math]. Это гарантирует выполнение условия [math]S_{m}(x)\in C_r[a,b][/math].

3. Алгоритм конструирования сплайна сочетает в себе два способа аппроксимации — кусочный и глобальный и поэтому является кусочно-глобальным. На кусочном способе основан процесс отыскания звеньев сплайна [math]S_{m,i}(x)[/math], т.е. коэффициентов [math]a_{k,i}[/math] многочлена (4.64), а на глобальном — соединение всех звеньев в единую (многозвенную) функцию [math]S_{m}(x),~ x\in[a,b][/math] со стыковкой звеньев так, чтобы обеспечивалась непрерывность [math]S_{m}^{(p)}(x)[/math] при [math]p=0,1,2,\ldots[/math], т.е. непрерывность как самой функции [math]S_{m}(x)[/math], так и ее производных.

4. Порядки производных [math]p_1[/math], выбранных для условий согласования (4.66), не должны совпадать с порядками производных, используемых в условиях стыковки (4.67), но вся совокупность порядков производных из [math]p_1[/math] и [math]p_2[/math] должна обеспечивать выполнение условия [math]S_{m}(x)\in C_r[a,b][/math].

5. Интегральное условие согласования (4.65) может быть переписано в виде [math]\textstyle{\int\limits_{x_i}^{x_{i+1}} S_{m,i}(x)dx= I_{i}^{i+1}}[/math], и тогда [math]I_{i}^{i+1}[/math] вместе со значениями [math]f^{(p_1)}(x_i)[/math], входящими в дифференциальное условие согласования, причисляется к параметрам ИД-сплайна. Данные параметры могут быть заданы при формулировке задачи аппроксимации либо определены при ее решении.

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

Интерполяционные дифференциальные кубические сплайны

Рассмотрим задачу восполнения заданной сеточной функции [math]y_i=f(x_i),~ i=\overline{0,n},~ x_i\in[a,b][/math] на базе интерполяционных глобальных кубических дифференциальных сплайнов дефекта [math]q=1[/math], то есть [math]S_3(x)\in C_2[a,b][/math]. При этом предположим, что восполняемая функция достаточно гладкая.

Решить эту задачу можно с помощью двух алгоритмов восполнения исходной сеточной функции, различающихся выбором порядков производных, на основе которых записываются условия согласования (4.66). Первый способ, наиболее широко распространенный, соответствует выбору вторых производных [math]m_i= f''_i~(i=\overline{0,n})[/math], а второй— выбору первых производных [math]\overline{m}_i= f'_i[/math].

Оба способа основаны на решении параметрически прямой задачи и весьма близки по характеру вычислительных процедур.

Основные этапы реализации первого способа.

Первый этап. Для алгебраического многочлена (полинома)

относящегося к одному i-му звену сплайна, для которого [math]x\in[x_i,x_{i+1}][/math], выбираются два дифференциальных условия согласования (4.66) с порядком [math]p_1=\{0;2\}[/math], каждое из которых относится к концам отрезка [math][x_i,x_{i+1}]\colon[/math]

Выбор значения [math]p=0[/math] обусловлен тем, что отыскивается [math]S_3(x)[/math] — интерполяционный многочлен, а выбор производных порядка [math]p=2[/math] зависит от выбора способа. Условия (4.69), (4.70) задают параметры сплайна [math]f_i,~ m_i= f''(x_i),~ i=\overline{0,n}[/math], первые из которых (функциональные) являются определенными, а вторые (дифференциальные) — неопределенными. Подстановка в соотношения (4.69), (4.70) полинома (4.68) приводит к системе четырех линейных алгебраических уравнений относительно неизвестных коэффициентов [math]a_{k,i}~ (k=0,1,2,3)[/math]. Разрешая их, подставляя полученные формулы для [math]a_{k,i}[/math] в (4.68) и распространяя этот результат на все частичные отрезки, имеем:

где [math]h_{i+1}= x_{i+1}-x_i,~ \Delta f_i= f_{i+1}-f_i,~ \Delta m_i= m_{i+1}-m_i[/math].

Таким образом, найдена общая формула одного звена сплайна, которая выражается через определенные [math](f_i)[/math] и неопределенные [math](m_i)[/math] параметры [math](i= \overline{0,n})[/math]. Эти параметры обеспечивают непрерывность сплайна: [math]\Bigl.{S_{3,i-1}(x)}\Bigr|_{x=x_i}= \Bigl.{S_{3,i}(x)}\Bigr|_{x=x_i},~ i=\overline{0,n-1}[/math] и его вторых производных [math]\Bigl.{S''_{3,i-1}(x)}\Bigr|_{x=x_i}= \Bigl.{S''_{3,i}(x)}\Bigr|_{x=x_i},~ i=\overline{0,n-1}[/math], для всех внутренних узлов.

Второй этап. При вычислении неопределенных параметров [math]m_i~ (i=\overline{0,n})[/math], входящих в (4.71), с учетом условия стыковки (4.67) для всех внутренних узлов устанавливается связь между [math]f_i[/math] и [math]m_i[/math],. Эта связь получается путем записи правой части (4.71) для звена [math]S_{3,i-1}(x)[/math] при [math]x\in [x_{i-1},x_i][/math] ee дифференцирования и приравнивания производной правой части (4.71) для звена [math]S_{3,i}(x)[/math] при [math]x\in[x_i,x_{i+1}][/math] в соответствии с условием стыковки (4.67) [math](p_2=1)[/math]. Проведя несложные выкладки и распространив найденное соотношение на все внутренние узлы, получим параметрическую систему — трехдиагональную систему линейных алгебраических уравнений, устанавливающую линейную связь комбинаций определенных [math]f_i[/math] и неопределенных параметров [math]m_i[/math], во внутренних узлах:

Данная система относительно [math]m_i~ (i=\overline{0,n})[/math] является незамкнутой, так как не хватает двух уравнений. Для ее замыкания используют различные пути выбора аппроксимаций производных на концах (граничных условий):

1. В простейшем случае вторые производные на концах принимаются нулевыми, т.е. [math]m_0=0;~ m_n=0[/math]. Такие условия называются условиями натурального сплайна.

Из рис. 4.11 видно, что задание условий [math]m_0=0[/math] и [math]m_n=0[/math] приводит к разрыву вторых производных [math]\Bigl.{f''(x)}\Bigr|_{x=x_0}[/math] и [math]\Bigl.{f''(x) }\Bigr|_{x=x_n}[/math] на концах отрезка [а,Ь], что может вызвать возрастание погрешности интерполирования.

2. Для первых двух и последних двух отрезков можно применить условие равенства третьей производной. Тогда получается соотношение [math]\frac{\Delta m_{k-1}}{h_k}= \frac{\Delta m_k}{h_{k+1}}[/math] [math](k=1,n-1)[/math]. Отсюда следуют два трехточечных уравнения:

Соотношения (4.73) позволяют замкнуть систему (4.72) при [math]h= \text{var}[/math] (неравномерная сетка). При [math]h=\text{const}[/math] вторые производные на концах вычисляются по явным формулам второго порядка аппроксимации:

Отметим, что как (4.73), так и формулы (4.74) для [math]m_0[/math] и [math]m_n[/math] имеют порядки аппроксимации, соответствующие порядку сходимости [math]S_3(x)[/math] к [math]f(x)[/math].

Третий этап. Определяются значения [math]m_i~ (i=\overline{0,n})[/math] путем решения систем (4.72), (4.73) (параметрически прямая задача). Для этого используется метод прогонки, причем необходимо предварительно из первых двух и последних двух уравнений исключить слагаемые с [math]m_2[/math] и [math]m_{n-2}[/math] соответственно. В случае равномерной сетки, как отмечено выше, система (4.72) замыкается аппроксимацией (4.74).

Четвертый этап. Формируются массивы коэффициентов [math]a_{0,i}, a_{1,i}, a_{2,i}, a_{3,i}[/math], для всех звеньев сплайна [math](i=\overline{0,n-1})[/math] и определяется глобальный интерполяционный сплайн [math]\textstyle{S_3(x)= \bigcup\limits_{i=0}^{n-1} S_{3,i}(x)}[/math] путем составления многозвенной функции.

Пятый этап. При необходимости полученная сплайн-функция [math]S_3(x)[/math] используется для вычисления значений функции и производных порядка [math]p\colon[/math] [math]f^{(p)}(x)= S_3^{(p)}(x)+ O(h_{i+1}^{4-p})[/math] [math](p=0,1,2,\ldots)[/math] в произвольных точках [math]x\in [x_i,x_{i+1}],~ i=\overline{0,n-1}[/math], или определенных интегралов [math]\textstyle{I_a^b= \int\limits_{a}^{b} S_3(x)dx}[/math] (в приложениях теории приближений).

Основные этапы реализации второго способа.

Первый этап. Для алгебраического многочлена (4.68), обозначаемого [math]\overline{S}_{3,i}(x)[/math], вместо вторых производных в (4.70) включаются первые производные. Тогда эти условия определяют параметры сплайна [math]f_i,~ \overline{m}_i= f'(x_i),~ i=\overline{0,n}[/math], первые из которых (функциональные) являются определенными, а вторые (дифференциальные) — неопределенными [math](p=1)[/math]. Проделав выкладки (аналогично первому этапу первого способа), получим общую формулу i-го звена искомого сплайна:

которая зависит от определенных [math](f_i)[/math] и неопределенных [math](\overline{m}_i= f'_i)[/math] параметров сплайна. Данные параметры обеспечивают непрерывность всего сплайна и его производных

Второй этап. С использованием условий стыковки (4.67) по вторым производным [math](p_2=2)[/math] во внутренних узлах получается система линейных алгебраических уравнений трехдиагонального типа, выражающая параметрическую связь между [math]\overline{m}_{i}[/math] и [math]f_i\colon[/math]

Замыкание системы (4.76) осуществляется аналогично п.2 на втором этапе реализации первого способа. При [math]h=\text{var}[/math] получаются два соотношения (граничные условия), не нарушающие трехдиагональность системы (где [math]H_{k-1}^k= h_{k-1}+ h_k[/math]):

Соотношения (4.77) аппроксимируют производные на концах с третьим порядком точности, что соответствует четвертому порядку сходимости [math]\overline{S}_3(x)[/math] к [math]f(x)[/math] (см. п.2 замечаний ниже). При [math]h=\text{const}[/math] для определения значений производных на концах можно использовать соответствующие аппроксимационные формулы.

Третий этап. Трехдиагональная система линейных алгебраических уравнений (4.76), (4.77) при конструировании глобального сплайна решается методом прогонки (параметрически прямая задача).

Четвертый и пятый этапы аналогичны изложенным для первого способа.

Замечания

1. Существование и единственность глобальных сплайнов [math]S_3(x)[/math] и [math]\overline{S}_3(x)[/math] следуют из рассмотренного выше построения общих формул их звеньев [math]S_{3,i}(x),\, \overline{s}_{3,i}(x)[/math] и единственности решения трехдиагональных систем (4.72), (4.73) или (4.74) и (4.76), (4.77) в силу выполнения условия диагонального преобладания.

2. Сходимость процесса построения глобальных сплайнов [math]S_{3,i}(x),\, \overline{s}_{3,i}(x)[/math] доказывается для непрерывных (формульных) функций, таких, что [math]f(x)\in C_4[a,b][/math]. Процесс является сходящимся, если при неограниченном увеличении числа [math]n[/math] узлов сетки соответствующая последовательность сплайн-функций, построенных на этих сетках, сходится к исходной функции [math]f(x)[/math]. Доказано, что на равномерной сетке сплайн-функции [math]S_{3,i}(x),\, \overline{s}_{3,i}(x)[/math] сходятся к [math]f(x)\in C_4[a,b][/math] с четвертым порядком, причем справедливы оценки

Таким образом, при использовании сплайн-функций для вычисления производных также реализуется сходимость, но ее порядок понижается на величину [math]p[/math].

3. В соответствии с принципом соответствия порядков аппроксимации точность вычисления интегралов на частичном отрезке [math][x_{i},x_{i+1}][/math], то есть [math]\textstyle{I_{i}^{i+1}= \int\limits_{x_i}^{x_{i+1}} f(x)dx \approx \widehat{I}_{i}^{i+1}= \int\limits_{x_i}^{x_{i+1}} S_3(x)dx}[/math], на основе сплайнов [math]S_{3,i}(x),\, \overline{s}_{3,i}(x)[/math] составляет [math]O(h_{i+1}^5)[/math], а интеграла [math]\textstyle{I_a^b= \int\limits_{a}^{b} f(x)dx}[/math] на отрезке [math][a,b][/math] составляет [math]O(H^4),~ H=\max_{i}h_{i+1}[/math] (потеря порядка происходит из-за суммирования значений [math]\widehat{I}_{i}^{i+1}[/math] и остаточного слагаемого при переходе от [math]\widehat{I}_{i}^{i+1},~ i=\overline{0,n-1}[/math] к [math]\textstyle{\widehat{I}_{a}^{b}= \sum\limits_{i=0}^{n-1} \widehat{I}_{i}^{i+1}}[/math].

4. Выражая из системы (4.72) комбинацию значений функций [math]f_{i-1},\, f_{i},\, f_{i+1}[/math] через комбинацию [math]m_{i-1},\, m_{i},\, m_{i+1}[/math], получаем трехдиагональную систему линейных алгебраических уравнений относительно [math]f_i\colon[/math]

Если считать, что вторые производные для некоторой функции [math]y_i=f(x_i)[/math] заданы или каким-либо способом определены, а [math]f_0[/math] и [math]f_n[/math] на концах [math][a,b][/math] известны, то из (4.79) можно восстановить значения [math]\overline{f}_i[/math] во внутренних узлах [math]x_i~ (i=\overline{1,n-1})[/math] путем решения параметрически обратной задачи. В этом случае сплайн [math]S_3(x)[/math], построенный по [math]\overline{f}_i,~ m_i~ (i=\overline{0,n})[/math] и (4.71), является в строгом понимании не интерполяционным, а восстанавливающим функцию [math]y_i=f(x_i)[/math]. Единственность решения задачи обеспечивается преобладанием диагональных коэффициентов в системе (4.79) относительно значений [math]f_i~ (i=\overline{0,n})[/math].

5. На основе соотношения (4.76) при заданных двух значениях функций [math]f_0,\,f_i[/math] или [math]f_{n-1},\, f_n[/math] и всех значениях [math]\overline{m}_i~ (i=\overline{0,n})[/math] можно восстановить функцию [math]\widetilde{f}(x_i)[/math] во всех точках сетки [math]x_i~ (i= \overline{2,n})[/math]. С этой целью система (4.76) разрешается относительно приращений функций [math]\Delta f_i[/math], которые затем вычисляются рекуррентно. При задании неравномерной сетки устойчивость процесса достигается при регулярном сгущении сетки в направлении слева направо [math]\left(\delta_{i+1}= \frac{h_{i+1}}{h_i}\leqslant 1\right)[/math] или при ее регулярном разрежении [math](\delta_{i+1}>1)[/math]. В первом случае восстановление осуществляется слева направо при заданных [math]f_0,\, f_1\colon[/math]

где [math]L_h= \frac{1}{3} \left[\frac{1}{h_i} \overline{m}_{i-1}+ 2\! \left(\frac{1}{h_i}+ \frac{1}{h_{i+1}}\right)\!\overline{m}_{i}+ \frac{1}{h_{i+1}} \overline{m}_{i+1}\right][/math], а во втором — справа налево при заданных [math]f_{n},\, f_{n-1}[/math]. Такой алгоритм обеспечивает устойчивость процесса обратной прогонки.

6. Локальный способ построения интерполяционных кубических сплайнов для заданной функции [math]y_i= f(x_i)[/math] основан на многочлене (4.75) и состоит в определении [math]\overline{m}_{i},~ i=\overline{0,n}[/math], по аппроксимационным формулам третьего порядка, в дальнейшей их подстановке в (4.75) и в формировании [math]\textstyle{\overline{S}_3(x)= \bigcup\limits_{i=0}^{n-1} \overline{S}_{3,i}(x)}[/math]. Необходимость аппроксимации третьего /"о порядка следует из принципа соответствия порядков, изложенного ранее, а также из параметрической связи (4.76).

Алгоритм построения интерполяционного кубического сплайна

1. Используя систему (4.72) , условия [math]m_0=0,~m_n=0[/math] или (4.73) (при [math]h= \text{const}[/math] можно использовать условия (4.74)), сформировать замкнутую систему относительно неопределенных параметров.

2. Вычислить коэффициенты системы (4.72), зависящие от шагов сетки [math]\Omega_n[/math], и решить ее методом прогонки. В результате найти неопределенные параметры [math]m_0,m_1,\ldots,m_n[/math].

3. Записать выражения для звеньев сплайна по формуле (4.71):

вычислить коэффициенты звеньев и сформировать многозвенную сплайн-функцию [math]\textstyle{S_3(x)= \bigcup\limits_{i=0}^{n-1} S_{3,i}(x)}[/math], являющуюся сплайном дефекта [math]q=1[/math].

Интерполяционные дифференциальные параболические сплайны

Выше отмечено, что применение параболических сплайнов более оправданно при аппроксимационной обработке широкого класса сеточных функций, например, рассчитанных методами сеток второго и третьего порядка. Сплайн-методы развиваются и как самостоятельный класс методов решения задач с обыкновенными дифференциальными уравнениями. Следует отметить, что широко распространенные на практике кубические дифференциальные сплайны теоретически обеспечивают четвертый порядок сходимости. Однако для их корректного применения необходимо, чтобы исходные данные (аппроксимируемая функция [math]y_i= f(x_i)[/math]) были заданы с порядком точности, не меньшим третьего. Это условие при обработке рассчитанных на основе методов сеток функций [math]y_i= f(x_i)[/math], как правило, не выполняется, и поэтому применять кубические сплайны не всегда целесообразно. Это обусловливает необходимость использования параболических сплайн-функций. Однако полученные традиционным способом сплайны [math]S_2(x)[/math] неустойчивы в смысле отсутствия их сходимости к исходной функции, что связано с несимметричным характером задания условий согласования и, как следствие, с аналогичным характером граничных условий. Поэтому для регуляризации сплайнов [math]S_2(x)[/math] Шенберг в построил такой алгоритм аппроксимации с помощью Д-сплайнов, в котором узлы сплайна сдвинуты относительно узлов сеточной функции. Для нахождения параметров этого сплайна глобальным способом получается трехдиагональная система линейных алгебраических уравнений с преобладанием диагональных элементов. Однако такая модификация существенно усложняет расчетные формулы алгоритма, и в ней используется дополнительное, обычно в задачах не требующееся условие непрерывности производных в серединах отрезков между узлами интерполяции. В ряде случаев это дополнительное требование может привести к непредсказуемому снижению точности аппроксимации и к отсутствию сходимости.

Рассмотрим принцип построения параболических дифференциальных сплайнов со сдвигом узлов сплайна. Пусть исходная функция [math]y_i= f(x_i),~ i=\overline{0,n}[/math], задана в узлах [math]x_i\in [a,b][/math], называемых узлами интерполяции и принадлежащими некоторой сетке [math]\Delta\colon\, a=x_0<x_1< \ldots< x_n< b~ (n\geqslant2)[/math]. Введем также в рассмотрение другую сетку [math]\overline{\Delta}\colon\, a=\overline{x}_0<\overline{x}_1< \ldots< \overline{x}_n< \overline{b}= \overline{x}_{n+1}[/math], узлы которой сдвинуты по отношению к узлам сетки [math]\Delta[/math] влево, так что [math]x_{i-1}< \overline{x}_i< x_i~ (i=\overline{1,n})[/math].

Функция [math]\textstyle{S_2(x)= \bigcup\limits_{i=0}^{n-1} S_{2,i}(x)}[/math] называется интерполяционной параболической сплайн-функцией, интерполирующей [math]f(x_i)[/math] на сетке [math]\Delta[/math] и принадлежащей классу гладкости [math]C_1[a,b][/math], если она составлена из совокупности звеньев [math]S_{2,i}(x)~ (i=\overline{0,n+1})[/math], определенных на каждом частичном отрезке [math][\overline{x}_i, \overline{x}_{i+1}][/math] сетки [math]\overline{\Delta}[/math] и представленных в виде алгебраических многочленов второй степени:

где [math]a_{k,i}[/math] — коэффициенты, которые находятся из дифференциальных условий согласования

где [math]p_1=0[/math] или [math]p_1=1[/math]. Звенья затем стыкуются по условию непрерывности (стыковки):

с порядками [math]p_2[/math], отличными от [math]p_1[/math], но такими, что [math]0\leqslant p_2 \leqslant2[/math].

В соответствии с этим определением приведем два сплайна [math]S_2(x)[/math] и [math]\overline{S}_2(x)[/math] и рассмотрим алгоритм построения одного из них , а именно [math]S_2(x)[/math]. Для некоторого отрезка [math][\overline{x}_i, \overline{x}_{i+1}][/math] формула [math]S_{2,i}(x)[/math] имеет вид (4.80)

На отрезке [math][\overline{x}_i, \overline{x}_{i+1}][/math] для определения [math]a_{0,i},\, a_{1,i},\, a_{2,i}[/math] используются условия согласования (4.81):

Для простоты предположим, что узлы сплайна расположены в серединах отрезков [math][x_i, x_{i+1}],~ i=\overline{0,n-1}[/math], составленных по узлам интерполяции.

На рис. 4.12 изображена схема расположения узлов сетки [math]\overline{\Delta}[/math] по отношению к узлам сетки [math]\Delta[/math]. Между шагами сеток [math]\Delta[/math] и [math]\overline{\Delta}[/math] в этом случае существует простая связь [math]\overline{h}_{i+1}= \frac{h_i}{2}+ \frac{h_{i+1}}{2}= \frac{1}{2}H_{i}^{i+1}[/math].

Из условия [math]\delta S'_{2,i}(\overline{x}_i)=0[/math] следует равенство [math]a_{1,i}= \overline{m}_i~ (\overline{m}_i= f'(\overline{x}_i))[/math], а из условия [math]\delta S'_{2,i} (\overline{x}_{i+1})=0[/math] вытекает, что [math]\overline{m}_i+ 2a_{2,i} \frac{1}{2}H_{i}^{i+1}= \overline{m}_{i+1};~ a_{2,i}= \frac{\Delta \overline{m}_i}{H_{i}^{i+1}}[/math]. Из условия интерполяции [math]\delta S_{2,i}(x_i)=0[/math] определяется коэффициент [math]a_{0,i}= f_i-\frac{h_i}{2}\overline{m}_i-\frac{h_i^2 \Delta \overline{m}_i}{4H_{i}^{i+1}}[/math].

Итак, звено сплайна [math]S_{2,i}(x)[/math] на произвольном отрезке [math][\overline{x}_i, \overline{x}_{i+1}][/math] выражается формулой

где [math]\delta_k= h_{k+1}\!\!\not{\phantom{|}}\, h_k~ (k=1,2,\ldots,i,i+1,\ldots)[/math].

Подставив (4.83) в (4.82), получим параметрическую систему относительно [math]\overline{m}_i\colon[/math]

где [math]H_{2(i-1)}^{i}= 2h_{i-1}+h_i[/math]. Система (4.84) обладает свойством преобладания диагональных элементов, но она незамкнута, так как число неопределенных параметров равно [math]n+2[/math]. В качестве двух дополнительных уравнений (граничных условий) могут быть приняты значения производных, аппроксимированных на концах со вторым порядком аппроксимации:

где [math]\overline{H}_{n}^{n+1}= \overline{h}_n+ \overline{h}_{n+1},~ \overline{f}_i= f(\overline{x}_i)[/math]. Здесь шаги [math]\overline{h}_1, \overline{h}_2, \overline{h}_{n-1}, \overline{h}_n, \overline{h}_{n+1}[/math] и их комбинации взяты на сетке [math]\overline{\Delta}[/math]. После определения всех параметров [math]\overline{m}_i[/math] их подставляют во все звенья [math]S_{2,i}(x),~ i=\overline{0,n}[/math], из которых затем формируется результирующая многозвенная сплайн-функция [math]S_2(x)[/math], имеющая минимальный дефект [math]q=1[/math].

При [math]h=\text{const}[/math] формулы (4.83), (4.84), (4.85) упрощаются:

Здесь вместо сетки [math]\overline{\Delta}[/math] используется равномерная сетка [math]\overline{\Delta}_{\mathsf{d}}= \bigl(\overline{x}_{0 \mathsf{d}}, \overline{x}_1, \ldots, \overline{x}_n, \overline{x}_{n+1 \mathsf{d}}\bigr)[/math], где [math]\overline{x}_{0 \mathsf{d}}=a-\frac{h}{2},[/math] [math]\overline{x}_{n+1 \mathsf{d}}= b+\frac{h}{2}[/math] — дополнительные узлы; [math]\overline{f}_{0 \mathsf{d}}= f(\overline{x}_{0 \mathsf{d}}),~ \overline{f}_{n+1 \mathsf{d}}= f(\overline{x}_{n+1 \mathsf{d}})[/math].

Другой параболический сплайн [math]\overline{S}_{2,i}(x)[/math] интерполяционного типа приведен некоторыми авторами, и он построен для более общего случая, когда узлы [math]\overline{x}_i[/math] не совпадают с серединами отрезков [math][x_i, x_{i+1}][/math]. При этом [math]\overline{S}_{2,i} (x)[/math] записывается для отрезка [math][x_i, x_{i+1}][/math] с узлами интерполяции [math]x_i,\, x_{i+1}\colon[/math]

где, как и раньше, [math]x_{i-1}< \overline{x}_i< x_i~ (i=\overline{1,n})[/math], а

Здесь коэффициенты при [math](x-x_0)^0[/math] и [math](x-x_i)^1[/math] найдены с учетом условий согласования [math]\delta \overline{S}_{2,i}(x_i)=0[/math] и [math]\delta \overline{S'}_{2,i} (x_i)=0[/math]. Коэффициенты [math]c_i[/math] и [math]d_i[/math] определяются из аналогичных условий согласования, но заданных в точке [math]x_{i+1}[/math] отрезка [math][x_i, x_{i+1}]\colon[/math]

Тогда для получения связи определенных и неопределенных параметров используется условие стыковки (4.82) при [math]p_2=2[/math], т.е. условие непрерывности второй производной в точке [math]x_i\colon~ \Bigl.{\overline{S}''_{2,i-1}(x)}\Bigr|_{x=x_i}= \Bigl.{\overline{S}''_{2,i}(x)}\Bigr|_{x=x_i}[/math]. Учтем, что [math]\Bigl.{\overline{S}''_{2,i-1}(x)}\Bigr|_{x=x_i}= c_{i-1}+d_{i-1}[/math], a [math]\Bigl.{\overline{S}''_{2,i}(x) }\Bigr|_{x=x_i}=c_i[/math] в силу (4.88).

Подставляя в выражение [math]c_{i-1}+d_{i-1}=c_i[/math] соотношения (4.89), получаем систему линейных алгебраических уравнений относительно параметров [math]\overline{m}_i~ (i=\overline{1,n-1})[/math].

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

Отметим, что алгоритм построения рассмотренного сплайна является достаточно сложным, так как [math]\overline{S}_{2,i}(x)[/math] имеет дополнительное слагаемое [math]d_i(x-\overline{x}_{i+ 1})_{+}^2[/math], а коэффициенты сплайна и системы (4.90) зависят от шагов [math]h[/math] и [math]\overline{h}[/math], характеризующих две сетки [math]\Delta[/math] и [math]\overline{\Delta}[/math], которые также достаточно сложны. Кроме того, требование непрерывности второй производной, присущее сплайну второй степени, может привести не к улучшению качества аппроксимации, а к ее ухудшению, так как эта производная для [math]\overline{S}_{2,i}(x)[/math] является константой.

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

Замечание. Использование интегрального подхода (ИД-сплайна) позволяет строить устойчивые алгоритмы аппроксимации функций параболическими сплайнами без смещения узлов сплайна относительно узлов интерполяции. Кроме того, эти сплайны являются консервативными, в том смысле, что интегралы от [math]S_2(x)[/math] и [math]y_i= f(x_i)[/math] равны как на [math][a,b][/math], так и на всех частичных отрезках [math][x_i, x_{i+1}]~ (i=\overline{0,n-1})[/math]. С помощью интегрального подхода могут быть построены замкнутые алгоритмы аппроксимации сплайнами произвольной четной степени.

Восстанавливающие, интерполяционные и сглаживающие интегрально-дифференциальные параболические сплайны

Ранее были приведены два интегрально-дифференциальных многочлена (4.37) (4.41), которые могут использоваться для интерполяции и сглаживания некоторых функций, заданных своими значениями и интегралами, либо интегралами и производными.

Здесь рассмотрим основанные на этих многочленах глобальные способы конструирования различных параболических ИД-сплайнов дефекта [math]q=1[/math], позволяющие восстанавливать некоторую функцию либо ее сглаживать. Рассмотрим применение многочлена (4.37).

Построение восстанавливающих сплайнов (задача 1)

Пусть некоторая сеточная функция определена не значениями [math]f_i[/math] в узлах [math]x_i[/math] в общем случае неравномерной сетки [math]\Omega_n[/math], а значениями интегралов [math]I_{i}^{i+1}~ (i=\overline{0,n-1})[/math].Требуется восстановить эту функцию путем построения сплайн-функции [math]\widetilde{S}_{2,\mathsf{id1}}(x)[/math] дефекта [math]q=1[/math].

Для решения этой задачи воспользуемся параболическим многочленом (4.37). Для получения [math]\widetilde{S}_{2,\mathsf{id1}}(x)[/math] при известных параметрах [math]I_{i}^{i+1}~ (i=\overline{0,n-1})[/math] необходимо сначала определить функциональные параметры [math]f_i~ (i=\overline{0,n})[/math]. Сплайн, составленный из звеньев (4.37), при любых значениях [math]f_i[/math] обеспечивает непрерывность [math]\widetilde{S}_{2,\mathsf{id1}}(x)[/math]. Поэтому для получения параметрической связи между [math]f_i[/math] и [math]I_{i}^{i+1}[/math] используются условия стыковки [math]\widetilde{S}_{2,\mathsf{id1}}(x)[/math] по первым производным:

Преобразование (4.91) после подстановки многочлена (4.37) приводит к параметрической системе линейных алгебраических уравнений относительно [math]f_i\colon[/math]

Из сопоставления систем (4.92) и (4.76), являющихся следствием условий стыковки параболических ИД-сплайнов и кубических Д-сплайнов, видно, что они являются подобными. Действительно, обе системы имеют одинаковые коэффициенты, а параметры в них, входящие в левую и правую части, подобны друг другу, так как различаются только порядком [math]p[/math]. Если их представить в виде [math]f_k^{(p)}~ (k=i-1,i,i+1)[/math] и считать [math]p=-1[/math] для [math]I_{i}^{i+1}= \Delta F_i= F_{i+1}-F_i[/math] ([math]F(x)[/math] — первообразная, [math]F_i= F(x_i)[/math]), то все порядки параметров, входящих в систему (4.76), будут на единицу больше соответствующих порядков параметров системы (4.92). Рассмотренные соотношения соответствуют также подобным звеньям (4.75) и (4.37), относящимся к кубическим и параболическим сплайнам. Действительно, звено параболического сплайна может быть получено путем дифференцирования звена кубического сплайна (4.75) и замены параметров [math]\overline{m}_i,\, \Delta f_i,\, \Delta \overline{m}_i[/math] на [math]f_i,\, I_{i}^{i+1}= \Delta F_i,\, \Delta f_i[/math].

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

являющиеся подобными (4.77) и полученные разрешением системы (4.92) и соотношений [math]a_{2,0}= a_{2,1},~ a_{2,n}= a_{2,n-1}[/math], где [math]a_{2,k},~ k=0;1[/math] и [math]k=n;n-1[/math] —коэффициенты при [math](x-x_i)^2[/math] (последние равенства для коэффициентов справедливы в силу постоянства второй производной на двух соседних отрезках). Подчеркнем, что соотношения (4.93) не нарушают трехдиагональности системы (4.92) и имеют порядок аппроксимации, равный трем, что соответствует порядку сходимости выстраиваемых ИД-сплайнов.

Разрешение системы (4.92), (4.93) методом прогонки, достаточное условие устойчивости которого для данной системы выполняется безусловно, т.е. при любых разбиениях, позволяет вычислить [math]f_i~ (i=\overline{0,n})[/math] (параметрически прямая задача). Подстановка [math]f_i,\, I_{i}^{i+1}[/math] в формулу (4.37) для звеньев [math]i=\overline{0,n-1}[/math] определяет интегрально-дифференциальную сплайн-функцию [math]S_{2,\mathsf{id1}}(x)[/math] восстанавливающего типа.

Из подобия параболического ИД-сплайна [math]S_{2,\mathsf{id1}}(x)[/math] Д-сплайну [math]\overline{S}_3(x)[/math] и возможности его получения путем дифференцирования кубического Д-сплайна, а также на основе соотношения (4.78), устанавливающего порядок сходимости Д-сплайна [math]\overline{S}_3(x)[/math], можно получить оценку порядка сходимости параболического ИД-сплайна [math]S_{2,\mathsf{id1}}(x)\colon[/math]

свидетельствующую о сходимости [math]S_{2,\mathsf{id1}}(x)[/math] с третьим порядком [math](p=0)[/math], а производных [math]S_{2,\mathsf{id1}}^{(p)}(x)[/math] с порядком [math](3-p),~ p=1;2[/math].

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

Построение интерполяционных сплайнов (задача 2)

Пусть на равномерно сгущающейся или разрежающейся сетке [math]\Omega_n[/math] заданы значения сеточной функции [math]f_i=f(x_i),~ i=\overline{0,n}[/math]. Требуется восполнить эту функцию глобальным интерполяционным интегрально-функциональным сплайном [math]S_{2,\mathsf{id1}}^{(\mathsf{int})}(x)[/math] дефекта [math]q=1[/math].

Решение этой задачи осуществляется путем объединения многочленов-звеньев (4.37) после определения интегралов [math]I_{i}^{i+1}[/math], выраженных из (4.92) по рекуррентным право- и левосторонним соотношениям (параметрически обратная задача):

где [math]L_i(h,f)= \frac{1}{h_i}f_{i-1}+ 2\! \left(\frac{1}{h_i}+ \frac{1}{h_{i+1}}\right)\!f_i+ \frac{1}{h_{i+1}}f_{i+1}[/math] — функция, равная сумме "удельных интегралов" по каждой паре смежных отрезков, [math]\delta_{i+1}= h_{i+1}\!\!\not{\phantom{|}}\,h_i[/math] — параметр, характеризующий неравномерность сетки, который называется параметром нерегулярности. Соотношения (4.95), (4.96) записаны в виде, учитывающем характер задания сеточной функции на неравномерной сетке, которая, по предположению, либо равномерно сгущается вправо, либо равномерно разрежается. Тогда (4.95) соответствует правостороннему обратному ходу метода прогонки, а (4.96) — левостороннему обратному ходу. Необходимость такого порядка действий обусловлена достаточным условием устойчивости прогонки [math]\delta_{i+1}\leqslant1[/math] для (4.95) и [math]\delta_{i+1}^{-1}\leqslant1[/math] для (4.96). Заметим, что в данном случае выполняется только обратный ход метода прогонки. Соотношения (4.95), (4.96) при [math]f(x)\in C_m[a,b]~ (m \geqslant 3)[/math] обеспечивают четвертый порядок точности вычисления интегралов.

В качестве краевых (граничных) условий для право- и левостороннего соотношений значения интегралов [math]I_0^i[/math] и [math]I_{n-1}^n[/math] вычисляются по одноинтервальным, трехточечным, лево- и правосторонним квадратурным формулам:

При [math]h=\text{const}[/math] (сетка [math]\Omega_n[/math] равномерная) данные формулы упрощаются:

В некоторых работах показано, что квадратурные формулы (4.97), (4.98) обеспечивают четвертый порядок аппроксимации интегралов на трехточечном шаблоне [math](x_{i-1}, x_i, x_{i+1})[/math], и поэтому порядок аппроксимации [math]I_{i}^{i+1}[/math] соотношением (4.95) или (4.96) при фиксированном [math]i[/math] за счет краевых условий не понижается. Показано также, что в данном способе построения интерполяционного сплайна достигается третий порядок сходимости. Однако способ имеет недостаток, состоящий в том, что при аппроксимации негладких функций ошибки, возникающие при расчете интегралов по формулам (4.95) или (4.96), переносятся на последующие значения. Данный недостаток устранен в задаче построения слабо сглаживающего параболического сплайна.

Построение сглаживающих сплайнов (задача 3)

Пусть в общем случае на неравномерной сетке Пп задана сеточная функция, определенная, например, в физическом эксперименте с достаточно большими погрешностями, определяющими разброс значений [math]f(x_i)[/math], то есть [math]f(x_i)\pm \varepsilon_i~ (i=\overline{0,n})[/math], где е, намного превышают величину [math]O(h_{i+1}^3)[/math].

Требуется восполнить эту функцию параболическим сглаживающим сплайном [math]S_{2, \mathsf{id1}}(x)[/math] дефекта [math]q=1[/math].

Поставленная задача решается глобальным способом путем сведения ее к задаче 1. Рассмотрим его особенности.

Сначала путем визуального изучения исходной функции [math]y_i= f(x_i)\pm \varepsilon_i[/math] исследователь-вычислитель формирует "полосу разброса" данных, офаниченную двумя функциями — "верхней" [math]y_i^{\mathsf{v}}= \widehat{f}^{(\mathsf{v})}(x_i)[/math] и "нижней" [math]y_i^{\mathsf{n}}= \widehat{f}^{(\mathsf{n})}(x_i)[/math] (на рис. 4.13 квадратиками условно показаны значения, полученные, например, в результате физического эксперимента). Если на некоторых участках разброс данных не заметен, функции [math]y_i^{\mathsf{v}}= \widehat{f}^{(\mathsf{v})}(x_i),~ y_i^{\mathsf{n}}= \widehat{f}^{(\mathsf{n})}(x_i)[/math] могут быть совмещены друг с другом. В случае отсутствия значений [math]f(x_i)[/math] на концах отрезка [math][a,b][/math] они определяются приближенно путем линейной или параболической экстраполяции.

Далее на каждом отрезке [math][x_i,x_{i+1}],~ i=\overline{0,n-1}[/math], по формулам (4.97) или (4.98) вычисляются интегралы [math]I_{i}^{i+1(\mathsf{v})}[/math] и [math]I_{i}^{i+ 1(\mathsf{n})}[/math] для функций [math]y_i^{\mathsf{v}}= \widehat{f}^{(\mathsf{v})}(x_i)[/math] и [math]y_i^{\mathsf{n}}= \widehat{f}^{(\mathsf{n})}(x_i)[/math] соответственно. Недостающие значения функций на линиях, ограничивающих "полосу разброса", определяются путем дополнительной интерполяции. Следует также учесть, что интеграл [math]I_0^1[/math] на отрезке [math][x_0,x_1][/math] вычисляется по левосторонней формуле (4.97), интеграл [math]I_{n-1}^n[/math] на отрезке [math][x_{n-1},x_n][/math] — по правосторонней, а интегралы по всем внутренним отрезкам могут рассчитываться по любой из формул (4.97), (4.98).

После вычисления интегралов [math]I_{i}^{i+1(\mathsf{v})}[/math] и [math]I_{i}^{i+ 1(\mathsf{n})}~ (i=\overline{0,n-1})[/math] по всем частичным отрезкам находятся их средние значения на каждом частичном отрезке:

Полученные значения интегралов [math]I_{i}^{i+1}[/math] используются далее в качестве исходных данных для решения задачи 1.

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

Слабо сглаживающие интерполяционные интегрально-дифференциальные параболические сплайны

Пусть в общем случае на неравномерной сетке [math]\Omega_n[/math] задана сеточная функция [math]y_i= f(x_i),~ i=\overline{0,n}[/math], определенная либо в численном эксперименте, либо путем измерений с небольшими погрешностями [math]\varepsilon_i[/math], не превышающими [math]O(h_{i+1}^3)[/math].

Применительно к данной задаче вводится понятие слабого сглаживания, в котором учитывается характер задания исходной функции с некоторой погрешностью, не превышающей [math]O(h_{i+1}^3)[/math]. При этом допускается возможность отклонения получаемой в алгоритме сплайн-функции в узлах [math]\widetilde{S}_{2, \operatorname{id}k}^{(\text{sgl})}(x_i)[/math] от [math]f(x_i)[/math] на величину не более [math]O(h_{i+1}^3)[/math], т.е. реализуется отклонение [math]\max \bigl|\widetilde{S}_{2, \operatorname{id}k}^{(\text{sgl})}(x)-f(x)\bigr|= O(H^3)[/math] по всему отрезку [math][a,b][/math], включая и узлы интерполяции. Здесь индекс [math]k[/math] ([math]k=1[/math] или [math]k=2[/math]) соответствует способу построения сплайна. Подчеркнем, что для традиционных интерполяционных сплайнов данное отклонение в узлах отсутствует.

Требуется восполнить данную функцию параболическим интерполяционным слабо сглаживающим сплайном [math]\widetilde{S}_{2, \operatorname{id}k}^{(\text{sgl})}(x)~ (k=1;2)[/math], компонующимся на основе ИД-многочленов (4.37) или (4.41).

Обоснованность такого расширения понятия интерполяции (дополнение его понятием слабого сглаживания) обусловлена следующими причинами.

1. Интерполяционные сплайн-функции, как и классические интерполяционные многочлены, часто применяются для восполнения и обработки сеточных функций, рассчитанных методами сеток второго и третьего порядка, и потому данные функции (исходные для задачи аппроксимации) неточные, и их погрешность обычно составляет [math]O(h_{i+1}^3)[/math]. Поэтому применение условий интерполяции, обеспечивающих нулевые невязки в узлах, вообще говоря, не является оправданным.

2. Если даже значения исходной функции заданы точно (хотя класс задач, в которых значения [math]f(x_i)[/math] являются точными, указать трудно), то точность сплайн-интерполяции из-за этого не повышается. Для параболических сплайн-функций порядок сходимости оценивается как третий, а для кубических — как четвертый.

Построение слабо сглаживающего параболического ИД-сплайна может быть осуществлено так же, как и построение кубического дифференциального сплайна двумя способами, различающимися порядками производных, заданных в условиях согласования и условиях стыковки. Это соответствует отмеченному выше принципу подобия дифференциальных и интегрально-дифференциальных сплайнов, описывающихся формулами (4.75) и (4.37). Аналогично легко показать подобие Д-сплайнов (4.71) и (4.41).

Поэтому, так же как и для Д-сплайнов, здесь рассматриваются два сглаживающих интегрально-дифференциальных сплайна [math]\widetilde{S}_{2, \operatorname{id}1}^{(\text{sgl})}(x)[/math] и [math]\widetilde{S}_{2, \operatorname{id}2}^{(\text{sgl})}(x)[/math].

Первый способ решения задачи слабого сглаживания (построение сплайна [math]\widetilde{S}_{2, \operatorname{id}1}^{(\text{sgl})}[/math]). Алгоритм восполнения функции [math]y_i= f(x_i)[/math] методом слабого сглаживания аналогичен алгоритму восстановления сеточной функции в задаче 1 и основан на формуле ИД-многочлена (4.37), принимаемой в качестве общей формулы звена искомого сплайна и параметрической системы (4.92), которые имеют вид

Волнистая черта указывает на отличие используемых значений [math]\widetilde{f}_i[/math] от значений [math]f_i[/math] исходной функции. В правой части равенства (4.100) обе группы параметров [math]\widetilde{f}_i~ (i=\overline{0,n})[/math] и [math]I_{i}^{i+1}~ (i=\overline{0,n-1})[/math] считаются неопределенными. При этом сначала вычисляются интегралы [math]I_{i}^{i+1}[/math] на всех частичных отрезках [math][x_i, x_{i+1}]~ (i=\overline{0,n-1})[/math] по одной из явных квадратурных формул (4.97) или (4.98). После этого проводится пересчет функциональных параметров из системы (4.101), которая замыкается связями [math]\widetilde{f}_0, \widetilde{f}_1,\, I_0^1,I_1^2[/math] и [math]\widetilde{f}_{n-1}, \widetilde{f}_n,\, I_{n-2}^{n-1}, I_{n-1}^n[/math] (граничными условиями) соответственно на левом и правом краях отрезка [math][a,b]\colon[/math]

совпадающими с граничными условиями, принятыми в задаче 1. Отметим, что оба типа неопределенных параметров находятся путем решения смешанной параметрической задачи, так как сначала определялись интегралы [math]I_{i}^{i+1}~ (p=-1)[/math], а затем значения функции [math]\widetilde{f}_i~ (p=0)[/math].

Определение [math]\widetilde{f}_i~ (i=\overline{0,n})[/math] из (4.101), (4.102), подстановка [math]\widetilde{f}_i~ (i=\overline{0,n})[/math] и [math]I_{i}^{i+1}~ (i=\overline{0,n-1})[/math] в (4.100) для всех звеньев и их объединение в общую формулу [math]\widetilde{S}_{2, \operatorname{id}1}^{(\text{sgl})}(x)[/math] завершают решение задачи слабого сглаживания исходной приближенно заданной сеточной функции, погрешность определения которой в узлах не превышает [math]O(h_{i+1}^3)[/math]. Если же погрешность задания исходной сеточной функции больше [math]O(h_{i+1}^3)[/math], to погрешность вычисления интегралов превысит величину [math]O(h_{i+1}^3)[/math] и третий порядок точности аппроксимации [math]\widetilde{S}_{2, \operatorname{id}1}^{(\text{sgl})}(x)[/math] не будет достигнут.

Первый способ решения задачи аппроксимации слабо сглаживающим сплайном является сходным с решением задачи 2 о восполнении сеточной функции интерполяционным интегрально-функциональным сплайном [math]\widetilde{S}_{2, \operatorname{id}1}^{(\text{int})}(x)[/math]. Различие состоит в том, что для [math]\widetilde{S}_{2, \operatorname{id}1}^{(\text{int})}(x)[/math] неопределенными параметрами являются только [math]I_{i}^{i+1}~ (i=\overline{0,n-1})[/math], а определенными — значения [math]f_i~ (i=\overline{0,n})[/math], и они не изменяются при решении задачи. Таким образом, задача 2 решается проще, однако изложенный способ построения сглаживающего сплайна характеризуется большей устойчивостью, так как он основан на решении трехдиагональной системы линейных алгебраических уравнений, которая удовлетворяет условию преобладания диагональных элементов.

Теорема 4.3 (о сходимости глобального параболического ИД-сплайна [math]\widetilde{S}_{2, \operatorname{id}1}^{(\text{sgl})}(x)[/math]). Пусть функцию [math]f(x),~ x\in[a,b][/math], определенную в общем случае на неравномерной сетке [math]\Omega_n~ \left(h_i=\text{var},~ H=\max_{i=1,n}h_i\right)[/math] с параметром неравномерности [math]D= \max_{i=\overline{1,n}}h_i \!\!\not{\phantom{|}}\, \max_{i=\overline{1,n}}h_i[/math], аппроксимирует глобальный параболический слабо сглаживающий ИД-сплайн [math]\widetilde{S}_{2, \operatorname{id} 1}^{(\text{sgl})}(x)[/math]. Тогда, если [math]f(x)\in C_3[a,b][/math] и интегралы [math]I_{i}^{i+1}~ (i=\overline{0,n-1})[/math] находятся по квадратурным формулам (4.97) или (4.98), a [math]\widetilde{f}_i~ (i=\overline{0,n})[/math] — из трехдиагональной системы (4.101), (4.102), то справедливы оценки

где [math]p[/math] — порядок производной, [math]M_3= \max_{x\in[a,b]} \bigl|f'''(x)\bigr|[/math], а константы имеют значения

Таким образом, при [math]f(x)\in C_3[a,b][/math] слабо сглаживающие сплайны [math]\widetilde{S}_{2, \operatorname{id}1}^{(\text{sgl})}(x)[/math] равномерно сходятся к функции [math]f(x)[/math] на последовательности сеток

с возрастающим разбиением [math]k[/math] отрезка [math][a,b]\colon\, k_1=1, k_2=2, \ldots, k_n=n[/math] и т.д. по крайней мере с кубичной скоростью, а их производные [math]\widetilde{S}'_{2, \operatorname{id}1}^{(\text{sgl})} (x)[/math] с ростом [math]n[/math] равномерно сходятся к [math]f'(x)[/math] по крайней мере с квадратичной скоростью.

Второй способ решения задачи слабого сглаживания (построение сплайна [math]\widetilde{S}_{2, \operatorname{id}2}^{(\text{sgl})}(x)[/math]). Алгоритм восполнения функции [math]y_i= f(x_i)[/math] в данной задаче основан на формуле ИД-многочлена (4.41), принимаемой в качестве общей формулы звена искомого сплайна, и на трехдиагональной параметрической системе, вытекающей из условия непрерывности сплайна [math]\widetilde{S}_{2, \operatorname{id}2, i}^{(\text{sgl})}(x)\colon[/math]

где [math]\overline{m}_i= f'(x)[/math]. Подчеркнем, что формула (4.104) для [math]\widetilde{S}_{2, \operatorname{id}2, i}^{(\text{sgl})}(x)[/math] и параметрическая система являются подобными соответствующим соотношениям кубического дифференциального сплайна [math]S_3(x)[/math] (4-70 и (4.72), в которых порядки всех параметров повышены на единицу.

Аналогично первому способу будем считать, что все параметры, входящие в (4.104), являются неопределенными. Их значения находятся из смешанной параметрической задачи, в которой вначале по формулам (4.97) или (4.98) определяются интегралы [math]I_{i}^{i+1}~ (i=\overline{0,n-1})[/math], а затем с использованием системы (4.105) трехдиагонального типа методом прогонки — производные [math]\overline{m}_i~ (i=\overline{0,n})[/math]. Для замыкания (4.105) можно использовать два различных подхода.

В рамках первого подхода производные на концах отрезка [math][a,b][/math] определяются по аппроксимационным формулам второго порядка аппроксимации:

Во втором подходе при [math]h_{i+1}= \text{var}[/math] система (4.105) замыкается двумя трехточечными соотношениями:

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

В случае, когда [math]h=\text{const}[/math] (сетка равномерная), производные на концах могут быть определены также и путем интегральных аппроксимаций:

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

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

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

В заключение отметим, что для сглаживающего сплайна [math]\widetilde{S}_{2, \operatorname{id} 2}^{(\text{sgl})}(x)[/math] как и для сплайна [math]\widetilde{S}_{2, \operatorname{id} 1}^{(\text{sgl})}(x)[/math], имеет место третий порядок сходимости.