Сплайн-аппроксимация и сглаживание данных

Цель работы: Получение навыков проведения сплайн-аппроксимации и сглаживания экспериментальных данных при помощи MathCAD

Используемые материалы: ЭВМ с установленным п/о «MathCAD»

Теоретические сведения

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

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

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

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

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

Линейная интерполяция. Для построения линейной интерполяции служит встроенная функция linterp

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

х - вектор действительных данных аргумента;

у - вектор действительных данных значений того же размера;

t - значение аргумента, при котором вычисляется интерполирующая функция.

Элементы вектора х должны быть определены в порядке возрастания, т. е. Х₁<Х₂<Х₃<. . . <X_N. Например:

Сплайн-интерполяция. Сплайн-интерполяция в MathCAD реализована чуть сложнее линейной. Перед применением функции interp необходимо предварительно определить первый из ее аргументов - векторную переменную s. Делается это при помощи одной из трех встроенных функций тех же аргументов (х,у).

- lspline(x,y) - вектор значений коэффициентов линейного сплайна;

- pspline(x,y) - вектор значений коэффициентов квадратичного сплайна;

- cspline(x,y) - вектор значений коэффициентов кубического сплайна;

- х,у - векторы данных.

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

Полиномиальная сплайн-интерполяция. Более сложный тип интерполяции - так называемая интерполяция В-сплайнами. В отличие от обычной сплайн-интерполяции, сшивка элементарных В-сплайнов производится не в точках x_i; а в других точках u_i, координаты которых предлагается ввести пользователю. Сплайны могут быть полиномами 1, 2 или 3 степени (линейные, квадратичные или кубические). Применяется интерполяция В-сплайнами точно так же, как и обычная сплайн-интерполяция, различие состоит только в определении вспомогательной функции коэффициентов сплайна.

- interp(s,x,y,t) - функция, аппроксимирующая данные векторов х и у с помощью В-сплайнов;

- bspline(x,y,u,n) - вектор значений коэффициентов В-сплайна;

- s - вектор вторых производных, созданный функцией bspline; - х - вектор действительных данных аргумента, элементы которого расположены в порядке возрастания;

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

- u - вектор значений аргумента, в которых производится сшивка В-сплайнов;

- n - порядок полиномов сплайновой интерполяции (l, 2 или 3). Например:

Сглаживание. При анализе данных часто возникает задача их фильтрации, заключающаяся в устранении одной из составляющих зависимости y(x_i). Наиболее часто целью фильтрации является подавление быстрых вариаций y(x_t), которые чаще всего обусловлены шумом. В результате из быстро осциллирующей зависимости y(x_i) получается другая, сглаженная зависимость, в которой доминирует более низкочастотная составляющая.

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

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

- medsmooth(y,b) - сглаживание алгоритмом "бегущих медиан";

- ksmooth(x,y,b) - сглаживание на основе функции Гаусса; - supsmooth(x,y) - локальное сглаживание адаптивным алгоритмом, основанное на анализе ближайших соседей каждой пары данных;

- х - вектор действительных данных аргумента (для supsmooth его элементы должны быть расположены в порядке возрастания); - у - вектор действительных значений того же размера, что и х;

- b - ширина окна сглаживания.

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

Часто бывает полезным совместить сглаживание с последующей интерполяцией или регрессией. Например:

Регрессия. Задачи математической регрессии имеют смысл приближения выборки данных (x_i,y_i) некоторой функцией f(x), определенным образом минимизирующей совокупность ошибок |f(x_i)-Y_i|. Регрессия сводится к подбору неизвестных коэффициентов, определяющих аналитическую зависимость f(х).

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

Полиномиальная регрессия. Полиномиальная регрессия означает приближение данных (x_i,y_i) полиномом k-й степени A(x)=a+b-x+c-x²+d-x³+.. .+h-x^k. При k=1 полином является прямой линией, при k=2 - параболой, при k=3 - кубической параболой и т. д. Как правило, на практике применяются k<5

В MathCAD полиномиальная регрессия осуществляется комбинацией встроенной функции regress и полиномиальной интерполяции.

- regress (х, у, k) - вектор коэффициентов для построения полиномиальной регрессии данных;

- interp(s,x,y,t) - результат полиномиальной регрессии; s:=regress(х,у,k); - х - вектор действительных данных аргумента, элементы которогорасположены в порядке возрастания; - у - вектор действительных данных значений того же размера; - k - степень полинома регрессии (целое положительное число); - t - значение аргумента полинома регрессии.

Для построения графика полиномиальной регрессии после функции regress используется interp. Например:

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

- loess (x, у, span) - вектор коэффициентов для построения регрессии данных отрезками полиномов;

- interp(s,x,y,t) - результат полиномиальной регрессии;

- s:=loess(x,y,span);

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

- у - вектор действительных данных значений того же размера; - span - параметр, определяющий размер отрезков полиномов (положительное число, хорошие результаты дает значение порядка span=0.75).

Параметр span задает степень сглаженности данных. При больших значениях span регрессия практически не отличается от регрессии одним полиномом (например, span=2 дает почти тот же результат, что и приближение точек параболой).

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

В MathCAD можно осуществить регрессию в виде линейной комбинации F(x)=с₁f₁(x)+c₂f₂(x) + ..., где f_i(x) - любые функции пользователя, с_i - подлежащие определению коэффициенты. Такая регрессия осуществляется при помощи функции linfit. Рассмотрим пример, регрессии, приведенный на рис. 1.

Рис. 1. Регрессия линейной комбинацией функций пользователя.

В начале задаются два вектора значений x и у экспериментальной зависимости, которая будет исследоваться. Вообще-то данные должны представлять собой матрицу из одного столбца, но для экономии места на странице данные были написаны в строчку и поставлен оператор транспонирования. Значения в векторе x обязательно должны располагаться либо в порядке возрастания, либо в порядке убывания. Далее располагается функция F(x), которая представляет собой матрицу из одного столбца, каждый элемент которой представляет собой функцию f_i(x), при которой будет найден коэффициент с_i. Затем располагается функция linfit(x,y,F), вычисляющая вектор С - параметров линейной комбинации функций пользователя, осуществляющей регрессию данных. Полученный результат приведен ниже, там же приведена полученная функция f(x) для графической проверки качества регрессии.

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

Рис. 2. Регрессия общего вида.

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

Далее располагается функция G(x,C), которая представляет собой матрицу из одного столбца, первым элементом которой является сама функция пользователя, а остальными – частные производные этой функции по искомым коэффициентам, то есть по C₀ и C₁. (Следует обратить внимание, что искомые коэффициенты C₀ и C₁ являются элементами матрицы, и чтобы ввести матричные индексы используется клавиша с изображением квадратной скобки [.) Таким образом функция G(x,C) имеет размерность N+1. Затем располагается функция genfit(x,y,g,G), вычисляющая вектор параметров С, реализующих регрессию данных с помощью функций пользователя общего вида. Полученный результат приведен ниже, там же приведена полученная функция f(x) для графической проверки качества регрессии.

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

Сведения и приемы, необходимые для выполнения данной лабораторной работы, описаны в /4/.

Порядок выполнения работы

1. При помощи оператора linterp провести интерполяцию функции, заданной табл. 2, и определить её значения в точках 1, 2, 3, 4, 5, 6, 7, 8, 9, 10.

Таблица 2.

Х	У
0,44	0,160
1,51	0,394
2,52	0,776
3,46	1,059
5,43	1,027

Продолжение табл.2

6,54	0,752
7,85	0,501
9,22	0,402
10,42	0,120

2. При помощи кубической, квадратичной и линейной сплайн-аппроксимации (операторы cspline, lspline, pspline) обработать функцию, заданную таблицами 2 и 3. Сравнить результаты и сделать выводы по применимости того или иного вида аппроксимации.

Таблица 3.

Х	У
1	0,5
2	1,0
3	7,5
4	9,0
5	7,0
6	1,0
7	1,5
8	5,0
9	6,0

3. При помощи аппроксимации В-сплайном (оператор bspline) обработать функцию, заданную таблицей 2. Сравнить результаты в зависимости от n – порядка полинома и сделать выводы по применимости.

4. При помощи функций (medsmooth(y,b), ksmooth(x,y,b), supsmooth(x,y)) провести сглаживание данных, представленных в таблице 4. Сделать вывод о наилучшем методе сглаживания для данного набора данных.

Таблица 4.

Х	У
1	1,47
2	3,25
3	6,55
4	7,44
5	8,60
6	7,53
7	7,20
8	4,54
9	4,05
10	4,80
11	5,39
12	5,20
13	7,53
14	7,93
15	9,67
16	12,07
17	13,41
18	14,43
19	16,70
20	17,15
21	17,60
22	16,57
23	16,17
24	15,32
25	12,80
26	11,85
27	10,82

Продолжение табл. 4

28	9,44
29	6,64
30	5,20

5. Провести регрессию кривой из задания 4 при помощи функций regress(X,Y,n) при различных значениях n, и и loess(X,Y,span) при различных значениях span (рекомендуется начать со span=0,75).

6. Провести линейную регрессию общего вида данных, представленных в таблице 5, функцией a0*x+a1*x^3-a2/x при помощи функции linfit

Таблица 5.

Х	У
1	-99,5
2	-48,9
3	-31,6
4	-22,4
5	-16,3
6	-11,5
7	-7,4
8	-3,4
9	0,7
10	5,0
11	9,7
12	14,9
13	20,8
14	27,3
15	34,6
16	42,7

Продолжение табл. 5

17	51,7
18	61,8
19	72,8
20	85,0
21	98,3
22	112,9
23	128,8
24	146,1
25	164,8
26	184,9
27	206,6
28	229,9
29	254,9
30	281,7

7. Провести общую регрессию данных, представленных в таблице 6, функцией вида а0*exp(a1*x)+a2*x^a3 при помощи функции genfit

Таблица 6.

Х	У
1	-7,300
2	-36,356
3	-85,081
4	-153,360
5	-241,037
6	-347,901
7	-473,668
8	-617,954
9	-780,241
10	-959,829

Продолжение табл. 6

11	-1155,775
12	-1366,804
13	-1591,195
14	-1826,627
15	-2069,966
16	-2316,979
17	-2561,956
18	-2797,187
19	-3012,265
20	-3193,142
21	-3320,856
22	-3369,810
23	-3305,451
24	-3081,138
25	-2633,915
26	-1878,796
27	-701,064
28	1054,133
29	3595,824
30	7206,168

Требования к оформлению отчета

Отчет должен содержать распечатку результатов выполнения заданий.

Лабораторная работа №3