% orthclas.m
%
% polynomes orthogonaux classiques
%
class=1;while class > 0,
 class=input(' tapez 1 si Jacobi, 2 si Laguerre, 3 si Hermite, 0:stop ');
 if class<=0 ; break ; end;
 if class>3 ; break ; end;
 clear al;clear be;
 if class < 3 ; alp=input(' alpha ? ');end;
 if class == 1 ; bet = input(' beta ?');end;
 n=input(' degre ? ');
 tit1=sprintf(' degre %g ',n);tit2=' ';tit3= ' ';
 if class==1,
  clear zn;zn=(0:n-1);clear den;den=2*zn+alp+bet;
  al=(bet*bet-alp*alp)*ones(1,n)./(  den.*(den+2) );
  clear zn;zn=(2:n-1);clear den;den=2*zn+alp+bet;
  be(3:n)=4*zn.*(zn+alp).*(zn+bet).*(zn+alp+bet)./((den-1).*den.*den.*(den+1));
  be(1)=0;be(2)=4*(alp+1)*(bet+1)/((alp+bet+2)^2 *(alp+bet+3) );
  a=-1;b=1;w='(1-x).^alp .* (1+x).^bet';
  tit0=' Jacobi ';tit2=sprintf(' alpha= %g ',alp);tit3=sprintf(' beta= %g',bet);
 end;
 if class==2,
  al=2*(0:n-1)+alp+1;be=(0:n-1).*((0:n-1)+alp);
  a=0;b=n+1;w='x.^alp.*exp(-x)';
  tit0=' Laguerre ';tit2=sprintf(' alpha= %g ',alp);
 end;
 if class==3,
  al=0*ones(1,n);be=(0:n-1)/2;
  a=-sqrt(n/2)-1;b=-a;w='exp(-x.^2)';
  tit0=' Hermite ';
 end;
 disp(al);
 disp(be);
 nx=20;x=a+0.001+(0:nx)*0.997*(b-a)/nx;
 p1=ones(1:nx+1);p0=0*p1;
 for m=1:n;
  % degre m
  p2=(x-al(m)).*p1-be(m)*p0;
  p0=p1;p1=p2;
 end;
 tit=[tit0 tit1 tit2 tit3];
 wx=eval(w);
 y1=p1/norm(p1,inf);y2=wx/norm(wx,inf);y3=sqrt(wx).*p1;y3=y3/norm(y3,inf);
 plot(x,y1,x,y2,x,y3);
 title(tit);
 meta orthclassg;
 % gpp orthclassg.met -dps   -> fichier ps
end
