MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Spline Toolbox

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

В этом разделе мы рассмотрим зависимость ошибки интерполяции при интерполировании кубическим фундаментальным сплайном непрерывных функций. Соответствующие теоретические результаты приведены, например, в классической книге Carl de Boor. A Practical Guide to Splines, Springer-Verlag 1978, перевод которой на русский язык был опубликован издательством "Радио и связь" в 1985г. В ней приведены ссылки на работы, в которых были получены соответствующие результаты.

Пусть нам дана некоторая функция ƒ(x), определенная на отрезке [a,b], и задана возрастающая последовательность узлов интерполирования

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

Фундаментальным кубическим сплайном, интерполирующим функцию ƒ(x) в этих точках, называется кубический сплайн s(x) такой, что:

  1. он совпадает с интерполируемой функцией в узлах интерполяции, т.е. s(xk) = ƒ(xk) для всех k = 1,2,...,n;
  2. и выполняются следующие граничные условия

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

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

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

Ошибка интерполяции функции

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

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

Для демонстрации данного утверждения проинтерполируем, к примеру, функцию ƒ(x) = sinx на отрезке [0.3π/2] фундаментальным кубическим сплайном по равномерно отстоящим узлам, уменьшая расстояние между соседними узлами и следя за поведением ошибки интерполяции.

Для построения фундаментального кубического сплайна используем функцию csape, входящую в состав Spline Toolbox. Эта функция позволяет конструировать кубический сплайн с произвольными граничными условиями. В нашем случае граничные условия таковы:

В функции csape задаются три входных аргумента:

  1. вектор с абсциссами узлов интерполяции;
  2. вектор, первый элемент которого является значением производной сплайна на левой границе отрезка, следующие элементы являются значениями интерполируемой функции в узлах, а последний элемент - значение производной сплайна на левой границе отрезка;
  3. вектор из двух чисел - порядков производных, для которых задаются значения, в нашем примере производные первого порядка, следовательно надо указать вектор [1, 1].

Функция csape возвращает сплайн в кусочно-полиномиальной форме в виде структуры, поля которой содержат коэффициенты, узлы, порядок (см. разд. "Интерполяция кубическими сплайнами, функция csape"). С этой структурой можно работать дальше, выполняя различные операции со сплайнами при помощи функций Spline Toolbox.

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

Будем менять расстояние h между узлами интерполяции в цикле, полагая на каждом k-ом шаге цикла h = 2k, и накапливать в массивах H и Err значения шага и соответствующие значения ошибки интерполяции.

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

Ниже приведен текст соответствующего m-файла SplineErr.m

% задаем функцию, границы отрезка и значения первых производных в них
f = inline('sin(x)');
a=0; 
b=3*pi/2; 
dfa = 1; 
dfb = 0;
%Инициализация массивов для значений шагов и ошибки
H=[];
Err=[];
% построение сплайнов в цикле и вычисление ошибки интерполяции
for k=2:12
    % задаем шаг, узлы интерполяции и значения функции в них
    n = pow2(k);  
    h = (b-a)/(n-1);    
    x=a:h:b;
    y=f(x);
    % строим фундаментальный кубический сплайн
    pp=csape(x,[dfa y dfb],[1 1]);
    % находим максимальное уклонение по большому набору точек
    xx=linspace(a,b,30*n);
    yy=fnval(pp,xx);
    err=max(abs(yy-f(xx)));    
    % пополняем массивы со значениями шага и ошибки
    H = [H h];
    Err = [Err err];
end
% создаем графическое окно и оси
figure
axes
% строим зависимость ошибки и оценки от шага в логарифмической шкале
loglog(H, Err, 'LineWidth', 2, 'Color', 'r')
hold on
loglog(H, 5/384*H.^4, 'LineWidth', 1, 'Color', 'b')
% задаем шаги сетки по координатным осям
set(gca, 'YTick',logspace(-15,1,17),'XTick',logspace(-3,1,5),...
    'XGrid','on','YGrid','on', 'XMinorGrid','off','YMinorGrid','off',...
    'XMinorTick','off','YMinorTick','off')
% выводим легенду
legend('max|f-s|', '(5/384)h^4')

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


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

Заметим, что существенную роль играет то, что значения первых производных сплайна в граничных точках отрезка интерполирования совпадают со значениями первых производных интерполируемой функции. Если построить сплайн, у которого, к примеру, первая производная в правой точке отрезка будет равна не 0, а 1, то скорость сходимости силно уменьшается и становится порядка h. Для демонстрации этого факта достаточно в m-файле SplineErr.m переменной dfb присвоить значение 1. Получаем тогда следующее поведение ошибки интерполяции кубическим сплайном.


Ошибка интерполяции кубическим сплайном функции sinx при несогласующихся граничных условиях

Гладкость интерполируемой функции так же существенно влияет на скорость убывания ошибки при интерполяции фундаментальными кубическими сплайнами. Рассмотрим, например, интерполяцию функции x2|x| на отрезке [-1;1]. В нуле эта функция имеет разрывную третью производную.

Для интерполяции этой функции фундаментальными сплайнами достаточно в файле SplineErr.m заменить первые пять строчек на следующие

f = inline('x.^2.*abs(x)');
a=-1; 
b=1; 
dfa = -3; 
dfb = 3;
(значения первых производных функции x2|x| в граничных точках отрезка интерполирования равны -3 и 3). В результате получаем график из которого видно, что скорость сходимости ошибки интерполяции имеет порядок h3.


Ошибка интерполяции и ее оценка при интерполяции кубическим фундаментальным сплайном функции x2|x| на отрезке [-1;1].

Дальнейшее снижение гладкости интерполируемой функции приведет к снижению скорости убывания ошибки интерполирования, например при интерполяции функции x|x| (у которой вторая производная разрывна в точке 0) на отрезке [-1;1] кубическим фундаментальным сплайном получим следующий результат.


Ошибка интерполяции и ее оценка при интерполяции кубическим фундаментальным сплайном функции x|x| на отрезке [-1;1].

Аппроксимация производных

Кубические сплайны обеспечивают хорошую аппроксимацию производных интерполируемой функции. Приведем оценки для ошибки аппроксимации производных четырежды непрерывно дифференцируемой функции ƒ(x) при ее интерполяции фундаментальным кубическим сплайном s(x).

Как следует из приведенной ниже оценки, аппроксимация первой производной приводит к ошибке порядка h3.

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

Для демонстрации этого утверждения можно воспользоваться m-файлом DerErr.m, текст которого приведен ниже. В m-файле DerErr.m определяются inline-функции для вычисления интерполируемой функции и ее первой производной. Далее в цикле уменьшается шаг между равномерными узлами интерполяции и для каждого значения шага находится максимальное уклонение производной фундаментального кубического сплайна от производной интерполируемой функции по 30 промежуточным точкам между соседними узлами. Для построения фундаментального кубического сплайна, как и в примере SplineErr.m, используется функция csape, а производная от сплайна находится при помощи функции fnder, входящей в Spline Toolbox, которая возвращает сплайн, составленный из квадратичных полиномов.

% задаем функцию, производную, границы отрезка и значения первых производных в них
f = inline('sin(x)');
df = inline('cos(x)');
a=0; 
b=3*pi/2; 
dfa = df(a); 
dfb = df(b);
%Инициализация массивов для значений шагов и ошибки
H=[];
Err=[];
% построение сплайнов в цикле и вычисление ошибки аппроксимации первой производной
for k=2:12
    % задаем шаг, узлы интерполяции и значения функции в них
    n = pow2(k);  
    h = (b-a)/(n-1);    
    x=a:h:b;
    y=f(x);
    % строим фундаментальный кубический сплайн
    pp=csape(x,[dfa y dfb],[1 1]);
    % вычисляем производную сплайна
    dpp = fnder(pp);
    xx=linspace(a,b,30*n);
    yy=fnval(dpp,xx);
    % находим максимальное уклонение производных по большому набору точек
    err=max(abs(yy-df(xx)));    
    % пополняем массивы со значениями шага и ошибки
    H = [H h];
    Err = [Err err];
end
% создаем графическое окно и оси
figure
axes
% строим зависимость ошибки и оценки от шага в логарифмической шкале
loglog(H, Err, 'LineWidth', 2, 'Color', 'r')
hold on
loglog(H, 1/24*H.^3, 'LineWidth', 1, 'Color', 'b')
% задаем шаги сетки по координатным осям
set(gca, 'YTick',logspace(-15,1,17),'XTick',logspace(-3,1,5),...
    'XGrid','on','YGrid','on', 'XMinorGrid','off','YMinorGrid','off',...
    'XMinorTick','off','YMinorTick','off')
% выводим легенду
legend('max|df-ds|', '(1/24)h^3')

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


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

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

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

    ds=fnval(dpp,x);
    err=max(abs(ds-df(x)));    

Разумеется, для получения графика оценки следует написать

loglog(H, 1/60*H.^4, 'LineWidth', 1, 'Color', 'b')

и соответствующим образом модифицировать легенду.

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


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

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

Для демонстрации этого утверждения достаточно немного изменить m-файл DerErr.m. В начало следует добавить inline-функцию для задания второй производной от интерполируемой функции

d2f = inline('-sin(x)');

Для вычисления второй производной сплайна снова можно воспользоваться fnder, вызвав ее со вторым входным аргументом 2 (порядок производной)

d2pp=fnder(pp,2);

Для нахождения максимального уклонения используем

xx=linspace(a,b,30*n);
yy=fnval(d2pp,xx);
err=max(abs(d2f(xx)-yy));   

Соответственно с оценкой следует изменить ее график

loglog(H, 3/8*H.^2, 'LineWidth', 1, 'Color', 'b')

и легенду

В результате получаем


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

Для третьей производной в случае равноотстоящих узлов интерполирования справедлива следующая оценка

которую легко продемонстрировать, внеся соответствующие изменения в m-файл DerErr.m. Получаем следующий график


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


Поиск по сайту:

Система Orphus

Яндекс.Метрика