БлогNot. Сглаживающий сплайн в MathCAD

Сглаживающий сплайн в MathCAD

В отличие от интерполяционного сплайна, проходящего через известные точки {xi,fi} и уже реализованного в MathCAD с помощью встроенных функций interp и cspline, сглаживающий сплайн через измерения проходить не обязан. Зато он гораздо гибче и способен, при соответствующей настройке, компенсировать "аномальные" измерения или просто "загладить" данные не хуже МНК (точней, лучше - единичный "выброс" fi не "потянет" вверх или вниз весь сглаживающий сплайн, в отличие от полинома МНК со степенным базисом).

Постановка задачи следующая: пусть на сетке x1<x2<...<xN заданы измеренные значения f1, f2, ... , fN некоторой функции f(x). Требуется найти дважды непрерывно дифференцируемую функцию s(x), которую на каждом интервале [xi,xi+1] можно представить в виде

s(x) = c3,i*h3 + c2,i*h2 + c1,i*h + yi ,

где h=x-xi - расстояние до ближайшего слева узла сетки, c3,i , c2,i , c1,i , i = 1, 2, ..., N - 1 - коэффициенты сглаживающего сплайна на i-ом интервале по оси x, а yi - значения сглаживающего сплайна, вычисленные в узлах сетки xi , i = 1, 2, ..., N.

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

  • вещественный массив df длины N, dfi>0, задаёт значения весов в узлах интерполяции, чем больше вес узла, тем сильнее будет заглаживаться функция в его окрестностях. При отсутствии дополнительной априорной информации о восстанавливаемой зависимости можно принять значения dfi одинаковыми, например, равными 1;
  • "глобальный" параметр сглаживания - вещественное значение sm (обозначается также α), sm>0, задаёт общую меру сглаживания данных - с его ростом она увеличивается.

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

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

Функция SmoothSpline(x,f,df,sm) возвращает объект, состоящий из матрицы c коэффициентов сглаживающего сплайна размерностью 3*(N-1) и вектора y значений сглаживающего сплайна в узлах сетки (размерностью N). Код функции здесь не привожу из-за её громоздкости - полторы страницы. Его можно увидеть, скачав приложенный файл, а вызвать функцию можно так:

вызов функции SmoothSpline
вызов функции SmoothSpline

Обратите внимание, что третье измерение "аномально", а остальные представляют собой значения f(xi)=xi2. Соответственно, завышен вес df3 (здесь и далее элементы матриц и векторов нумеруются с единицы, поэтому файл MathCAD начинается определением системной переменной ORIGIN:=1).

Также показан неплохой (для небольших значений N) способ выбора значения параметра сглаживания sm.

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

Функция дополнительно использует только служебный вещественный массив r размерностью 7*(N+2), то есть, достаточно экономична по памяти.

Тест и результат по нему сверен с материалами научно-исследовательского вычислительного центра МГУ, всё совпало :)

Вторая функция - GetSmoothSpline(x,dx,y,c) занимается построением сглаживающего сплайна по известной матрице коэффициентов c и значениям сплайна yi в узлах сетки xi, аргумент x внутри функции меняется от x1 до xlength(x) с шагом dx. Функция вернёт матрицу из двух вектор-столбцов - в первом будут содержаться значения аргумента, во втором - вычисленные значения сплайна. Воспользоваться функцией и затем наглядно вывести результаты сглаживания можно, например, так:

вызов функции GetSmoothSpline
вызов функции GetSmoothSpline

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

Данные на этот раз сгенерируем программно, просто "сняв" их с кусочка синусоиды, затем пару измерений сделаем сильно отличающимися от остальных:

данные с "аномальными" измерениями
данные с "аномальными" измерениями

Вот как справился сглаживающий сплайн с этими двумя геями "выбросами":

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

Видно, что отклонений он почти не заметил.

Ну а хотите ещё лучше - обоснованно выбирайте веса узлов и параметр сглаживания.

Разумеется, сглаживающий сплайн не всесилен - например, если "перемельчить" сетку и переборщить с количеством "аномалий", он может пойти "вразнос", тогда для обработки таких данных придётся привлекать аппарат, способный учесть качественные априорные ограничения ("а вот на этом интервале функция должна возрастать!"). Такой аппарат известен как дескриптивные сплайны, но статей о них я не писал уже давно, может, когда-нибудь :)

 Скачать функции для построения сглаживающего сплайна + простой тест, формат .xmcd (209 Кб)

02.03.2014, 09:27 [15078 просмотров]


теги: числа mathcad

показать комментарии (1)