المبدأ

نحسب حدودية الاستكمال P_n(x) تحت صيغة تآلفية الخطية لحدودية:

P_n(x) = D_0,0 + (x-x₀)D_0,1 + (x-x₀)(x-x₁)D_0,2 +... + (x-x₀)...(x-x_n-1)D_0,n

عوامل D_i,j هي الفوارق المقسمة، محسوبة بالترجع.

الخوارزم

حساب الفوارق المقسمة

نهيئ أولا ما يلي:

 D_i,i = y_i,     i = 0,... n

تكرار

حساب حدوديات الاستكمال

b₀(x) = D_0,n

b₁(x) = D_0,n-1 + b₀(x) (x-x_n-1)

b₂(x) = D_0,n-2 + b₁(x) (x-x_n-2)

....
....

b_n-1(x) = D_0,1 + b_n-2(x) (x-x₁)

b_n(x) = D_0,0 + b_n-1(x) (x-x₀)             P_n(x) = b_n(x)

مثال

نعتبر حدودية الاستكمال P(x)  بحيث:

P(0) = 0, P(1)=3, P(3)=5, P(9/2)=4, P(5)=3

احسب قيمةP(x)  في 20 نقطة متباعدة بنفس المسافة من المجال [0, 5]

حساب P(2)

x0=0, x1=1, x2=3, x3=9/2, x4=5
y0=0, y1=3, y2=5, y3=4, y4=3 
D0,0 = y0 = 0
D1,1 = y1 = 3
D2,2 = y2 = 5
D3,3 = y3 = 4
D4,4 = y4 = 3 
D0,1 = (3-0)/(1-0) = 3
D1,2 = (5-3)/(3-1) = 1
D2,3 = (4-5)/(9/2-3) = -2/3
D3,4 = (3-4)/(5-9/2) = -2 
D0,2 = (1-3)/(3-0) = -2/3
D1,3 = (-2/3-1)/(9/2-1) = -10/21
D2,4 = (-2-(-2/3))/(5-3) = -2/3 
D0,3 = (-10/21-(-2/3))/(9/2-0) = 8/189
D1,4 = (-2/3-(-10/21))/(5-1) = -1/21 
D0,4 = (-1/21-8/189)/(5-0) = -17/945 
b0(2) = -17/945
b1(2) = (-17/945)*(2-9/2)+8/189 = 11/126
b2(2) = (11/126)*(2-3)-2/3 = -95/126
b3(2) = (-95/126)*(2-1)+3 = 283/126
b4(2) = (283/126)*(2-0)+0 = 283/63 
P(2) = b4(2) = 283/63 

البرمجة

#include <math.h>
#include <conio.h>
#define NMAX 5
 
void it_tabdifdiv(int n,double x[NMAX],double f[NMAX],double d[NMAX][NMAX]);
double it_pol_difdiv(int n,double d[NMAX][NMAX],double x[NMAX],double alpha);

main()
{
 double alpha,ord;
 int i,j;
 double x[NMAX],f[NMAX];
 double d[NMAX][NMAX];
 clrscr();
 printf("Méthode de Newton\n");
 x[0]=0;x[1]=1;x[2]=3;x[3]=4.5;x[4]=5;
 f[0]=0;f[1]=3;f[2]=5;f[3]=4;f[4]=3;
 it_tabdifdiv(4,x,f,d);
 printf("Tableau des différences divisées\n");
 for(i=0;i<=4;i++)
 {
  for(j=0;j<=4-i;j++)
   printf("%14.5e",d[i][j]);
  printf("\n");
 }
 printf("Calcul du polynôme d'interpolation\n");
 for(j=0;j<=20;j++)
 {
  alpha=0.25*j;
  ord=it_pol_difdiv(4,d,x,alpha);
  printf("x =%6.2f     y =%24.17e\n",alpha,ord);
 }
}

نتيجة البرنامج

البرنامج يظهر جدول الفوارق المقسمة ثم حساب: قيمة المتغير وقيمة الملائمة لحدودية الاستكمال.

جدول الفوارق المقسمة
    0.0000e+00    3.0000e+00   -6.6667e-01    4.2328e-02   -1.7989e-02
    3.0000e+00    1.0000e+00   -4.7619e-01   -4.7619e-02
    5.0000e+00   -6.6667e-01   -6.6667e-01
    4.0000e+00   -2.0000e+00
    3.0000e+00

حساب  دالة الاستكمال
x =  0.00     y =  0.0000000000000000e+00
x =  0.25     y =  9.3624751984126986e-01
x =  0.50     y =  1.7380952380952381e+00
x =  0.75     y =  2.4213169642857144e+00
x =  1.00     y =  3.0000000000000000e+00
x =  1.25     y =  3.4865451388888893e+00
x =  1.50     y =  3.8916666666666666e+00
x =  1.75     y =  4.2243923611111107e+00
x =  2.00     y =  4.4920634920634921e+00
x =  2.25     y =  4.7003348214285721e+00
x =  2.50     y =  4.8531746031746028e+00
x =  2.75     y =  4.9528645833333336e+00
x =  3.00     y =  5.0000000000000000e+00
x =  3.25     y =  4.9934895833333339e+00
x =  3.50     y =  4.9305555555555562e+00
x =  3.75     y =  4.8067336309523814e+00
x =  4.00     y =  4.6158730158730163e+00
x =  4.25     y =  4.3501364087301599e+00
x =  4.50     y =  4.0000000000000009e+00
x =  4.75     y =  3.5542534722222219e+00
x =  5.00     y =  3.0000000000000004e+00

الدالتين it_pol_difdiv و it_tabdifdiv

يجب إستدعاء هاتين الدالتين بشكل متتابع:

• n : عدد نقاط الإستكمال
  x : جدول يحتوي على أفاصيل الإستكمال
  f : جدول يحتوي على قيم f(x)i لنقاط الإستكمال
  d : جدول يضم معاملات الحدودية في البداية
  
  تكون النتائج في الجدول d بحيث يتم حساب الفوارق المقسمة المرتبطة بالدالة f(x)i، التي بدورها يتم تزويدها بقيم المعيار f، من أجل n+1 قيمة مختلفة للمتغير المعطى من خلال المعيار x.

معايير الدالة it_pol_difdiv:

• n :عدد نقاط الإستكمال
  d : جدول يحتوي على معاملات حدوديات القاعدة
  x : جدول يحتوي على أفاصيل الإستكمال
  alpha : قيمة المتغير

ترجع هذه الدالة، بخصوص القيمة المعطاة للمتغير، قيمة حدودية الإستكمال المحسوبة انطلاقا من عناصر جدول الفوارق المقسمة.

الثابتة الصحيحة NMAX تساوي العدد القصوي لنقاط الإستكمال.

void it_tabdifdiv(int n,double x[NMAX],double f[NMAX],double d[NMAX][NMAX]) 
{
 int i,j;
 for(j=0;j<=n;j++) d[j][0]=f[j];

 for(j=1;j<=n;j++)
  for(i=0;i<=n-j;i++)
   d[i][j]=(d[i+1][j-1]-d[i][j-1])/(x[i+j]-x[i]);
 }

 double it_pol_difdiv(int n,double d[NMAX][NMAX],double x[NMAX],double alpha)
 {
  int i;
  double p;

  p=d[0][n];
  for(i=1;i<=n;i++)
    p=p*(alpha-x[n-i])+d[0][n-i];

  return(p); 
}