Материал: Matlab. Практический подход. Самоучитель

Глава 8 Обработка данных

- На этом и построен мой метод дедукции. - А как это выглядит на практике?

К/ф "Приключения Шерлока Холмса и доктора Ватсона. Знакомство"

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

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

Интерполяционные полиномы

- Нашел!

- Что именно? - Формулу!

К/ф "Приключения Шерлока Холмса и доктора Ватсона. Знакомство"

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

Предположим, что имеется набор узловых точек x1 , x2 , ..., xn и значения y1 , y2 , ..., yn некоторой функции в этих узловых точках. Если речь идет о задаче интерполирования, то она может быть сформулирована так: необходимо построить (в соответствии с определенными критериями) некоторую функцию f (x), такую, что эта функция в узловых точках x1 , x2 , ..., xn принимает значения y1 , y2 , ..., yn , то есть f (xk ) = yk для всех индексов k от 1 до n .

Этого критерия совершенно недостаточно, чтобы однозначно определить интерполяционную функцию f (x). Необходимо уточнить ее вид. Очень часто функцию f (x) ищут в виде полинома. Степень полинома выбирается в соответствии с количеством узловых точек. Если узловых точек n , то степень интерполяционного полинома равняется n −1. Такой полином име-

311

Самоучитель Matlab

ет вид f(x) = a1xn−1 +a2xn−2 +a3xn−3 +... +an−1x +an и содержит n параметров a1 , a2 , ..., an , благодаря которым удается добиться того, чтобы

интерполяционный полином проходил через базовые точки (xk,yk ) (индекс k = 1,2,...,n ). Фактически задача сводится к определению коэффициентов интерполяционного полинома на основании значений базовых точек. Существуют различные алгоритмы решения такой задачи, но в нашем случае все существенно упрощается, поскольку в Matlab есть специальные функции для построения интерполяционного полинома и выполнения сопутствующих вычислений. Тем не менее, кратко рассмотрим математическую суть проблемы.

торую следует решать относительно коэффициентов a1 , a2 , ..., an . Система

то есть для элементов W (индексы i, j = 1,2,...,n ) этой матрицы имеет место соотношение Wij =ij xin−j . Для вычисления матрицы Вандермонда в Matlab предназначена функция vander(). Аргументом функции указыва-

полинома можно вычислить в любой точке – для этого в Matlab тоже есть специальные функции. Рассмотрим документ, представленный в документе на рис. 8.1.

Командами x=[0,0.2,0.4,0.7,0.9,1.2,1.5] и y=cos(2*pi*x)./ (1+x) задаются узловые точки и значения функции в узловых точках соот-

ветственно. В данном случае табулируется функция y(x) = . Ма-

трицу Вандермонда вычисляем командой W=vander(x). Коэффициенты

312

Глава 8. Обработка данных

Рис. 8.1. Вычисление интерполяционного полинома

интерполяционного полинома вычисляются командой a=y*inv(W.'). На этом процесс вычисления интерполяционного полинома собственно и заканчивается. Дальнейший код предназначен для иллюстрации полученного результата. В частности, командами z=0:0.01:1.5 и f=polyval(a,z) создаем набор точек, по которым будем строить полиномиальную кривую. Командой plot(x,y,'ro') на графике отображаем базовые точки, затем с помощью инструкции hold on переходим в режим удержания графики и отображаем интерполяционную кривую с помощью команды plot(z,f,'k-'). Инструкция grid on нужна для отображения координатной сетки. На рис. 8.2 представлен график с базовыми точками и кривая, построенная на основе интерполяционного полинома.

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

Для выполнения полиномиальной интерполяции в Matlab может использоваться функция polyfit().

313

Самоучитель Matlab

Рис. 8.2. Базовые точки и интерполяционная кривая

На заметку

Строго говоря, функция polyfit() предназначена для выполнения аппроксимации. Частным случаем задачи аппроксимации является задача интерполяции.

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

На заметку

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

314

Глава 8. Обработка данных

Пример использования функции polyfit() для создания интерполяционного полинома приведен в документе на рис. 8.3.

Рис. 8.3. Создание интерполяционного полинома с помощью функции polyfit()

Здесь мы использовали такую последовательность команд:

>>x=-5:1:5;

>>y=sin(pi*x./2)./(1+x.^2);

>>P=polyfit(x,y,length(x)-1);

>>z=-5.1:0.01:5.1;

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

>>hold on

>>grid on

>>plot(z,polyval(P,z),'r-')

>>title('Интерполяционный полином')

>>P2=polyfit(x,y,5);

>>plot(z,polyval(P2,z),'k-.')

>>legend('базовые точки','интерполяция','аппроксимация')

Первые две команды (x=-5:1:5 и y=sin(pi*x./2)./(1+x.^2)) создают наборы точек, на основе которых будет строиться интерполяционный полином. Командой P=polyfit(x,y,length(x)-1) создается интерполяционный полином. Собственно, под созданием полинома здесь имеется в виду вычисление коэффициентов этого полинома. Коэффициенты записываются в переменную P. Командой z=-5.1:0.01:5.1 создаем список узловых точек, по которым будет строиться график интерполяционного полинома. Диапазон значений аргумента, для которого отображается график интерполяционного полинома, несколько шире диапазона, определяемого базовыми точками.Базовыеточкиотображаемкомандойplot(x,y,'ks'),акомандой plot(z,polyval(P,z),'r-') строим график интерполяционного полинома. Для сравнения, кроме интерполяционного полинома, выполняется еще

315