//SECTION 1

// Question 1
// P0=1 puis par recurrence :
// Orth de Schmidt appliqué à (P0,...,Pn,X*Pn)   
// En remarquant que Pn est ortho à X P{n-2}
//Question 2
// On commence par définir une fonction omega que l'on pourra
// modifier au besoin. On commence par la plus simple possible.
// La fonction constante définie sur l'intervalla [0,1]
function S=omega(x)
    S=1
endfunction
a=-1
b=1

 // Il nous faut ensuite définir le produit scalaire.
 
function S=scal(P,Q)
    h=0.0001 // ce h correspond à un pas de précision que l'on peut ajuster pour 
            // obtenir la précision voulue du résultat final.
    x=a:h:b // On définit le vecteur x=[a,a+h,a+2h,...,b-h]= (a+ZZ*h) /\ [a,b[
    // Pour évaluer un polynome en un point, on utilise la commande horner
    // du nom de la méthode de horner qui permet de minimiser le nombre de 
    //multiplications nécessaires. 
    // Si P(x)=x^3+6x^2+2*x+5, on écrit  P(x)=x(x(x +6) +2)+5 
    // puis horner(P(x),3)=3*(3*(3+6)+2)+5
    S=h*omega(x).*horner(P*Q,x) //evalue omega et P*Q pour toutes les valeurs de
                                // x simultanément. On obtient l'aire de chaque       
                                //rectangle.
    S=sum(S)                    // Calcule la somme des valeurs obtenue c'est à
                                // dire la somme de l'aire des rectange 
endfunction

X=poly(0,"X")   // Définit une constante X générateur de l'algèbre des polynômes

scal(X,X)  // On fait un test pour s'assurer que l'on obtient bien la valeur 
           // escomptée.

function Q=ortho(n)
    P=1
    Q=1
    if n>0 then
        Q=X-scal(P,X*P)/scal(P,P)*P
    end
    for i=2:n
        R=X*Q-scal(P,X*Q)/scal(P,P)*P--scal(Q,X*Q)/scal(Q,Q)*Q
        P=Q
        Q=R
    end
    
endfunction

//Question 3
x=(-1:0.001:1)';
scf(0) // Sélectionne la fenetre de plot 0 comme fenetre courante
clf()  // efface la fenetre courante

// représentation graphique des trois polynômes orthogonaux
plot(x,[horner(ortho(10),x),horner(ortho(9),x),horner(ortho(8),x)]) 


scf(1) // La représentation graphique suivante sera dans la fentre 1.
//Question 4 
// Lafonction roots de scilab donne les racines d'un polynôme
roots(ortho(10))

//Question 5
roots(X^6-14*X^5+70*X^4-140*X^3+49*X^2+154*X-120)
//Question 6
function S=orthoroots(n)
    S=roots(ortho(n))'
endfunction


//SECTION 2

orthoroots(10)
function S=foo(x)
    S=x
endfunction

//Question 7

// Fonction d'interpolation par la méthode de Newton
// Prenant en entrée les points d'interpolation
function S=Newton(x)
    X=poly(0,"X")  // Le générateur de l'algèbre des polynômes
    n=size(x,'*')  // le nombre de points d'interpolation
    F=zeros(n,n)   // On initialise la matrice des coefficients de Newton
    F(:,1)=foo(x)' // La première colonne est donnée par les valeurs de la fonction
                   // au points d'interpolation
                   // Puis par récurrence on applique la formule du cours
                   // /!\   le dénominateur est un traitre car x0 et xk deviennent 
                   //       xi et xi+k
    for j=2:n
        F(j:n,j)=   (F(j:n,j-1)-F((j-1):(n-1),j-1))./((x(j:n)-x(1:(n - j+1)))')
    end
    P=diag(F)       //Seule la diagonale des coef de Newtons nous intéresse
    S=zeros(n,1)    // S contiendra les polynômes d'interpolation 
    Q=zeros(n,1)    // Q contient les polynômes de Newtons
    Q(1)=1          // Le premier polynôme de Newton est 1
    S(1)=sum(P.*Q)  // Pour otbenir le polynome d'interpolation, 
                    // Il suffit de multiplier terme à terme la liste 
                    //des polynomes de Newton et la liste des coefficients de Newton
                    
    for i=2:n       // On créé la liste des polynômes de newton et d'interpolation 
                    // par récurrence.
        Q(i)= prod(X*ones(i-1)-x(1:(i-1))) 
        S(i)=sum(P.*Q)
    end
    S //On renvoie la liste des polynômes d'interpolations
endfunction

// En utilisant les polynômes de lagrange
function P=lagrange(x)
    n=size(x,'*')
    X=poly(0,"X")
    P=ones(n,1)
    for i=1:n
        for j=1:n
            if j<>i then
              P(i)=P(i)*(X-x(j))/(x(i)-x(j))  
            end
        end
    end
    
endfunction

//On peut tester sur un exemple
//Question 8 


function S=foo(x)
    S=ones(x)./(1+25*x.^2)
endfunction
a=-1
b=1
pas=0.2
y=a:pas:b
n=size(y,'*')
P=Newton(y);
Q=lagrange(y);
scf(1)
clf()
PN=P(n);
PL=PN;
x=-1:0.01:1;
plot2d(x',[horner(PN(1),x'),foo(x')],rect=[-1,-1,1,1])





