Полиномиальная аппроксимация.

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

q=polyfit(x,y,n)

где y – вектор значений функции; x – вектор значений аргумента; n – порядок аппроксимирующего полинома; p – полученный в результате вектор коэффициентов аппроксимирующего полинома длиной n+1.

Пусть имеется массив значений аргументов

x=[1 2 3 4 4 5 3 6 10 8 9 7];

и массив соответствующих им значений измеряемой величины:

y=[2 2.5 2 4 4.5 5 4 5 9 7 8 7];

Примерим к этим данным функцию polyfit при n=1, n=2 и n=3:

>> x=[1 2 3 4 4 5 3 6 10 8 9 7];

>> y=[2 2.5 2 4 4.5 5 4 5 9 7 8 7];

>> p1=polyfit(x,y,1)

p1 =

0.7918 0.9089

>> p2=polyfit(x,y,2)

p2 =

0.0107 0.6724 1.1593

>> p3=polyfit(x,y,3)

p3 =

-0.0017 0.0392 0.5421 1.3132

В результате получили три вектора p1, p2, и p3, первый из которых содержит коэффициенты аппроксимирующей прямой , второй вектор – коэффициенты аппроксимирующей квадратной параболы , а третий – коэффициенты аппроксимирующей кубической параболы . Аналогичным образом можно аппроксимировать данные параболой четвертой степени, пятой, шестой, седьмой и т.д. На рис. 2.24.1 представлены построенные для наглядности графики исходной зависимости и аппроксимирующих полиномов (для вычисления полиномов использовалась функция polyval):

>> x=[1 2 3 4 4 5 3 6 10 8 9 7];

>> y=[2 2.5 2 4 4.5 5 4 5 9 7 8 7];

>> title('Полиномиальная аппроксимация')

>> xlabel('x'); ylabel('y'); hold on

>> plot(x,y,'ko')

>> p1=polyfit(x,y,1);

>> p2=polyfit(x,y,2);

>> p3=polyfit(x,y,3);

>> t=0:0.05:10;

>> y1=polyval(p1,t);

>> y2=polyval(p2,t);

>> y3=polyval(p3,t);

>> plot(t,y1,'k-',t,y2,'k--',t,y3,'k:'), grid on

>> legend('Данные эксперимента',...

'Аппроксимирующая прямая',...

'Аппроксимирующая квадратная парабола',...

'Аппроксимирующая кубическая парабола',2)

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

Интерполяция сплайнами.

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

Рис. 2.24.1. Пример аппроксимации полиномами.

Рис. 2.24.2. Интерполяция кубическими сплайнами.

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

yy=spline(x,y,xx)

В результате будет выполнена интерполяция значений вектора y, заданного для значений аргумента, представленного вектором x и вернет вектор yy, содержащий значения интерполирующей функции для значений аргумента, заданного вектором xx.

Выполним, например, интерполяцию кубическими сплайнами для значений, заданных векторами x и y, и отобразим полученный результат на рис. 2.24.2.

>> x=[1 2 3 4 5 6 7 8 9 10 11 12];

>> y=[2 3 6 4.5 7 7.5 9 8 10 11 5 1];

>> title('Аппроксимация кубическими сплайнами');

>> xx=0:0.05:12;

>> yy=spline(x,y,xx);

>> xlabel('x'); ylabel('y');

>> plot(x,y,'ko',xx,yy,'k-'), grid on;

>> legend('Экспериментальные данные','Кубический сплайн')

Для одномерной (кусочно-полиномиальной) интерполяции экспериментальных данных в системе MATLAB имеется функция interp1, обращение к которой имеет вид:

yy=interp1(x,y,xx,method)

где x – абсциссы аппроксимируемой функции; y – ординаты аппроксимируемой функции; xx – абсциссы контрольных точек, в которых вычисляются значения аппроксимирующих полиномов, возвращаемые в качестве вектора yy; method – способ аппроксимации, задаваемый в виде строки символов (в действительности, достаточно задать лишь первый символ).

Функция interp1 предоставляет на выбор следующие методы:

'nearest' – ступенчатая интерполяция (здесь значение в каждой промежуточной точке принимается равным ближайшему экспериментальному значению – аппроксимация кусочными полиномами нулевой степени (ступеньками));

'linear' – линейная интерполяция (соседние точки соединяются отрезками прямых в соответствии с данными эксперимента – аппроксимация кусочными полиномами первой степени) – используется по умолчанию;

'spline' – интерполяция кубическими сплайнами;

'pchip' или 'cubic' – интерполяция кусочными полиномами Эрмита третьей степени.

При аппроксимации по методу 'spline' в узловых точках непрерывны первая и вторая производные, а во внутренних узлах, соседних с концевыми также непрерывна и третья производная.

При использовании метода 'pchip' вторая производная может иметь разрывы. Сплайн, построенный по методу 'pchip', является монотонным на любом участке, на котором монотонны данные, этот сплайн воспроизводит локальные экстремумы данных.

Для гладких данных предпочтителен метод 'spline', для негладких целесообразно использовать 'pchip'.

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

Рис. 2.24.3. Сопоставление различных методов интерполяции.

В системе MATLAB можно интерполировать не только одномерные, но также двумерные (для этих целей предназначена функция interp2), трехмерные (для этих целей предназначена функция interp3) и многомерные (для этих целей предназначена функция interpn) данные.

Часть 5.

СИМВОЛЬНЫЕ ВЫЧИСЛЕНИЯ