%Es werden Lagrange-Koeffizienten, Lebesgue-Funktionen und -Konstanten 
%für zwei Haarsche Räume V, (i) und (ii) 
%und zwei Knotenverteilungen (1: äquidistant, 2: Tschebyscheff) berechnet.
%Programm: Gerhard Opfer, 4/1999
%(i)   V=<1,t,...,t^{n-1}>,
%(ii)  V=<1/(t-tau_1),1/(t-tau_2),...,1/(t-tau_n),
% t in [li,re],
%Knoten=(t_1,t_2,...,t_n) in [li,re].
%tau=(t_n+1,t_1-1,t_n+2,t_1-2,...), außerhalb [li,re]
%Zuerst Raum (i)
clear
n=4;                                             %Wahl der Dimension
t=linspace(-6,6,1000); li=t(1); re=t(length(t)); %Wahl des Intervalls
for Fall=1:2 
  if Fall==1   %Fall 1: äquidistante Knoten
    if n==1, 
      Knoten=0.5*(li+re); else
      Knoten=linspace(li,re,n); %äquidistant
    end; %if n
  end; %if Fall
  if Fall==2 %Fall 2: Tschebyscheffknoten
    j=1:2:2*n-1;
    Knoten=0.5*(li+re+(re-li)*cos(j/(2*n)*pi)); 
         %Tschebyscheffknoten von rechts nach links
    Knoten=fliplr(Knoten); %in aufsteigender Reihenfolge
  end; %if Fall
  for k=1:n
    Knotenk=Knoten; Knotenk(k)=[];
    for j=1:length(t)
      L(j,k)=prod(t(j)-Knotenk)/prod(Knoten(k)-Knotenk);
    end; %for j
    figure(1)
    subplot(4,n,k+2*n*(Fall-1));
    plot(t,L(:,k),[li,re],[0,0],Knoten,zeros(size(Knoten)),'o',...
    [Knoten(k),Knoten(k)],[0,1]);
    axis([li,re,0,2]); %ggf. weglassen oder anpassen
  end; %for k
  Lebesgue=sum(abs(L')); %Lebesgue-Funktion
  lambda=max(Lebesgue);  %Lebesgue-Konstante
  [n,lambda] 
  figure(2)
  subplot(1,4,Fall);
  plot(t,Lebesgue);
  axis([li,re,0,6]); %ggf. weglassen oder anpassen
  if Fall==1
    xlabel(['(i) äqui, \lambda_n = ',num2str(lambda)]), else 
    xlabel(['(i) Tsch, \lambda_n = ',num2str(lambda)]);
  end; %if

  %--------------------------------------------------------------------------
  %Jetzt Raum (ii)
  disp('Jetzt Raum (ii)');
  Knotentau=zeros(1,n);
  if floor(n/2)==n/2 %n gerade
    J=2:2:n;
    Knotentau(J)=-J/2+li;
    Knotentau(J-1)=re-1+[2:(n+2)/2]; else %n ungerade
    J=2:2:n-1;
    Knotentau(J)=-J/2+li;
    Knotentau([J-1,n])=re-1+[2:(n+3)/2];
  end; %if
  for j=1:length(t)
    Knotenpolynom(j)=prod(t(j)-Knotentau);
  end; %for j
  Lagrange=zeros(size(L));
  for k=1:n
    [mindiff,j0]=min(abs(t-Knoten(k)));
    for j=1:length(t)
      Lagrange(j,k)=Knotenpolynom(j0)*L(j,k)/Knotenpolynom(j);
    end; %for j
    figure(1)
    subplot(4,n,k+n+2*n*(Fall-1));
    plot(t,Lagrange(:,k),[li,re],[0,0],Knoten,zeros(size(Knoten)),'o',...
    [Knoten(k),Knoten(k)],[0,1])
    axis([li,re,0,2]); %ggf. weglassen oder anpassen
  end; %for k
  Lebesgueii=sum(abs(Lagrange')); %Lebesgue-Funktion
  lambdaii=max(Lebesgueii);       %Lebesgue-Konstante
  [n,lambdaii] 
  figure(2)
  subplot(1,4,2+Fall)
  plot(t,Lebesgueii);
  axis([li,re,0,6]); %ggf. weglassen oder anpassen
  if Fall==1
    xlabel(['(ii) äqui, \lambda_n = ',num2str(lambdaii)]), else 
    xlabel(['(ii) Tsch, \lambda_n = ',num2str(lambdaii)]);
  end; %if
end; %for Fall

%===============================================================
%Ergebnisse Lebesgue-Konstante n=5 (im Polynomfall also Grad 4):
%Intervall: [0,1] (überall 1000 Punkte)
%                Äquidistante Knoten   T-Knoten
%Raum (i)            2.2078             1.9889
%Raum (ii)           2.5028             2.2184

%Intervall: [-1,1]
%                Äquidistante Knoten   T-Knoten
%Raum (i)            2.2078             1.9889
%Raum (ii)           2.9218             2.5869

%Intervall: [-5,5]
%                Äquidistante Knoten   T-Knoten
%Raum (i)            2.2078             1.9889
%Raum (ii)           8.5501             9.8837


%Ergebnisse Lebesgue-Konstante n=4 (im Polynomfall also Grad 3):
%Intervall: [-5,5]
%                Äquidistante Knoten   T-Knoten
%Raum (i)           1.6311              1.8478
%Raum (ii)          2.5326              5.0102

%In Raum (ii) sind die tau-Knoten immer rechts+j, links-j, j ganzzahlig
%bei [0,1]  also 2,-1,3,-2, etc,
%bei [-5,5] also 6,-6,7,-7, etc. 
