tel/TPmatlab/interpolation/TP2f/diff_div.m

   1 function D=diff_div(X,Y)\r
   2 \r
   3 \r
   4 %       diff_div : calcul des différences divisées pour des points de supports quelconques\r
   5 %\r
   6 % *********************************************************\r
   7 %\r
   8 %       D=diff_div(X,Y)\r
   9 %     les différences divisées sont calculées par l'algorithme \r
  10 %     pyramidal donné en TD.\r
  11 %\r
  12 %       variables d'entrées : \r
  13 %     * X : contient les valeurs x_i, pour 0 <=i<=n (dans l'ordre croissant)\r
  14 %     * Y : contient les valeurs f^(r)(x_i), pour 0 <=i<=n, de telle sorte que :\r
  15 %             si x_k est multiple avec x_{k-1}~=x_k=x_{k+1}=...=x_{k+r}~=x_{k+r+1} alors\r
  16 %                y_k=f(x_k), y_{k+1}=f'(x_k), ..., y_{k+r}=f^(r)(x_k).\r
  17 %\r
  18 %       variables de sortie\r
  19 %     * D : contient les différences divisées (généralisée à support quelconque)         \r
  20 %           f[x_0], f[x_0,x_1], ...., f[x_0,x_1,...,x_n]\r
  21 %\r
  22 %\r
  23 %\r
  24 % ************ Fonctions auxiliaires utilisées ************\r
  25 %\r
  26 %   val_derivee_multiple     \r
  27 %\r
  28 % *********************************************************\r
  29 \r
  30 \r
  31 % Contrôles d'entrée\r
  32 %\r
  33 % nombre d'arguments\r
  34 if nargin~=2\r
  35    error('nombre d''arguments de la fonction incorrect');\r
  36 end\r
  37 % taille des deux variables \r
  38 n=length(X)-1;\r
  39 np=length(Y)-1;\r
  40 if n~=np\r
  41    error('les deux tableaux n''ont pas la même taille');\r
  42 end\r
  43 % vérification de l'ordre croissant des x_i \r
  44 % (si X non formel)\r
  45 if isnumeric(X)\r
  46    if max(abs(sort(X)-X))~=0\r
  47       error('support non ordonné');\r
  48    end\r
  49 end\r
  50 \r
  51 % Corps d'algorithme\r
  52 i=0;\r
  53 while (i<=n)\r
  54    [alpha,valfder]=val_derivee_multiple(i,X,Y);\r
  55    D(i+1:i+alpha)=valfder(1)*ones(1,alpha);\r
  56    i=i+alpha;\r
  57 end\r
  58 factoi=1;\r
  59 for i=1:n;\r
  60    factoi=factoi*i;\r
  61    for j=n:-1:i;\r
  62        dx=X(j+1)-X(j-i+1);\r
  63        if (dx==0)\r
  64           [alpha,valfder]=val_derivee_multiple(j-i,X,Y);\r
  65           D(j+1)=(valfder(i+1))/factoi;\r
  66        else\r
  67           D(j+1)=(D(j+1)-D(j))/(dx);\r
  68        end\r
  69     end\r
  70 end\r