#include <math.h>
#include <conio.h>
#define ITERMAX 7

double qu_f(double x)
{
 return(1/(x+1));
}

void qu_trapezes_adapt(double a,double b,double eps,double t[ITERMAX]);
void qu_romberg(double a,double b,double c[ITERMAX][ITERMAX]);
 
main()
{
 int i,k;
 long n;
 double a,b,r,vr_val;
 double c[ITERMAX][ITERMAX];
 clrscr();
 a=0;b=1;
 vr_val=log(2);
 printf("Méthode de Romberg\n");
 printf("Vraie valeur = %.17e\n",vr_val);
 n=1;
 qu_romberg(a,b,c);
 printf("\n  n              Trapèzes              Erreur\n");
 for(i=0;i<ITERMAX;i++)
 {
  printf("%4ld      %.17e      %.3e  \n",n,c[0][i],fabs(c[0][i]-vr_val));
  n*=2;
 }
 for(k=1;k<ITERMAX;k++)
 {
  printf("\n  n          Extrapolation%3d          Erreur\n",k);
  n=pow(2,k-1);
  for(i=0;i<ITERMAX-k;i++)
  {
   n *= 2;
   r=c[k][i];
   printf("%4ld      %.17e      %.3e\n",n,r,fabs(r-vr_val));
  }
 }
}

القيمة الحقيقية= 6.9314718055994529e-01

   n           شبه المنحرف             الارتياب
   1      7.5000000000000000e-01      5.69e-02  
   2      7.0833333333333326e-01      1.52e-02  
   4      6.9702380952380949e-01      3.88e-03  
   8      6.9412185037185026e-01      9.75e-04  
  16      6.9339120220752681e-01      2.44e-04  
  32      6.9320820826924900e-01      6.10e-05  
  64      6.9316243888340345e-01      1.53e-05  

   n          الاستكمال                 الارتياب                   
   1  
   2      6.9444444444444431e-01      1.30e-03
   4      6.9325396825396823e-01      1.07e-04
   8      6.9315453065453048e-01      7.35e-06
  16      6.9314765281941904e-01      4.72e-07
  32      6.9314721028982307e-01      2.97e-08
  64      6.9314718242145490e-01      1.86e-09

   n       الاستكمال                    الارتياب                   
   2  
   4      6.9317460317460311e-01      2.74e-05
   8      6.9314790148123462e-01      7.21e-07
  16      6.9314719429707827e-01      1.37e-08
  32      6.9314718078784998e-01      2.28e-10
  64      6.9314718056356373e-01      3.62e-12

   n              الاستكمال             الارتياب                   
   3 
   8      6.9314747764483198e-01      2.97e-07
  16      6.9314718307193290e-01      2.51e-09
  32      6.9314718057341773e-01      1.35e-11
  64      6.9314718056000357e-01      5.83e-14

  n            الاستكمال               الارتياب                   
  4  
  16      6.9314718191674507e-01      1.36e-09
  32      6.9314718056361968e-01      3.67e-12
  64      6.9314718055995095e-01      5.66e-15

   n              الاستكمال            الارتياب                   
   5  
  32      6.9314718056229696e-01      2.35e-12
  64      6.9314718055994740e-01      2.11e-15

  n               الاستكمال            الارتياب                   
  6  
  64      6.9314718055994684e-01      1.55e-15

الدالة التكميلية qu_trapezes_adapt

معايير هذه الدالة هي كالتالي:

a و b: حدي المجال [a, b]
eps : الدقة المطلوبة
t : جدول ذي الدرجة ITERMAX.

تستدعي الدالة الرئيسة دالة أخرى هي كالتالي:

 double qu_f(double x)

 التي يجب أن تعرف في البرنامج الأساسي وهي تمثل الدالة الرياضية f.

ترجع الدالة في الجدول t التكرارات المتتالية لطريقة شبه المنحرف.

الثابتة الصحيحة ITERMAX تساوي العدد القصوي للتكرارات المطلوبة (+1).

void qu_trapezes_adapt(double a,double b,double eps,double t[ITERMAX])
{
 int j;
 long i,n;
 double h,s,err,x;

 n=1;
 h=b-a;
 j=0; err=2*eps;
 t[0]=h*(qu_f(a)+qu_f(b))/2;

 while((fabs(err)>eps) && (j < ITERMAX-1))
 {
  j++;
  n *= 2;
  h /= 2;
  x=a-h;
  s=0;

  for(i=1;i<=n/2;i++) 
  {
   x += 2*h;
   s += qu_f(x);
  }

  t[j]=t[j-1]/2+h*s;
  err=t[j]-t[j-1];
 } 
}  

 الدالة التكميلية qu_romberg

معايير هذه الدالة هي كالتالي:
a و b: حدي المجال [a, b]
c: جدول ذي الدرجة ITERMAX*ITERMAX.

تستدعي الدالة الرئيسة دالة أخرى هي كالتالي:

void qu_trapezes_adapt(double a,double b,double eps,double t[ITERMAX])

وهي بدورها تستدعي الدالة التالية:

double qu_f(double x)

التي يجب أن تعرف في البرنامج الأساسي وهي تمثل الدالة الرياضية f وترجع في سطور الجدول c قيم الإستكمالات المتتابعة المولد من طرف طريقة رومبرغ.

الثابتة الصحيحة ITERMAX تساوي عدد الإستكمالات المطلوبة (+1).

void qu_romberg(double a,double b,double c[ITERMAX][ITERMAX])
{
 int i,k;
 double eps;
 eps=1e-14;
 qu_trapezes_adapt(a,b,eps,c);
 for(k=1;k<ITERMAX;k++)
   for(i=0;i<ITERMAX-k;i++)
      c[k][i]=(pow(4,k)*c[k-1][i+1]-c[k-1][i])/(pow(4,k)-1); 
}