TP Scilab MP* N°2 & 3

// ne pas oublier funcprot(0) pour les fonctions récursives
funcprot(0);
 
 
// horner pour calculer p(t) récursion terminale
// les opérations sont effectuées dans l'&ccumulateur acc
function h = dhorner(p,t,acc)
  if p == [] then h = acc
  else h = dhorner(p(1:$-1),t,p($) + t * acc) //b_i = a_i + t*b_{i+1}
  end
endfunction
 
// on part d'un accumulateur nul
function h = horner(p,t)
  h = dhorner(p,t,0)
endfunction
 
// horner en gardant la trace des coefficients bi
// renvoie la liste des coefficients bi
// récursion terminale
// on rajoute une liste qui stocke les bi
function h = lhorner(p,t,acc_num,acc_liste)
  if p == [] then h = acc_liste
  else
    ai = p($)
    bi = ai + t * acc_num;
    h  = lhorner(p(1:$-1),t,bi,[bi,acc_liste])
  end
endfunction
 
// on part de 0 et d'une liste vide 
function h = liste_horner(p,t)
  h = lhorner(p,t,0,[])
endfunction
 
 
 
function y = quotient_horner(poly,t)
  // quotient de poly par X-t
  z = liste_horner(poly,t)
  y = [z(2:$)]
endfunction
 
function y = reste_horner(poly,t)
  // reste de poly par T-t
  z = liste_horner(poly,t)
  if z == [] then y=0
  else y = z(1)
  end
endfunction
 
function y = jeme_poly(poly,j,t)
  // calcul de poly^(j)(t) avec Horner
  if j == 0 then y = reste_horner(poly,t)
  else y = j*jeme_poly(quotient_horner(poly,t),j-1,t)
  end
endfunction
 
 
 
function y = jeme_poly_mat(poly,k,t)
  n = length(poly) + 1;
  T = zeros(k+1,length(poly)+1);
  for j = 2:n do
    T(1,j) = poly(j-1);
  end
  for i = 2:(k+1) do
    T(i,n-i+1) = T(i-1,n-i+2);
    for j = (n-i):-1:1 do
      T(i,j) = T(i-1,j+1) + t * T(i,j+1)
    end
  end
  y = T
endfunction
 
 
// N E W T O N 
 
 
function y = prod_mono(poly,t) 
  y = [0,poly] - t*[poly,0]
endfunction
 
function y = add_poly(p1,p2)
  l1 = length(p1)
  l2 = length(p2)
  l = max(l1,l2)
  if l1 < l then 
    for k = 1:(l-l1) do p1 = [p1,0];end
  else
    for k = 1:(l-l2) do p1 = [p2,0];end
  end
  y = p1+p2
endfunction
 
 
function P = newton(X,Y,t)
  n = length(Y)
  if n <> length(X) then 
    y = "Il manque des abscisses ou des ordonnées"
  end
  D = zeros(n,n)
  for i = 1:n do D(i,1) = Y(i); end 
  for j = 2:n do
    for i = 1:(n-j+1) do
      D(i,j) = (D(i+1,j-1) - D(i,j-1)) / (X(i+j-1) - X(i))
    end
  end
  P = D(1,n)
  for k = (n-1):-1:1 do
    P = D(1,k) + (t-X(k))*P
  end
endfunction
 
function P = newton_poly(X,Y)
  n = length(Y)
   if n <> length(X) then 
     y = "Il manque des abscisses ou des ordonnées"
  end
  D = zeros(n,n)
  for i = 1:n do D(i,1) = Y(i); end 
  for j = 2:n do
    for i = 1:(n-j+1) do
      D(i,j) = (D(i+1,j-1) - D(i,j-1)) / (X(i+j-1) - X(i))
    end
  end
  P = [D(1,n)]
  for k = (n-1):-1:1 do
    P = add_poly([D(1,k)] , prod_mono(P,X(k)))
  end
endfunction
 
 
function y  = poly2x(poly,x)
  n = length(poly)
  f = 0
  for k = 1:n do
    f = f + poly(k)*x^(k-1)
  end
  y = f
endfunction
 
 
for k = [0.1:0.1:1] do
  x = [0:k:2];
  p = newton_poly(x,1/(1+x^2));
  plot(x,poly2x(p,x))
end
 
 
//////////////// LAGRANGE
 
function y = lagrange(X,Y)
n=length(X);
P=[0];
for i=1:n, L=[1];
  for j=[1:i-1,i+1:n] L=prod_mono(L,X(j))/(X(i)-X(j));end
  P=add_poly(P,L*Y(i)); 
end
y=P
endfunction