Как сделать интерполяцию в матлабе

Добавил пользователь Евгений Кузнецов
Обновлено: 04.10.2024

Двумерная интерполяция существенно сложнее, чем одномерная, рассмотренная выше, хотя смысл ее тот же – найти промежуточные точки некоторой зависимости z(x, у) вблизи расположенных в пространстве узловых точек. Для двумерной табличной интерполяции используется функция interp2:

  • ZI = interp2(X,Y.Z,XI.YI) – возвращает матрицу ZI, содержащую значения функции в точках, заданных аргументами XI и YI, полученные путем интерполяции двумерной зависимости, заданной матрицами X, Y и Z. При этом X и Y должны быть монотонными и иметь тот же формат, как если бы они были получены с помощью функции meshgrid (строки матрицы X являются идентичными; то же можно сказать о столбцах массива Y). Матрицы X и Y определяют точки, в которых задано значение Z. Параметры XI и YI могут быть матрицами, в этом случае interp2 возвращает значения Z, соответствующие точкам (XI(i,j),YI(i.j)). В качестве альтернативы можно передать в качестве параметров вектор-строку xi и вектор-столбец yi. В этом случае interp2 представляет эти векторы так, как если бы использовалась команда mesh-grid(xi.yi);
  • ZI = interp2(Z,XJ.YI) – подразумевает, что Х=1:n и Y=l:m, где [m.n]=size(Z);
  • ZI = interp2(Z,ntimes) – осуществляет интерполяцию рекурсивным методом с числом шагов ntimes;
  • ZI = interp2(X,Y,Z.XI,YI.method) – позволяет с помощью опции method задать метод интерполяции:
    • 'nearest' – интерполяция по соседним точкам;
    • 'linear' – линейная интерполяция;
    • 'cubic' – кубическая интерполяция (полиномами Эрмита);
    • 'spline' – интерполяция сплайнами.

    Все методы интерполяции требуют, чтобы X и Y изменялись монотонно и имели такой же формат, как если бы они были получены с помощью функции meshgrid. Когда X и Y – векторы равномерно распределенных точек, для более быстрой интерполяции лучше использовать методы '*1inear', '*cubic', или '*nearest'.


    Рис. 17.13. Применение функции interpZ

    Рисунок 17.13 иллюстрирует применение функции interp2 для двумерной интерполяции (на примере функции peaks).

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

    Не будем вдаваться в математические определения и термины, перейдём сразу к сути:


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

    В Matlab такой способ реализован с помощью команды
    interp1(x,y, xi, 'linear') или просто interp1(x,y, xi), где xи y массивы из табличных данных (координаты точек), xi — массив промежуточных точек, координаты которых требуется найти.

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

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

    Разберём для самого часто встречающегося полинома реализацию в Matlab. Вот пример использования:

    Проинтерполировать функцию sin x на отрезке [1, 9] с шагом 2 и построим графики sin x и полученного интерполяционного полинома.

    function yy=lagrange(x,y,xx)% вычисление интерполяционного полинома в форме Лагранжа% x - массив координат узлов% y - массив значений интерполируемой функции% xx - массив значений аргумента, для которых надо вычислить значения полинома% yy - массив значений полинома в точках xx % узнаем число узлов интерполяции (N=n+1)N=length(x);% создаем нулевой массив значений интерполяционного полиномаyy=zeros(size(xx));% в цикле считаем сумму по узламfor k=1:N % вычисляем произведения, т.е. функции Psi_k t=ones(size(xx)); for j=[1:k-1, k+1:N] t=t.*(xx-x(j))/(x(k)-x(j)); end % накапливаем сумму yy = yy + y(k)*t;end

    Теперь создайте ещё один файл и запишем в него само решение поставленной задачи:

    % задание узлов интерполяцииx=1:2:9;y=sin(x);% задание точек, в которых требуется найти значения интерполяционного полиномаxx=linspace(1,9,1000);% нахождение значений интерполяционного полиномаyy=lagrange(x,y,xx);% построение графиковfigure('Color','w')% вывод графика sin(x)fplot(@sin,[1 9])hold on% вывод графика полиномаplot(xx,yy,'r')% вывод узлов интерполяцииplot(x,y,'bo')% размещение легендыlegend('sin\itx','<\itL_n (интерполяция)>','узлы',-1)

    Ссылки на файлы — исходники сможете найти в конце статьи. Более подробную информацию о полиноме Лагранжа вы сможете найти на официальном сайте Matlab.

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

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

    Вычисление сплайн интерполяции в Matlab осуществляется с помощью команды spline(x, y, xx), где х и у — массивы табличных данных, а хх — промежуточные значения по оси абцисс (Х). Вот небольшой пример:

    Рассмотрим построение полинома n -й степени в Matlab на примере следующей задачи.

    Задача 1 . Найти приближенное значение функции при заданных значениях аргумента с помощью интерполяционных полиномов n -й степени в точках x 1=0,702, x 2=0,512, x 3=0,608. Функция задана таблично с постоянным шагом:

    В Matlab отсутствуют встроенные функции для построения интерполяционного полинома n -й степени, однако средств пакета достаточно для написаний функций, реализующих полином. Ниже представлена функция на М-языке для вычисления в точке t по экспериментальным значениям ( xi,yi,i=1,2,…,n) значения:

    1) Вычисление значения канонического полинома

    Function s=kanon(x, y, t)

    % вычисление количества точек в массивах x и y

    % формирование коэффициентов системы уравнений :

    % Решение системы уравнений :

    % вычисление значения канонического полинома по формуле :

    2) Вычисление значений полинома Ньютона

    Function s = newton(x,y,t)

    % Вычисление количества точек в массивах x и y для полинома Ньютона:

    % запись в первый столбец матрицы разделенных разностей вектора y :

    % формирование матрицы разделенных разностей :

    % формирование массива коэффициентов полинома Ньютона :

    % расчет значения полинома Ньютона в точке t по формуле :

    3) Вычисление значения полинома Лагранжа

    Function s=lagrang (x, y, t)

    % вычисление количества точек в массивах x и y

    % расчет суммы произведений по формуле

    % для вычисления значения полинома Лагранжа в точке t

    Для вычисления значения функции в точках x 1=0,702, x 2=0,512, x 3=0,608 необходимо сформировать массив r =(0,702*0,512*0,608) и в цикле обратиться к одной из приведенных в программах функций. Далее в программе показано обращение к функции lagrang (так же можно обратиться к функциям kanon и newton ) для вычисления значения полинома в точках x 1=0,702, x 2=0,512, x 3=0,608.

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


    МНК (Метод Наименьших Квадратов)

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

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

    Аппроксимация в Matlab по МНК осуществляется с помощью функции polyfit. Функция p = polyfit(x, y, n) находит коэффициенты полинома p(x) степени n, который аппроксимирует функцию y(x) в смысле метода наименьших квадратов. Выходом является строка p длины n+1, содержащая коэффициенты аппроксимирующего полинома.

    Пример (аппроксимация с помощью встроенных MATLAB функций):

    Осуществить аппроксимацию в Matlab табличных данных x = [0, 0.1 , 0.2, 0.3, 0.5] и y = [3, 4.5, 1.7, 0.7, -1]. Применяя метод наименьших квадратов, приблизить ее многочленами 1-ой и 2-ой степени. Для каждого определить величину среднеквадратической ошибки. Построить (на одном листе) графики и заданной таблично функции (ломанная линия) и приближающих ее многочленов 1-ой и 2-ой степени.


    Как построить график по n точкам? Самое простое — отметить их маркерами на координатной сетке. Однако для наглядности их хочется соединить, чтобы получить легко читаемую линию. Соединять точки проще всего отрезками прямых. Но график-ломаная читается довольно тяжело: взгляд цепляется за углы, а не скользит вдоль линии. Да и выглядят изломы не очень красиво. Получается, что кроме ломаных нужно уметь строить и кривые. Однако тут нужно быть осторожным, чтобы не получилось вот такого:

    Немного матчасти


    Восстановление промежуточных значений функции, которая в данном случае задана таблично в виде точек P1 .  Pn, называется интерполяцией. Есть множество способов интерполяции, но все они могут быть сведены к тому, что надо найти n – 1 функцию для расчёта промежуточных точек на соответствующих сегментах. При этом заданные точки обязательно должны быть вычислимы через соответствующие функции. На основе этого и может быть построен график:

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

    Ставим опыты

    Самый простой пример — линейная интерполяция, в которой используются полиномы первой степени, а в итоге получается ломаная, соединяющая заданные точки.
    Давайте добавим немного конкретики. Вот набор точек (взяты почти с потолка):


    Результат линейной интерполяции этих точек выглядит так:

    Однако, как отмечалось выше, иногда хочется получить в итоге гладкую кривую.


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

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


    Другое традиционное решение, кроме кубических сплайнов дефекта 1 — полиномы Лагранжа. Это полиномы степени n – 1, принимающие заданные значения в заданных точках. То есть членения на сегменты здесь не происходит, вся последовательность описывается одним полиномом.
    Но вот что получается:

    Гладкость, конечно, присутствует, но наглядность пострадала так сильно, что… пожалуй, стоит поискать другие методы. На некоторых наборах данных результат выходит нормальный, но в общем случае ошибка относительно линейной интерполяции (и, соответственно, ложные экстремумы) может получаться слишком большой — из-за того, что тут всего один полином на все сегменты.



    Как это можно использовать для интерполяции? На основе этих кривых тоже можно построить сплайн. То есть на каждом сегменте сплайна будет своя кривая Безье k-й степени (кстати, k = 1 даёт линейную интерполяцию). И вопрос только в том, какое k взять и как найти k – 1 промежуточную точку.
    Здесь бесконечно много вариантов (поскольку k ничем не ограничено), однако мы рассмотрим классический: k = 3.
    Чтобы итоговая кривая была гладкой, нужно добиться дефекта 1 для составляемого сплайна, то есть сохранения непрерывности первой и второй производных в точках сочленения сегментов (Pi), как это делается в классическом варианте кубического сплайна.
    Решение этой задачи подробно (с исходным кодом) рассмотрено здесь.
    Вот что получится на нашем тестовом наборе:

    Стало лучше: ложные экстремумы всё ещё есть, но хотя бы не так сильно отличаются от реальных.

    Думаем и экспериментируем


    Можно попробовать ослабить условие гладкости: потребовать дефект 2, а не 1, то есть сохранить непрерывность одной только первой производной.
    Достаточное условие достижения дефекта 2 в том, что промежуточные контрольные точки кубической кривой Безье, смежные с заданной точкой интерполируемой последовательности, лежат с этой точкой на одной прямой и на одинаковом расстоянии:

    В качестве прямых, на которых лежат точки Ci – 1 (2) , Pi и Ci (1) , целесообразно взять касательные к графику интерполируемой функции в точках Pi. Это гарантирует отсутствие ложных экстремумов, так как кривая Безье оказывается ограниченной ломаной, построенной на её контрольных точках (если эта ломаная не имеет самопересечений).

    Методом проб и ошибок эвристика для расчёта расстояния от точки интерполируемой последовательности до промежуточной контрольной получилась такой:

    image


    Первая и последняя промежуточные контрольные точки равны первой и последней точке графика соответственно (точки C1 (1) и Cn – 1 (2) совпадают с точками P1 и Pn соответственно).
    В этом случае получается вот такая кривая:

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


    К текущему варианту мы пришли, уменьшив гладкость на один порядок. Можно сделать это ещё раз: пусть сплайн будет иметь дефект 3. По факту, тем самым формально функция не будет гладкой вообще: даже первая производная может терпеть разрывы. Но если рвать её аккуратно, визуально ничего страшного не произойдёт.
    Отказываемся от требования равенства расстояний от точки Pi до точек Ci – 1 (2) и Ci (1) , но при этом сохраняем их все лежащими на одной прямой:

    Эвристика для вычисления расстояний будет такой:


    Результат получается такой:

    Если абсцисса точки пересечения касательных в точках Pi(xi, yi) и Pi + 1(xi + 1, yi + 1) лежит в отрезке [xi; xi + 1], то l1 либо l2 полагаем равным нулю. В том случае, если касательная в точке Pi направлена вверх, нулю полагаем максимальное из l1 и l2, если вниз — минимальное.


    Результат следующий:

    На этом было принято решение признать цель достигнутой.
    Может быть, кому-то пригодится код.

    А как люди-то делают?

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

    MS Excel


    Это очень похоже на рассмотренный выше сплайн дефекта 1, основанный на кривых Безье. Правда, в отличие от него в чистом виде, тут всего два ложных экстремума — первый и второй сегменты (у нас было четыре). Видимо, к классическому поиску промежуточных контрольных точек тут добавляются ещё какие-то эвристики. Но ото всех ложных экстремумов они не спасли.

    LibreOffice Calc


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


    Есть там ещё один тип интерполяции, который мы тут не рассматривали: B-сплайн. Но для нашей задачи он явно не подходит, так как даёт вот такой результат :)

    Highcharts, одна из самых популярных JS-библиотек для построения диаграмм


    amCharts, ещё одна популярная JS-библиотека


    Картина очень похожа на экселевскую, те же два ложных экстремума в тех же местах.

    Coreplot, самая популярная библиотека построения графиков для iOS и OS X


    Есть ложные экстремумы и видно, что используется сплайн дефекта 1 на основе Безье.
    Библиотека открытая, так что можно посмотреть в код и убедиться в этом.

    aChartEngine, вроде как самая популярная библиотека построения графиков для Android


    Вместо заключения

    Читайте также: