Построение аппроксимационных полиномов (Лагранжа, Бернштейна, Чебышева)

@Clarisse · Регистрация: 02.11.2012

Author24 — интернет-сервис помощи студентам

Добрый вечер. Есть функция:
(-x(1)+1)^2 + (-x(1)/(1+0.03/x(2))+0.56)^2 + (-x(1)/(1+0.3/x(2))+0.16)^2 + (-x(1)/(1+2/x(2))+0.042)^2 + (-x(1)/(1+5/x(2))+0.005)^2 + (-x(1)/(1+10/x(2)))^2
С минимизацией в матлабе смогла разобраться,а с этим заданием для диплома нужна помощь. Нужно провести аппроксимацию полиномами (Лагранжа, Бернштейна, Чебышева).
Для Лагранжа вроде нашла функции: lagrange и lagrangef, но не совсем понимаю как это доработать. Подскажите,пожалуйста.

@vital792 · 04.03.2014, 17:04

что то странно выглядит ваша функция. Не понятно, одномерная функция это, поверхность или что. В том виде как у вас записано, если х вектор, то это просто число. Я попробовал для интереса аппроксимацию данными методами на примере простой функции sin(5*x). Не совсем уверен в своих манипуляциях с Лагранжем, я просто взял его интерполяционную формулу и использовал для задачи аппроксимации, не знаю, можно ли так делать). Интервал я взял [-1 1], что бы легче было Чебышеву. Но Бернштейну надо [0 1], поэтому сделал замену переменных. Получилось плоховато. Полностью уверен только в коде, реализующем аппроксимацию полиномами Чебышева. Используя их ортогональность, просто аппроксимирую рядом фурье-чебышева. Код:

Matlab M

clear all; close all; clc;
 
N = 40;
x = linspace(-1, 1, N);
y = sin(5*x);
 
% degree
n = 16;
 
%% lagrange
 
coeff = zeros(1, n);
for k=1:n
    num = 1; den = 1;
    for j=1:n
        if j == k; continue; end;
        num = conv(num, [1 x(j)]);
        den = den * (x(k)-x(j));
    end
    coeff = coeff + y(k)*num/den;
end
 
fLagr = coeff(n) * ones(1, N);
for k=1:n-1
    fLagr = fLagr + coeff(n-k)*x.^k;
end
figure; plot(x, y); hold on; plot(x, fLagr, 'r');
 
errLagr = std(y - fLagr)
 
%% bernstain
 
% maps interval 0:1 onto the interval -1:1
a = x(1);
b = x(end);
t = (x-a)/(b-a);
x1 = (0:n)/n;
y1 = sin(5*(a+(b-a)*x1));
 
fBern = 0;
for k=1:n+1
    fBern = fBern + y1(k)*nchoosek(n,k-1) * t.^(k-1).*(1-t).^(n-(k-1));
end
figure; plot(x, y); hold on; plot(x, fBern, 'r');
 
errBern = std(y - fBern)
 
%% chebyshev
 
T(1,:) = ones(size(x));
T(2,:) = x;
coeff = zeros(1, n);
fCheb = -sum(y)*(x(end)-x(1))/length(x)/pi;
for k=1:n
    T(k+2,:) = 2*x.*T(k+1,:) - T(k,:);
    coeff(k) = 2*sum(y.*T(k,:))*(x(end)-x(1))/length(x)/pi;
    fCheb = fCheb + coeff(k)*T(k,:);
end
fCheb = [1 1./sqrt(1-x(2:end-1).^2) 1].* fCheb;
figure; plot(x, y); hold on; plot(x, fCheb, 'r');
 
errCheb = std(y - fCheb)

@Clarisse · 04.03.2014, 21:05 **[ТС]**

Спасибо за пример,думаю, пригодится в дальнейшем. Уточнила задание: надо построить такой полином для заданных значений,а не для функции (например [1, 0.56022, 0.16226, 0.042648, 0.0052242, 0])
Как это можно сделать?

@vital792 · 05.03.2014, 09:36

в данном случае не важно что это, дискретная функция или просто массив данных. Например:

Matlab M

1 2	y = sin(10*(0:0.05:1)); y = y + rand(size(y));

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

Matlab M

1 2	x = linspace(-1, 1, length(y)); n = length(y) - 2;

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

Matlab M

T(1,:) = ones(size(x));
T(2,:) = x;
coeff = zeros(1, n);
fCheb = -sum(y)*(x(end)-x(1))/length(x)/pi;
for k=1:n
    T(k+2,:) = 2*x.*T(k+1,:) - T(k,:);
    coeff(k) = 2*sum(y.*T(k,:))*(x(end)-x(1))/length(x)/pi;
    fCheb = fCheb + coeff(k)*T(k,:);
end
fCheb = [1 1./sqrt(1-x(2:end-1).^2) 1].* fCheb;
 
plot(y); hold on; plot(fCheb, 'r');

Общий вид аппроксимирующей функции в данном случае: F = K0*T0 + K1*T1 + ... + Kn*Tn, где Ti - многочлены Чебышева, а Ki - соответствующие коэффициенты фурье.

@Clarisse · 05.03.2014, 12:34 **[ТС]**

vital792, спасибо! А почему не важна функция? И как мне использовать мои данные( цифры, которые я указала выше)?

@vital792 · 05.03.2014, 13:36

Сообщение от Clarisse

А почему не важна функция?

вид функции никак не используется. Просто замените y = f(x) на y = [ваши данные]; Единственно где используется вид функции, в аппроксимации многочленами Бернштейна, для приведения оси абсцисс к нужному интервалу, для совмещения результатов разных видов аппроксимаций. Но если у вас просто набор данных, вы же можете выбирать любой удобный вам интервал для каждого конкретного случая, как в моем втором примере...

@Clarisse · 05.03.2014, 15:01 **[ТС]**

vital792, а аппроксимация при помощи сплайнов и симплеск- методом можете что- нибудь подсказать?
N это количество узлов или что-то другое? Как выбрать это число?

@vital792 · 05.03.2014, 15:53

что такое симплекс- аппроксимация не слышал. Возможно, тоже самое что и кусочно-линейная или кусочно-полиномиальная? Попробуйте использовать встроенные функции interp1() и spline().

Не по теме:

Вы решили перепробовать все виды аппроксимаций?) Тогда в первую очередь начинать надо с простой полиномиальной при помощи мнк

@Clarisse · 05.03.2014, 16:23 **[ТС]**

vital792, нет, науч рук меняет задания) как выбрать n? В Ваш код для обоих методов можно просто вставить вместо функции значения?

@vital792 · 05.03.2014, 16:28

Сообщение от Clarisse

N это количество узлов или что-то другое?

N это количество значений исходных данных. А n - это в общем случае степень аппроксимирующего многочлена. В аппроксимации многочленами чебышева - это количество членов ряда, но и в этом случае термин "степень" уместно, т.к. старший полином чебышева как раз получается этой же степени

@Clarisse · 05.03.2014, 16:36 **[ТС]**

Раз у меня 6 значений, значит N=6, а чему тогда равно n?

@vital792 · 05.03.2014, 16:37

Сообщение от Clarisse

как выбрать n?

не знаю, наверно произвольно. Чем выше степень, чем меньше ошибка аппроксимации, но соответственно и больше вычислений. Естественно, для каждого метода ошибка будет разной и соотвественно разной будет оптимальная степень, если например в качестве критерия оптимальности задать порог среднеквадратической ошибки.

@Clarisse · 05.03.2014, 22:32 **[ТС]**

А Вы упоминали выше простую полиномиальную аппроксимацию при помощи мнк. Что это такое?

@vital792 · 06.03.2014, 15:39

Сообщение от Clarisse

Что это такое?

это самый очевидный способ аппроксимации, "в лоб". Аппроксимирующий полином выбирается минимизацией среднеквадратической ошибки. Естественно, такой способ будет иметь наименьшую ошибку по сравнению с остальными. Но он довольно затратный и на большом количестве данных довольно медленно работает.

Matlab M

clear all; close all; clc;
 
y = [1, 0.56022, 0.16226, 0.042648, 0.0052242, 0];
N = length(y);
% degree
n = N-2;
 
%% lagrange
 
x = linspace(0, 1, N);
coeff = 0;
for k=1:n
    num = 1; den = 1;
    for j=1:n
        if j == k; continue; end;
        num = conv(num, [1 -x(j)]);
        den = den * (x(k)-x(j));
    end
    coeff = coeff + y(k)*num/den;
end
 
fLagr = polyval(coeff, x);
 
fLagr = (fLagr);
figure; plot(x, y); hold on; plot(x, (fLagr), 'r');
 
errLagr = std(y - fLagr)
 
%% bernstain
 
x = linspace(0, 1, N);
fBern = 0;
for k=1:n
    fBern = fBern + y(k)*nchoosek(n,k-1) * x.^(k-1).*(1-x).^(n-(k-1));
end
figure; plot(x, y); hold on; plot(x, fBern, 'r');
 
errBern = std(y - fBern)
 
%% chebyshev
 
x = linspace(-1, 1, N);
% T(1,:) = ones(size(x));
% T(2,:) = x;
coeff = 0;
fCheb = -sum(y)*(x(end)-x(1))/length(x)/pi;
for k=1:n
%     T(k+2,:) = 2*x.*T(k+1,:) - T(k,:);
    T = cos((k-1)*acos(x));
    coeff(k) = 2*sum(y.*T)*(x(end)-x(1))/length(x)/pi;
    fCheb = fCheb + coeff(k)*T;
end
fCheb = [1 1./sqrt(1-x(2:end-1).^2) 1].* fCheb;
figure; plot(x, y); hold on; plot(x, fCheb, 'r');
 
errCheb = std(y - fCheb)
 
%% lms
x = 0:length(y)-1;
p = polyfit(x, y, n-1);
fPoly = polyval(p, x);
figure; plot(x, y); hold on; plot(x, fPoly, 'r');
 
errPoly = std(y - fPoly)

Кстати в моем первом коде для метода лагранжа у меня была ошибка - забыл знак "-". Здесь исправил.

@Clarisse · 06.03.2014, 16:35 **[ТС]**

vital792, спасибо Вам огромное! Еще один вопрос: почему Вы выбираете n=N-2?

@vital792 · 06.03.2014, 17:23

Сообщение от Clarisse

почему Вы выбираете n=N-2?

да в принципе наугад. Оптимально конечно выбрать n=N, но не показательно. Вообще 6 точек конечно маловато, например по вашим данным можно сказать, что аппроксимация полиномами чебышева плохая, а мнк хорошая. Это не всегда так, и ниже пример демонстрации плохих свойств всех выбранных методов. Я взял последовательность похожую на вашу, только побольше точек:

Matlab M

clear all; close all; clc;
 
% y = [1, 0.56022, 0.16226, 0.042648, 0.0052242, 0];
y = 2-exp(0:.01:1);
y = y+0.1*rand(size(y));
N = length(y);
% degree
n = N-15;
 
%% lagrange
 
x = linspace(0, 1, N);
coeff = 0;
for k=1:n
    num = 1; den = 1;
    for j=1:n
        if j == k; continue; end;
        num = conv(num, [1 -x(j)]);
        den = den * (x(k)-x(j));
    end
    coeff = coeff + y(k)*num/den;
end
 
fLagr = polyval(coeff, x);
 
fLagr = (fLagr);
figure; plot(x, y); hold on; plot(x, (fLagr), 'r');
 
errLagr = std(y - fLagr)
 
%% bernstain
 
x = linspace(0, 1, N);
fBern = 0;
for k=1:n
    fBern = fBern + y(k)*nchoosek(n,k-1) * x.^(k-1).*(1-x).^(n-(k-1));
end
figure; plot(x, y); hold on; plot(x, fBern, 'r');
 
errBern = std(y - fBern)
 
%% chebyshev
 
x = linspace(-1, 1, N);
% T(1,:) = ones(size(x));
% T(2,:) = x;
coeff = 0;
fCheb = -sum(y)*(x(end)-x(1))/length(x)/pi;
for k=1:n
%     T(k+2,:) = 2*x.*T(k+1,:) - T(k,:);
    T = cos((k-1)*acos(x));
    coeff(k) = 2*sum(y.*T)*(x(end)-x(1))/length(x)/pi;
    fCheb = fCheb + coeff(k)*T;
end
fCheb = [1 1./sqrt(1-x(2:end-1).^2) 1].* fCheb;
figure; plot(x, y); hold on; plot(x, fCheb, 'r');
 
errCheb = std(y - fCheb)
 
%%
x = 1:length(y);
p = polyfit(x, y, n-1);
fPoly = polyval(p, x);
figure; plot(x, y); hold on; plot(x, fPoly, 'r');
 
errPoly = std(y - fPoly)

Лагранж в моей кривой реализации быстро переполнился и начал показывать погоду. Бернштейну нужны биномиальные коэффициенты. Для больших чисел C(n,k) матлаб тоже начал ругаться, но в целом метод отработал неплохо. LMS начал жаловаться на плохую обусловленность и тоже не дал результата. Чебышев только ошибся на границах, что нормально для ряда фурье. (в цифровой обработке сигналов это известно как явление Гиббса). И это еще без учета скорости работы методов

@Clarisse · 06.03.2014, 17:37 **[ТС]**

Вообще я немного ошиблась, должно быть 8 точек, исправлю, спасибо Вам огромное за помощь!

Построение аппроксимационных полиномов (Лагранжа, Бернштейна, Чебышева)

Решение

Решение

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
		1
	Построение аппроксимационных полиномов (Лагранжа, Бернштейна, Чебышева) 03.03.2014, 22:58. Показов 7180. Ответов 16 Метки нет (Все метки) Добрый вечер. Есть функция: (-x(1)+1)^2 + (-x(1)/(1+0.03/x(2))+0.56)^2 + (-x(1)/(1+0.3/x(2))+0.16)^2 + (-x(1)/(1+2/x(2))+0.042)^2 + (-x(1)/(1+5/x(2))+0.005)^2 + (-x(1)/(1+10/x(2)))^2 С минимизацией в матлабе смогла разобраться,а с этим заданием для диплома нужна помощь. Нужно провести аппроксимацию полиномами (Лагранжа, Бернштейна, Чебышева). Для Лагранжа вроде нашла функции: lagrange и lagrangef, но не совсем понимаю как это доработать. Подскажите,пожалуйста. 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	04.03.2014, 21:05 [ТС]	3
	Спасибо за пример,думаю, пригодится в дальнейшем. Уточнила задание: надо построить такой полином для заданных значений,а не для функции (например [1, 0.56022, 0.16226, 0.042648, 0.0052242, 0]) Как это можно сделать? 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	05.03.2014, 12:34 [ТС]	5
	vital792, спасибо! А почему не важна функция? И как мне использовать мои данные( цифры, которые я указала выше)? 0

@vital792 2014 / 1286 / 61 Регистрация: 05.06.2010 Сообщений: 2,213
	05.03.2014, 13:36	6
	Сообщение от Clarisse А почему не важна функция? вид функции никак не используется. Просто замените y = f(x) на y = [ваши данные]; Единственно где используется вид функции, в аппроксимации многочленами Бернштейна, для приведения оси абсцисс к нужному интервалу, для совмещения результатов разных видов аппроксимаций. Но если у вас просто набор данных, вы же можете выбирать любой удобный вам интервал для каждого конкретного случая, как в моем втором примере... 1

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	05.03.2014, 15:01 [ТС]	7
	vital792, а аппроксимация при помощи сплайнов и симплеск- методом можете что- нибудь подсказать? N это количество узлов или что-то другое? Как выбрать это число? 0

@vital792 2014 / 1286 / 61 Регистрация: 05.06.2010 Сообщений: 2,213
	05.03.2014, 15:53	8
	что такое симплекс- аппроксимация не слышал. Возможно, тоже самое что и кусочно-линейная или кусочно-полиномиальная? Попробуйте использовать встроенные функции interp1() и spline(). Не по теме: Вы решили перепробовать все виды аппроксимаций?) Тогда в первую очередь начинать надо с простой полиномиальной при помощи мнк 1

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	05.03.2014, 16:23 [ТС]	9
	vital792, нет, науч рук меняет задания) как выбрать n? В Ваш код для обоих методов можно просто вставить вместо функции значения? 0

@vital792 2014 / 1286 / 61 Регистрация: 05.06.2010 Сообщений: 2,213
	05.03.2014, 16:28	10
	Сообщение от Clarisse N это количество узлов или что-то другое? N это количество значений исходных данных. А n - это в общем случае степень аппроксимирующего многочлена. В аппроксимации многочленами чебышева - это количество членов ряда, но и в этом случае термин "степень" уместно, т.к. старший полином чебышева как раз получается этой же степени 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	05.03.2014, 16:36 [ТС]	11
	Раз у меня 6 значений, значит N=6, а чему тогда равно n? 0

@vital792 2014 / 1286 / 61 Регистрация: 05.06.2010 Сообщений: 2,213
	05.03.2014, 16:37	12
	Сообщение от Clarisse как выбрать n? не знаю, наверно произвольно. Чем выше степень, чем меньше ошибка аппроксимации, но соответственно и больше вычислений. Естественно, для каждого метода ошибка будет разной и соотвественно разной будет оптимальная степень, если например в качестве критерия оптимальности задать порог среднеквадратической ошибки. 1

	06.03.2014, 17:37

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	05.03.2014, 22:32 [ТС]	13
	А Вы упоминали выше простую полиномиальную аппроксимацию при помощи мнк. Что это такое? 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	06.03.2014, 16:35 [ТС]	15
	vital792, спасибо Вам огромное! Еще один вопрос: почему Вы выбираете n=N-2? 0

@Clarisse 1 / 1 / 0 Регистрация: 02.11.2012 Сообщений: 63
	06.03.2014, 17:37 [ТС]	17
	Вообще я немного ошиблась, должно быть 8 точек, исправлю, спасибо Вам огромное за помощь! 0