المبدأ

• نعتبر أن A متماثلة وموجبة.

• لنحديد حل النظمة الذي يعود إلى تحديد النقطة التي من أجلها تصل  إلى قيمتها الدنوية.

• ننتقل من المتجهةx^(k)  إلى المتجهةx^(k+1)  بالعلاقة : 

هو u^(k) الاتجاه المرافق لـ u⁽¹⁾و u⁽²⁾و...و u^(k-1)، ونختار   بحيث تكون  دنوية.

الخوارزم

التهيئة : x⁽⁰⁾ اعتباطية

التكرار :

مثال

لنحل النظمة Ax=b مع:

  و 

(نعتبر اعتباطيا المتجهة    )

حل يدوي (باستعمال الكسر)

البرمجة

• #include <math.h>
  #include <conio.h>
  #define NMAX 5
  
  double al_norme_vect(double x[NMAX],int n);
  
  void al_prod_mat_vect(double a[NMAX][NMAX],int n,double x[NMAX],double y[NMAX]);
  
  void sl_grad_conj(double a[NMAX][NMAX],double b[NMAX],double x1[NMAX],
                    int n,double t[NMAX][NMAX]);
  main()
  {
   double x[NMAX],b[NMAX],sol[NMAX];
   double a[NMAX][NMAX];
   double t[NMAX][NMAX];
   int i,j,n;
   clrscr();
   printf("Méthode du Gradient conjugué\n");
   n=4;
   b[1]=1;b[2]=0;b[3]=0;b[4]=-1;
   sol[1]=1;sol[2]=-1;sol[3]=1;sol[4]=-1;
   printf("                            A                                    b\n");
   for(i=1;i<=n;i++)
   {
    for(j=1;j<=n;j++)
    {
     a[i][j]=0;
     if (j==i+1) a[i][i+1]=1;
     if (j==i-1) a[i][i-1]=1;
     if (j==i)   a[i][j]=2;
     printf(" %12.6e ",a[i][j]);
    }
   x[i]=0;
   printf(" I  %12.6e\n",b[i]);
   }
   x[1]=1;
   printf("Solution exacte\n");
   for(i=1;i<=n;i++)
   printf("%22.16e\n",sol[i]);
   sl_grad_conj(a,b,x,n,t);
   for(j=NMAX-n;j<=NMAX-1;j++)
   {
    printf("Itération %d \n",j);
    for(i=1;i<=n;i++)
    {
     x[i]=t[i][j]-sol[i];
     printf("   %22.16e\n",t[i][j]);
    }
    printf("                     Norme erreur = %.3e\n",al_norme_vect(x,n));
   }
  }
  

  ---

  حلول البرنامج

   

  في كل كرة البرنامج يظهر:

  • الكرة ؛x^(k)

  • الارتياب (الفرق بين الكرةx^(k)  والقيمة المضبوطة)

  			    A                                    b
    2.00000e+00   1.00000e+00   0.00000e+00   0.00000e+00  I   1.00000e+00
    1.00000e+00   2.00000e+00   1.00000e+00   0.00000e+00  I   0.00000e+00
    0.00000e+00   1.00000e+00   2.00000e+00   1.00000e+00  I   0.00000e+00
    0.00000e+00   0.00000e+00   1.00000e+00   2.00000e+00  I  -1.00000e+00
  Solution exacte
   1.000000000000000e+00
  -1.000000000000000e+00
   1.000000000000000e+00
  -1.000000000000000e+00
  Itération 1 
      6.250000000000000e-01
     -3.750000000000000e-01
      0.000000000000000e+00
     -3.750000000000000e-01
  		     Norme erreur = 1.39e+00
  Itération 2 
      5.454545454545454e-01
     -4.545454545454545e-01
      6.363636363636364e-01
     -7.727272727272727e-01
  		     Norme erreur = 8.29e-01
  Itération 3 
      9.375000000000000e-01
     -8.124999999999999e-01
      8.750000000000001e-01
     -1.062500000000000e+00
  		     Norme erreur = 2.42e-01
  Itération 4 
      1.000000000000000e+00
     -1.000000000000000e+00
      1.000000000000000e+00
     -1.000000000000000e+00
  		     Norme erreur = 0.00e+00
  

double al_norme_vect(double x[NMAX],int n) 
{
 double t;
 int i;
 t=0;
 for(i=1;i<=n;i++) 
   t += x[i]*x[i]; 
 t=sqrt(t);
 return(t);
}

void al_prod_mat_vect(double a[NMAX][NMAX],int n,double x[NMAX],double y[NMAX]) 
{
 int i,j;
 for(i=1;i<=n;i++)
 {
  y[i]=0;
  for(j=1;j<=n;j++)
    y[i] += a[i][j]*x[j]; 
  } 
}

الدالة الأساسية sl_grad_conj

• معايير الدالة

a : المصفوفة
b : العنصر الثاني
n : درجة المصفوفة
x1: متجهة البداية
t: جدول


القيمة الرجعية
 تستدعي الدالة دالة أخرى وهي:

void al_prod_mat_vect(double a[NMAX][NMAX],int n,double x[NMAX],double y[NMAX])

ترجع الدالة القيم المكررة لهذه الطريقة في الجدول t.

تساوي الثابتة NMAX البعد القصوي للمصفوفة زائد 1.

void sl_grad_conj(double a[NMAX][NMAX],double b[NMAX],double x1[NMAX],
 int n,double t[NMAX][NMAX])
{
 int i,j,k;
 double alfa0,alfa1,s,mu,lambda;
 double r[NMAX],u[NMAX],v[NMAX];
 for(i=1;i<=n;i++)
 for(j=1;j<=NMAX-1;j++)
 t[i][j]=0.0;
 al_prod_mat_vect(a,n,x1,r);
 alfa0=0;
 for(i=1;i<=n;i++)
 {
 r[i] -= b[i];
 u[i]=-r[i];
 alfa0 += r[i]*r[i];
 }
 for(k=0;k<=n;k++)
 {
 al_prod_mat_vect(a,n,u,v);
 s=0;
 for(i=1;i<=n;i++)
 s += u[i]*v[i];
 mu=alfa0/s;
 for(i=1;i<=n;i++)
 {
 x1[i] += mu*u[i];
 for(j=1;j<=NMAX-2;j++)
 t[i][j]=t[i][j+1];
 t[i][NMAX-1]=x1[i];
 }
 al_prod_mat_vect(a,n,x1,r);
 alfa1=0;
 for(i=1;i<=n;i++)
 {
 r[i] -= b[i];
 alfa1 += r[i]*r[i];
 }
 if (alfa1==0) break;
 lambda=alfa1/alfa0;
 alfa0=alfa1;
 for(i=1;i<=n;i++)
 u[i]=-r[i]+lambda*u[i];
 }
}