Теория и построение

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

Пусть на отрезке [a, b] требуется интерполировать некоторую функцию f (x) по известным значениям, r = 0,1,2, на сетке Д: a==b. Обычно для решения такой задачи используют сплайны пятой степени с узлами на сетке Д, но можно применить кубические сплайны с дополнительными узлами. Для этого введем на [a, b] еще одну сетку.

д = {,.

где .

Кубический сплайн S (x) дефекта 1 с узлами на сетке, удовлетворяющий условиям.

для (1).

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

Условия интерполяции (1) на концах промежутков [, дают N линейных систем, по 6 уравнений в каждой, для определения коэффициентов сплайна, решая которые получаем при.

с коэффициентами.

.

Оценка погрешности локальной интерполяции.

Приведем оценку точности приближения функции f (x) сплайном S (x) в зависимости от ее дифференциальных свойств. В данном случае константы в оценках будут зависеть от параметра. В силу сложности получаемых выражений приведены результаты только для фиксированного .

Теорема 1. Если S (x) интерполирует на сетке Д функцию f (x), то имеют место оценки.

(1).

где при даны в таблице 1.

Класс функций.

Класс функций.

;

;

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

Примеры.

Функция |sin (x)|, h = 0.2.

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

Та же функция, но с шагом h=0.025.

На графике видно, что изменение порядка шага дало значительно лучшие результаты, хотя ошибки вычисления всё равно могут присутствовать Код решения на C#.

class ExtraSpline: Builder.

{.

public double a1;

public double a2;

ExtraSplineSegment[] segments;

List points;

int n;

public ExtraSpline (int _n, List _points, double _a1, double _a2).

{.

n = _n;

a1 = _a1;

a2 = _a2;

points = _points;

segments = new ExtraSplineSegment[n-1];

for (int i = 0; i < n-1; ++i).

{.

segments[i]. x = points[i]. x;

segments[i]. f = points[i]. y;

segments[i]. f1 = points[i]. y1;

segments[i]. f2 = points[i]. y2;

//находим h.

segments[i]. h = points[i + 1]. x — points[i]. x;

//находим c1 и с2.

segments[i]. c1 = c_function (a1, a2, segments[i]. f, segments[i]. f1, segments[i]. f2, points[i + 1]. y, points[i + 1]. y1, points[i + 1]. y2, segments[i]. h);

segments[i]. c2 = c_function (a2, a1, segments[i]. f, segments[i]. f1, segments[i]. f2, points[i + 1]. y, points[i + 1]. y1, points[i + 1]. y2, segments[i]. h);

segments[i]. a = a_function (segments[i]. c1, segments[i]. c2, points[i + 1]. y2, segments[i]. f2, segments[i]. h);

}.

}.

public double Interpolate (double x).

{.

if (segments == null).

{.

return double. NaN; // Если сплайны ещё не построены — возвращаем NaN.

}.

int n = segments. Length;

ExtraSplineSegment s;

if (x <= segments[0]. x) // Если x меньше точки сетки x[0] - пользуемся первым эл-тов массива.

{.

s = segments[1];

}.

else if (x >= segments[n — 1]. x) // Если x больше точки сетки x[n — 1] - пользуемся последним эл-том массива.

{.

s = segments[n — 1];

}.

else // Иначе x лежит между граничными точками сетки — производим бинарный поиск нужного эл-та массива.

{.

int i = 0;

int j = n — 1;

while (i + 1 < j).

{.

int k = i + (j — i) / 2;

if (x <= segments[k]. x).

{.

j = k;

}.

else.

{.

i = k;

}.

}.

s = segments[j];

}.

double dx = x — s. x;

double res = 0;

res = (s.a * dx * dx * dx) / 6 + (s.f2 * dx * dx) / 2 + (s.f1 * dx) + s. f + (s.c1 * (dx — a1 * s. h) * (dx — a1 * s. h) * (dx — a1 * s. h)) + (s.c2 * (dx — a2 * s. h) * (dx — a2 * s. h) * (dx — a2 * s. h));

return (double)res;

}.

public double c_function (double a1, double a2, double fi, double fi1, double fi2, double fii, double fii1, double fii2, double hi).

{.

double c = 0;

c = (1/(a1*a2*(a2-a1)*hi*hi)) * (((fi-fii)/hi)+(((1+a2)*fi1+(1+a1)*fii1)/3) +(((a2*fi2) — (a1*fii2))*hi/6));

Console.Write («h:»);

Console.WriteLine (hi);

Console.Write («c:»);

Console.WriteLine (c);

return c;

}.

public double a_function (double c1, double c2, double fii2, double fi2, double h).

{.

double a = 0;

a = ((fii2 — fi2)/h — 6*a2*c1 — 6*a1*c2);

return a;

}.

}.

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

Теория и построение.

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

(1).

В узлах сетки.

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

Коэффициенты находятся методом прогонки приведенной выше системы, так как матрица системы трехдиагональна:

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

В частности, если третья производная функции терпит разрыв в некоторых точках отрезка [a, b], то для улучшения качества аппроксимации эти точки следует включить в число узлов интерполяции.

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

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

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

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

При наличии разрывов функции в некоторых точках отрезка [a, b] внесем следующие изменения в построение сплайна. Значение сплайна и/или его производных слева от точки разрыва принимаем равным значению интерполируемой функции в этом узле слева, а величину скачка относим к последующему сплайну. Таким образом, для нахождения коэффициентов в узлах, соответствующих разрывам, вместо (2−5) имеем следующую систему уравнений:

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

Примеры Для сравнения взята кусочно-заданная функция, которая при x2,4.

y = -x².

Это естественный кубический сплайн h = 0.2.

А это уже сплайн, с учетом скачка, по графикам видно, что осцилляции, хоть и не ушли, но заметно уменьшились.

Та же функция, но с шагом h = 0.025.

Естественный кубический сплайн.

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

Код решения на C#.

class SplineD: Builder.

{.

SplineSegment[] segments;

List points;

int n;

/// Конструктор SplineBuilder.

/// Количество узлов сетк.

/// Значения в узлах.

public SplineD (int _n, List _points, List _extremums).

{.

double D = 0;

n = _n;

points = _points;

foreach (PointExtremum p in _extremums).

{.

int k = n — 1;

for (int i = 0; i < k; ++i).

{.

if (p.p.x > points[i]. x && p.p.x < points[i + 1]. x).

{.

points.Insert (i+1,p.p);

k++;

break;

}.

}.

}.

n = points. Count;

segments = new SplineSegment[n];

for (int i = 0; i < n; ++i).

{.

segments[i]. x = points[i]. x;

segments[i]. a = points[i]. y + D;

}.

segments[0]. c = 0;

segments[n — 1]. c = 0;

double[] alpha = new double[n — 1];

double[] beta = new double[n — 1];

alpha[0] = beta[0] = 0.0;

for (int i = 1; i < n — 1; ++i).

{.

foreach (PointExtremum pex in _extremums).

{.

if (pex.p.x == points[i]. x).

{.

D = pex. D;

}.

else.

{.

D = 0;

}.

}.

double hi = points[i]. x — points[i — 1]. x;

double hi1 = points[i + 1]. x — points[i]. x;

double A = hi;

double C = 2.0 * (hi + hi1);

double B = hi1;

double F = 6.0 * ((points[i + 1]. y — points[i]. y + D) / hi1 — (points[i]. y — points[i — 1]. y) / hi);

double z = (A * alpha[i — 1] + C);

alpha[i] = -B / z;

beta[i] = (F — A * beta[i — 1]) / z;

}.

for (int i = n — 2; i > 0; —i).

{.

segments[i]. c = alpha[i] * segments[i + 1]. c + beta[i];

}.

for (int i = n — 1; i > 0; —i).

{.

double hi = points[i]. x — points[i — 1]. x;

segments[i]. d = (segments[i]. c — segments[i — 1]. c) / hi;

segments[i]. b = hi * (2.0 * segments[i]. c + segments[i — 1]. c) / 6.0 + (points[i]. y — points[i — 1]. y + D) / hi;

}.

}.

public List getPoints.

{.

get { return points; }.

}.

public SplineSegment[] getSpline.

{.

get { return segments; }.

}.

public double Interpolate (double x).

{.

if (segments == null).

{.

return double. NaN; // Если сплайны ещё не построены — возвращаем NaN.

}.

int n = segments. Length;

SplineSegment s;

if (x <= segments[0]. x) {.

s = segments[1];

}.

else if (x >= segments[n — 1]. x) {.

s = segments[n — 1];

}.

else.

{.

int i = 0;

int j = n — 1;

while (i + 1 < j).

{.

int k = i + (j — i) / 2;

if (x <= segments[k]. x).

{.

j = k;

}.

else.

{.

i = k;

}.

}.

s = segments[j];

}.

double dx = x — s. x;

return s. a + (s.b + (s.c / 2.0 + s. d * dx / 6.0) * dx) * dx;

}.

}.

Эффект Гиббса для сплайн — интерполяции.

Пусть — функция вида.

Рассмотрим кубический интерполяционный сплайн для. Пусть и — кубический сплайн с узлами), интерполирующий в точках, т. е.; в каждом является алгебраическим многочленом третей степени, удовлетворяющим краевым условиям или и (. Тогда в окрестности точки a имеет место эффект Гиббса, точнее.

Теорема. Для.

Следствие. Для имеемпри (число узлов).

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

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

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

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

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

Для кубического сплайна с дополнительными узлами эффект Гиббса отсутствует.

точка разрыва в середине интервала.

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

Примеры.

Естественный сплайн для параметрической функции, описанной выше, h=0.2.

А это сплайн с дополнительными узлами для той же функции.

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