RU: Curves. Monotone Interpolation

Геометрические примитивы

Кривые: монотонная интерполяция

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

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

Пусть на отрезке $[a,b]$ задана сетка узлов $\{t_i\}$: $\quad a=t_0< t_1<\dots <t_{n}= b$ и значения фунции $f_i = f(t_i)$, $0 \le i \le n$. Интерполяционный сплайн $S(t)$ третьего порядка строится на основании следующих условий:

• На каждом отрезке $[t_i,t_{i+1}],\quad 0\le i <n$ функция $S(t)$ является многочленом третьей степени:

(1)

\begin{align} S(t)=S_i(t)=a_{i0} + a_{i1}(t-t_i)+ a_{i2}(t-t_i)^2+ a_{i3}(t-t_i)^3, \qquad t_i\leq t \leq t_{i+1} \end{align}

• Функция $S(t)$ принимает в точках сетки $t_i$ заданные значения $f_i$.
• Функция $S(t)$ имеет непрерывные первые и вторые производные.

Чтобы найти значения коэффициентов $a_{i0}$, $a_{i1}$, $a_{i2}$, $a_{i3}$, выразим их через значения функции $f_i$ и ее производной $f'_i$ в узлах сетки. Значения производной нам неизвестны, и мы найдем их из условия непрерывности второй производной. Обозначим через $f'_{i+1/2}=(f_{i+1}-f_i)/(t_{i+1}-t_i)=\Delta f_{i+1/2}/\Delta t_{i+1/2}$ — оценку (разностную) производной на интервале $[t_i, t_{i+1}]$. Легко проверить, что

(2)

\begin{align} a_{i0} &=f_i & a_{i2}&=\frac{3f'_{i+1/2}-f'_{i+1}-2f'_i}{\Delta t_{i+1/2}}, &\nonumber\\ a_{i1}&= f'_i & a_{i3}&=-\frac{2f'_{i+1/2}-f'_{i+1}-f'_i}{(\Delta t_{i+1/2})^2} \end{align}

Приравнивая значения вторых производных в точке $t_i$, приходим к следущему условию для значений производных во внутренних узлах сетки $(1\leq i \leq n-1)$:

(3)

\begin{align} \frac{1}{\Delta t_{i-1/2}}f'_{i-1}+2\left(\frac{1}{\Delta t_{i-1/2}}+\frac{1}{\Delta t_{i+1/2}}\right)f'_{i} +\frac{1}{\Delta t_{i+1/2}}f'_{i+1}=3\left(\frac{f'_{i-1/2}}{\Delta t_{i-1/2}}+\frac{f'_{i+1/2}}{\Delta t_{i+1/2}}\right) \end{align}

Такого типа системы уравнений, связывающие значения функции в трех последовательных узлах сетки (в данном
случае $f'_{i-1}$, $f'_{i}$ и $f'_{i+1}$) называются трехточечными. Матрица системы уравнений имеет ненулевые элементы только на трех диагоналях: на главной и двух побочных. Обычно для систем с трехдиагональной матрицей используется следующая форма записи:

(4)

\begin{align} A_{i}Y_{i-1}+C_{i}Y_{i}+B_{i}Y_{i+1}=F_{i}, \qquad 1\leq i \leq n-1 \end{align}

Здесь $Y_i$ — искомые значения, в данном случае $Y_i=f'_{i}$, а $A_{i}$, $B_{i}$, $C_{i}$, $F_{i}$,— коэффициенты уравнений:

(5)

\begin{align} A_{i}&=\frac{1}{\Delta t_{i-1}} & C_{i}&=2\left(\frac{1}{\Delta t_{i-1}}+\frac{1}{\Delta t_{i}}\right)= 2(A_{i}+B_{i})\\ B_{i}&=\frac{1}{\Delta t_{i}} & F_{i} &=3\left(\frac{f'_{i-1/2}}{\Delta t_{i-1/2}}+\frac{f'_{i+1/2}}{\Delta t_{i+1/2}}\right) \nonumber \end{align}

Как уже было отмечено выше, соотношения во внутренних точках дают $n-1$ уравнение для $n+1$ неизвестных $f'_{i}$. Чтобы замкнуть систему, нужно добавить еще два уравнения, которые обычно получают из условий в граничных узлах $i=0$ и $i=n$. Эти граничные условия могут быть различных типов и определяются постановкой задачи. Например, если из некоторых дополнительных соображений известны значения первых производных в граничных точках [$a=x_0$]] и $b=x_n$, т.е. заданы $f'_{0}$ и $f'_{n}$, то условия имеют вид (4) при $i=0$ и $i=n$, если положить

(6)

\begin{align} A_{0}&=0 & C_{0}&=1 & B_{0}&=0 & F_{0}&=f'_{0} \\ A_{n}&=0 & C_{n}&=1 & B_{n}&=0 & F_{n}&=f'_{n} \nonumber \end{align}

Вернемся к обсужденю монотонности интерполяционной функции. Будем говорить, что заданные значения функции монотонны в точке $t_i$, если выполняются условия $f_{i-1}\le f_i \le f_{i+1}$ или $f_{i-1}\ge f_i \ge f_{i+1}$. В этом случае разностные производные на соседних интервалах имеют одинаковый знак $\text{sign}(f'_{i-1/2}) = \text{sign}(f'_{i+1/2})$. К сожалению, из монотонности заданных значений не следует монотонность
интерполяционной функции, как показывает следующий типичный пример.

Пример. Пусть в точках $t_i = i$, $0 \le i \le 5$ (их всего шесть) заданы следующии значения функции: $f_i=0$ при $i=0, 1, 2$ и $f_i=1$ при $i=3, 4, 5$. Значеня функции монотонны. Построим аппроксимирующий сплайн, считая производные в граничныъ точках нулевыми.Решая выписанную выше систему уравненй, найдем хначения прозводных:

(7)

\begin{align} f'_1 &= -0.6, & f'_2 &= 1.2, & f'_3 &= -1.2, & f'_4 &= 0.6. \end{align}

и уравнение сплайна имеет коэффициенты:

(8)

\begin{align} a_{00} &=0, & a_{01} &=0, & a_{02} &= 0.6, & a_{03} &= -0.6, \\ a_{10} &=0, & a_{11} &=-0.6, & a_{12} &= 0, & a_{13} &= 0.6, \\ a_{20} &=0, & a_{21} &= 1.2, & a_{22} &=-0,6, & a_{23} &= 0,4 \\ a_{30} &=1, & a_{31} &= 1.2, & a_{32} &= -1.8, & a_{33} &= 0.6, \\ a_{40} &=1, & a_{41} &=-0.6, & a_{42} &=1.2, & a_{43} &= -0.6 . \end{align}

Найдем условия, при выполнении которых интерполяционный многочлен будет монотонным на интервале $[t_i, t_{i+1}]$. Ясно, что необходимым условием служит совпадение знаков левой и правой производных со знаком разностной производной, т.е. соотношение

(9)

\begin{align} \text{sign}(f'_i) = \text{sign}(f'_{i+1/2}) = \text{sign}(f'_{i+1}) \end{align}

Если это соотношение не выполнено, то многочлен на интервале имеет экстремум. Но выполнение этого условия не является достаточным, как показывает следующий пример.

Пример. Пусть на интервале $[0,1]$ слева функция равна 0, а производная 4, а справа функция равна 1, а прозводная 0. График показывает немонотонность такой функци, хотя достаточные условия выполнены.

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

(10)

\begin{align} S'_i(t)&= a_{i1}+ 2a_{i2}(t-t_i)+ 3a_{i3}(t-t_i)^2,\\ S''_i(t)&= 2a_{i2}+ 6a_{i3}(t-t_i). \end{align}

мы можем выделить следующе ситуации.

• Если $f_i=f_{i+1}$ т.е. $f'_{i+1/2}=0$, то необходимым и достаточным условием монотонности на интервале является соотношение $f'_i=f'_{i+1}=0$. В ситуациях далее предполагается, что $f'_{i+1/2}\neq 0$.
• Пусть $a_{i3}=0$. В этом случае производная линейна, ее значения заключены между $f'_i$ и $f'_{i+1}$ и необходимые условия являются достаточными. Условие $a_{i3}=0$ означает, что $f'_i+f'_{i+1}-2f'_{i+1/2}=0$. Поскольку по предположению $f'_{i+1/2} \neq 0$, введем следующие обозначения: $\alpha_i=f'_i/ f'_{i+1/2}$, $\beta_i=f'_{i+1}/ f'_{i+1/2}$ и запишем данное условие как $\alpha_i+\beta_i-2=0$. На плоскости $(\alpha, \beta)$ (см. рисунок ниже) это отрезок $AB$.
• Пусть теперь $f'_i+f'_{i+1}-2f'_{i+1/2} \neq 0$. В этом случае производная является параболой. Она выгнута вниз, если $f'_i+f'_{i+1}-2f'_{i+1/2} > 0$ и вверх если $f'_i+f'_{i+1}-2f'_{i+1/2} < 0$. Если $f'_i$ и $f'_{i+1}$ положительны и парабола выгнута вверх, то производная не обращается на интервале в нуль, т.е. экстремума нет. Аналогично, если $f'_i$ и $f'_{i+1}$ отрицательны и парабола выгнута вниз, то производная на интервале не обращается в нуль, экстремума нет и кривая монотонна. Оба эти случая объединяются достаточным условием монотонности: $\alpha_i+\beta_i-2<0$. На плоскости $(\alpha, \beta)$ (см. рисунок ниже) это треугольник $OAB$.
• Пусть теперь $\alpha_i+\beta_i-2>0$. Найдем точку экстремума $t_*$ параболы $S'(t)$. Кривая $S(t)$ будет монотонна, если точка экстремума не лежит на интервале $[t_i,t_{i+1}]$, т.е. дибо $t_*<t_i$, либо $t_*>t_{i+1}$. Первое условие эквивалентно неравенству $2\alpha_i+\beta_i-3\le 0$, а второе $\alpha_i+2\beta_i-3\le0$. На плоскости $(\alpha, \beta)$ эти области представлены треугольниками $ACD$ и $BCE$.

• Наконец, если точка эктремума производной лежит на интервале, то ее значение может иметь тот же знак, что и производные на границах, так что парабола не пересекает ось абсцисс, производная не обращается в нуль, т.е. функция $S(t)$ монотонна. Экстремум призводной равен(11)
  \begin{align} S'(t_*) = a_{i1} -\frac{a_{i2}^2}{3a_{i3}}=f'_{i-1/2}\phi(\alpha,\beta), \quad \text{где} \quad \phi(\alpha,\beta)= \alpha - \frac{(2\alpha+\beta-3)^2}{3(\alpha+\beta-2)} \end{align}

  Таким образом, условие монотонности есть условие выполнения неравенства $\phi(\alpha,\beta)\ge 0$. На плоскости $(\alpha, \beta)$ эта область ограничена эллипсом.
  

Выписанные выше условия довольно громоздки, но, как видно из рисунка, достаточным является выполнение довольно простого условия принадлежности точки квадрату $( 0 \le \alpha_i \le 3, 0 \le \beta_i \le 3)$. Таким образом, для построения монотонной интерполяционной кривой по монотонным данным, можно поступить следующим образом:

• Построить интерполяционный сплайн.
• Подправить значения производных в узлах таким образом, чтобы выполнялсь условия
  • $0 \le f'_i \le 3\min( f'_{i-1/2},f'_{i+1/2})$ для монотонно возрастающей функции,
  • $0 \ge f'_i \ge 3\max( f'_{i-1/2},f'_{i+1/2})$ для монотонно убывающей функции.

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

Дальнейшее обсуждение этого вопроса и примеры можно найти в работах [1],[2].

Bibliography

1. F. N. Fritsch, R. E. Carlson, Monotone piecewise cubic interpolation, SIAM J. Num. Anal., 17 (1980), pp. 238-246.

2. James M. Hyman, Accurate Monotonicity Preserving Cubic Interpolation, SIAM J. Sci. Stat. Comput., 4 (1983), pp. 645-654.