main.m

用有理多项式做Chebyshev逼近的一种直接方法的matlab程序

clear;
%close all;
format long;
%------------------------------------------------
%in this program y(k)=x(i,k), k=0,1,...l+m+1;
%r(x)=p_{l-m}(x)+(p_{l-m+1}(x)-p_{l-m}(x))/psi(x);
%p_{l-m}(x)=e_0+e_1(x-x0)+...+e_{l-m}(x-x0)...(x-x_{l-m-1});
%p_{l-m+1}(x)=e_0+e_1(x-x0)+...+e_{l-m}(x-x0)...(x-x_{l-m});
%psi(x)=1+(x-x_{l-m+1})_|-e_{l-m+2}+...+(x-x_{l+m+1})_|-e_{l+m};
l=3;
m=3;
n=l+m+1;
MAX=5;
%f(x):C[aa,bb],|f|<ff;   
%g(x):C+[aa,bb],|g|>gg;

ff=log(2);%in here, f=log(x); g=1;
ff=1;%in here, f=1/(1+log(x)); g=1;
gg=1;
aa=1;
bb=2;
%eps=5*10^(-10);
eps1=ff/gg*eps;
%------------------------------------------------
%initial:
x=(cos(pi*[n:-1:0]/n)+1)/2*(bb-aa)+aa,
lb=0;

%-----------------------------------------------

% loops
for run=1:MAX
    run,
    %stage 1----------to calculate a,c,e, lb.
    F=f(x);
    G=g(x);
    H=h(lb,F,G,l,m);
    A=a(x,H,l,m);
    B=b(x,G,l,m);
    if m>0
        [C0,D0]=c0d0(x,A,B,H,lb,G,l,m);
        [C,D,c1,d1]=cd1(x,C0,D0,l,m);
    end
    if m==0
        if abs(A(l+2))>eps
            dlb=-A(l+2)/B(l+2);
            lb=lb+dlb;
            A=A+B*dlb;   
        end
      %  remain=A(l+2),
        E=A(1:l+1);
    else
        for i=1:2
            if abs(C(2*m)-c1)<eps
                break;
            end
            dlb=-(C(2*m)-c1)/(D(2*m)-d1);
            lb=lb+dlb,
            A=A+B*dlb;
            H=h(lb,F,G,l,m);
            [C0,D0]=c0d0(x,A,B,H,lb,G,l,m);
            [C,D,c1,d1]=cd1(x,C0,D0,l,m);
        end
%-----------------------
   %     H=h(0,F,G,l,m);
   %       A=a(x,H,l,m);
    %       lb=21.6*10^(-10);
      %   dlb=21.6*10^(-10),
   %       A=A+B*dlb;
   %       H=h(lb,F,G,l,m);
   %      [C0,D0]=c0d0(x,A,B,H,lb,G,l,m);
   %      [C,D,c1,d1]=cd1(x,C0,D0,l,m),
            %-------------------------------
    %    remain=abs(C(2*m)-c1),
        E=[A,C(2),C(3:2*m)-C(1:2*m-2)];
    end
    
  %  r(x,E,x,l,m),
    
    %stage 2-----------to calculate new x.
    n=length(x);
    h1=diff(x)/5;
    for i=2:n-1
%        [h2(i),index]=min([h1(i-1),h1(i)]);
        
        
      %  v=[x(i)+[-4:0]*h1(i-1),x(i)+[1:4]*h1(i)];%问题：xi越来越小
       v=[x(i)+[0:4]*h1(i)];
        delta1=(r(v,E,x,l,m)-f(v))./g(v);
        if mod(i,2)==0
            [r1(i),k(i)]=max(-sign(lb)*delta1);
         %   kk=k(i),
          %  dd=delta1(k(i)),
        else
            [r1(i),k(i)]=max(sign(lb)*delta1);
         %   kk=k(i),
         %   dd=delta1(k(i)),
        end
            v1(i)=v(k(i));

        if k(i)==1
            v0(i)=v(k(i))-h1(i);
        else
             v0(i)=v(k(i)-1);
        end
        if k(i)==5
             v2(i)=v(k(i))+h1(i);
        else
             v2(i)=v(k(i)+1); 
        end
        
    end
    t0=(r(v0(2:n-1),E,x,l,m)-f(v0(2:n-1)))./g(v0(2:n-1));
    t1=(r(v1(2:n-1),E,x,l,m)-f(v1(2:n-1)))./g(v1(2:n-1));
    t2=(r(v2(2:n-1),E,x,l,m)-f(v2(2:n-1)))./g(v2(2:n-1));
    xx(1)=x(1);
    xx(n)=x(n);
 %   v0,v1,v2

 xx(2:n-1)=v1(2:n-1);
 %xx(2:n-1)=v1(2:n-1)-h1(i)/2.*(t2-t0)./(t0+t2-2*t1);
    delta=(r(xx(2:n-1),E,x,l,m)-f(xx(2:n-1)))./g(xx(2:n-1)),
    if max(abs(abs(lb)-abs(delta)))<=eps1
    %if max(abs(abs(lb)-abs([delta, r([aa bb],E,x,l,m)-f([aa bb])./g([aa bb])])))<=eps1
        break;
    end    
    x=xx,
end

%% e, lb, x have all been calculated.
run,
lb,
E,
x,
tt=[1:.01:2];
%figure
%plot(tt,f(tt),tt,r(tt,E,x,l,m));
%legend('log(x)','the approximation curve')

figure
plot(tt,f(tt)-r(tt,E,x,l,m));
title('error')