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

симирующей зависимости в случае линейной функции f (x,a1,a2,...,am ) для определения оптимальных параметров a1 , a2 , ..., am достаточно методов линейной алгебры. В общем случае поиск параметров a1 , a2 , ..., am аппроксимирующей зависимости формулируется как задача поиска абсолютного

n

минимума для выражения ∑(yk − f(xk,a1,a2,...,am ))2 .

k =1

На заметку

В линейном случае задача решается проще и быстрее. В этом отношении лучше и приятнее иметь дело с линейной задачей. К счастью, нередко нелинейную задачу можно переформулировать так, что она станет линейной. Для линейной задачи на будущее договоримся называть параметр ak (k = 1,2,...,m ) коэффициентами аппроксимирующего выражения, а функции ϕk (x) (k = 1,2,...,m ) базисными функциями аппроксимирующего выражения.

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

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

>>x=0:0.01:10;

>>plot(k,y,'rs','LineWidth',2)

>>grid on

>>hold on

>>q=polyfit(k,y,2)

q =

>> plot(x,polyval(p,x),'b-','LineWidth',2)

321

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

Рис. 8.7. Использование функции polyfit() для выполнения полиномиальной аппроксимации

>>plot(x,polyval(q,x),'g--','LineWidth',2)

>>legend('базовые точки','полином 4-й степени','полином 2-й степени')

Узловые точки по аргументу создаем командой k=0:10 (целочисленные значения от 0 до 10 включительно). Значения табулированной функции в этих точках вычисляются командой y=cos(pi*k/5)./(1+k). Таким об-

значений индекса k = 0,1,...,10 . Аппроксимирующий полином вычисляется командой p=polyfit(k,y,4). Первым аргументом функции polyfit() передается список k со значениями узловых точек по аргументу. Список y содержит значения табулированной функции в узловых точках. Третий аргумент функции polyfit() определяет степень аппроксимирующего полинома (в данном случае создается полином четвертой степени). В результате возвращается список коэффициентов аппроксимирующего полинома. Этот список записывается в переменную p.

Базовые точки отображаются на графике с помощью команды plot(k,y,'rs','LineWidth',2). Командой q=polyfit(k,y,2) создается еще один аппроксимирующий полином, но на этот раз второй степени. Затем командами plot(x,polyval(p,x),'b-','LineWidth',2)

иplot(x,polyval(q,x),'g--','LineWidth',2) отображаются

322

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

аппроксимирующие полиномы. При этом переменная x предварительно определяется командой x=0:0.01:10. Для вычисления значения полинома использовалась функция polyval(). Легенда добавляется командой legend('базовые точки','полином 4-й степени','полином 2-й степени'). Результаты аппроксимации лучше всего исследовать в графическом виде. На рис. 8.8 представлены результаты построений, выполненные при создании аппроксимирующих полиномов.

Рис. 8.8. Результат аппроксимации полиномами разной степени

Относительно "качества" аппроксимации, то совершенно очевидно, что полином четвертой степени в этом отношении более чем хорош. Но такая идиллия имеет место далеко не всегда. Другими словами, полином не всегда является приемлемым типом аппроксимирующей функции. Составим простенький программный код, с помощью которого будем вычислять аппроксимирующие зависимости линейного вида для произвольных базисных функций. Мы будем использовать три функции. Вот их код:

function res=LSFit(x,y,phi,m) b=zeros(m);

c=zeros(m,1); for i=1:m

for j=1:m b(i,j)=sum(phi(x,i).*phi(x,j));

end

323

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

c(i)=sum(y.*phi(x,i)); end

res=b\c; end

Это функция для вычисления коэффициентов возле базисных функций в аппроксимирующем выражении. У функции LSFit() четыре аргумента. Первые два – это список узловых точек x и список y значений табулированной функции. По этим точкам строится аппроксимирующая зависимость (точнее, определяются параметры этой зависимости). Указатель phi на базисные функции, через линейную комбинацию которых определяется аппроксимирующее выражение, является третьим аргументом функции LSFit(). Четвертый аргумент m функции LSFit() определяет количество слагаемых в аппроксимирующем выражении. В качестве результата функция возвращает список коэффициентов, с которыми базисные функции входят в аппроксимирующее выражение (переменная-результат обозначена как res).

В теле функции командами b=zeros(m) и c=zeros(m,1) создается квадратная матрица b размера m, заполненная нулями, а также нулевой векторстолбец c. Матрица и вектор заполняются в процессе выполнения вложенных операторов цикла. В каждом цикле индексная переменная пробегает значения от 1 до m. Командой b(i,j)=sum(phi(x,i).*phi(x,j)) вычисляются элементы матрицы b. Элементы вектора c вычисляются командой c(i)=sum(y.*phi(x,i)).

На заметку

n

В переменную b записывается матрица с элементами bij = ∑ϕi(xk )ϕj (xk ), k =1

где базисные функции ϕp(x) ( p = 1,2,...,m ) формируют аппроксимирую-

щее выражение f (x) = ϕ1(x)a1 + ϕ2(x)a2 +... + ϕm(x)am . При фиксированных индексах i и j для вычисления суммы по индексу k в команде

b(i,j)=sum(phi(x,i).*phi(x,j)) использована встроенная функция sum(). Поскольку в инструкции phi(x,i).*phi(x,j) использован оператор .* поэлементного произведения, то результатом этой инструкции является список, который получается поэлементным произведением списков phi(x,i) и phi(x,j). При этом мы предполагаем, что функция, на которую указывает указатель phi, такая, что ее первым аргументом может быть список. Как раз для вычисления суммы по всем элементам списка и применялась функция sum(). Аналогичным образом вычисляется результат выражения

c(i)=sum(y.*phi(x,i)). В этом случае вычисляются элементы вектора n

ci = ∑ykϕi(xk )

324

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

Результат (вектор-столбец коэффициентов аппроксимирующего выражения) вычисляется командой res=b\c. Окно редактора m-файлов с кодом этой функции представлено в документе на рис. 8.9.

Рис. 8.9. Код функции для вычисления коэффициентов аппроксимирующего выражения

Если функция LSFit() позволяет вычислить коэффициенты в аппроксимирующем выражении, то функция GetFit() предназначена для вычисления по этим коэффициентам значения аппроксимирующего выражения. Она имеет следующий код:

function res = GetFit(a,x,phi) m=length(a);

res=0; for i=1:m

res=res+phi(x,i)*a(i); end

end

У функции три аргумента: список коэффициентов в линейной комбинации (переменная a), список со значениями аргумента, для которых вычисляется значение аппроксимирующего выражения (переменная x), а также указатель на базовую функцию (переменная phi). Функция GetFit() возвращает результат – список значений аппроксимирующего выражения в узловых точках, определяемых списком x (второй аргумент функции). Переменная-результат обозначена как res. В теле функции командой m=length(a) определяется количество элементов в списке коэффициентов. Параметр m автоматически определяет количество слагаемых в линейной комбинации базисных функций для аппроксимирующего выражения.

325