最佳逼近方法_大眾科普網

最佳逼近方法

remez算法的實現函數（轉）
%theremezalgorithm(remez1934),alsocalledtheremezexchangealgorithm,
%isanapplicationofthechebyshevalternationtheoremthatconstructs
%thepolynomialofbestapproximationtocertainfunctionsunderanumber
%ofconditions.theremezalgorithmineffectgoesastepbeyondtheminimax
%approximationalgorithmtogiveaslightlyfinersolutiontoanapproximation
%problemconstructsthepolynomialofbestapproximationtocertainfunctions
%--hereis1/(x+2)underanumberofconditions.
%byhuangsheng.
%square
function[result]=remez_com(n)
x=ones(n+2,1);
forj=0:n+1
x(j+1)=cos(pi*(n+1-j)/(n+1));
end
f=1./(x+2);%theconstantofthelinerfunction
one=ones(n+2,1);
two=one;
forj=2:n+2
two(j)=(-1)*two(j-1);
end
a=ones(n+2,n+2);
a(:,1)=one;
a(:,n+2)=two;
forj=2:n+1
fort=2:j
a(:,j)=x.*a(:,j);
%thex.*xisthecoefficientofa2,xisthecoefficientofa1,andoneisfora0andtwoform.
end
end
e=(1e-10)*one;
result=a\f;%result[1]isa0,result[2]isa1,result[3]isa2,andresult[n+2]ism
%thefollowingistosolvetheextremumofthefunction.
while(1)
%thefollowingistoconstructthedifferentialequation.
a_tmp='';
fori=1:n
fort=1:i
if(t==1)
a_tmp=strcat('(',num2str(result(t+1)),')');
elseif(t>1)
x_tmp='t';j=2;
whilej<t
x_tmp=strcat('t*',x_tmp);j=j+1;
end
a_tmp=strcat(num2str(t),'*(',num2str(result(t+1)),')*',x_tmp,'+',a_tmp);
end
end
end
%thefollowingistogettheapproachingfunction.
a_para1=result(1:n+1,[1]);
fori=1:n+1
a_para(i)=a_para1(n+2-i);
end
y_tmp=strcat('-1/(t+2)^2=',a_tmp);
[y]=solve(y_tmp,'t');
y=numeric(y);
%thefollowingistogettingthenewinterleavingsetofpoints.
fori=1:n+1
ify(i)<1&y(i)>-1
forj=2:n+2
ify(i)<x(j)&y(i)>x(j-1)
if(1/(x(j-1)+2)-polyval(a_para,x(j-1)))*(1/(y(i)+2)-polyval(a_para,y(i)))>0
x(j-1)=y(i);
elseif(1/(x(j-1)+2)-polyval(a_para,x(j-1)))*(1/(y(i)+2)-polyval(a_para,y(i)))<0
x(j)=y(i);
end
end
end
end
end
a=ones(n+2,n+2);
a(:,1)=one;
a(:,n+2)=two;
forj=2:n+1
fort=2:j
a(:,j)=x.*a(:,j);
end
end
f=1./(x+2);
result_new=a\f;
if(abs(result_new-result)<e)
break;
end
result=result_new;
end
fplot('1/(x+2)',[-1,1],'k');
holdon;
a_sym=poly2sym(a_para);
ezplot(a_sym,[-1,1]);
num=num2str(n);
legend('1/(x+2)',strcat(num,'次最佳一致逼近'));
title('精度比較');
end