كتبه نظمات المعادلات الخطية 7254
طريقة تفكيك QR طريقة تفكيك QR
قيم الموضوع
(0 أصوات)

 

 
 
مصطلحات
العربية: طريقة التفكيك QR
الإنجليزية: QR decomposition (QR factorization)i
الفرنسية: Décomposition QR

المبدأ

  • في كل مرحلة  من طريقة هوزهولد، ننشئ مصفوفة H(k)  متماثلة ومتعامدة، بحيث A(k+1) = H(k)A(k). 

المصفوفة R = A(n) مثلثية علوية بحيث A = H(1)H(2)... H(n-1)R =QR

المصفوفة Q متعامدة.

  • طريقة تفكيك QR  تنبني على إنشاء المصفوفة Q بالإضافة إلى المصفوفة R، ثم حل النظمة المثلثية Rx = QTb.

الخوارزم

المرحلة 1: حساب Q و R  بحيث  A = QR ، حيث Q متعامدة و R مثلثية علوية

A(1) = A, b(1) = b.

Q = H(1) H(2)... H(n-1)   و  R = A(n)

 

المرحلة 2: حل النظمة المثلثية  Rx = QTb

y = QTb

مثال

لنحل النظمة  Ax=b بحيث:

و

الحلول الصحيحة هي : x1 = 32/7,  x2 = 47/14,  x3 = 16/7,  x4 = -18/7.


بداية التفكيك (باستعمال الكسور و الأصول)

 المرحلة 1: تحويل العمود الأول للمصفوفة A، ثم حساب H(1)

  • = -2
    = -2 =(-2- 8) = 372 + 16
    = 8 + 2


             = (8 + 2)(-4) + (4)(2) + (-16)(6) + (6)(10) = -60 - 8

             = (8 + 2)(3) + (4)(-6) + (-16)(-2) + (6)(-15) = -58 + 6

             = (8 + 2)(7) + (4)(4) + (-16)(-15) + (6)(10) = 372 + 14










    H = G + I


البرمجة

#include <math.h>
#include <conio.h>
#define NMAX 5

void al_prod_mat(double a[NMAX][NMAX],double b[NMAX][NMAX],double c[NMAX][NMAX],int n);
int sl_decomp_qr(double a[NMAX][NMAX],double q[NMAX][NMAX],int n);
void sl_resol_qr(double q[NMAX][NMAX],double r[NMAX][NMAX],double b[NMAX],double x[NMAX],int n);

main()
{
 double a[NMAX][NMAX],q[NMAX][NMAX];
 double b[NMAX],x[NMAX];
 int i,j,n;
 int err;
 n=4;
 clrscr();
 printf("Méthode de la Décomposition QR\n");
 a[1][1]=8;a[1][2]=-4;a[1][3]=3;a[1][4]=7;b[1]=12;
 a[2][1]=4;a[2][2]=2;a[2][3]=-6;a[2][4]=4;b[2]=1;
 a[3][1]=-16;a[3][2]=6;a[3][3]=-2;a[3][4]=-15;b[3]=-19;
 a[4][1]=6;a[4][2]=10;a[4][3]=-15;a[4][4]=10;b[4]=1;
 printf("                             A(1)\n");
 for(i=1;i<=n;i++)
 {
  for(j=1;j<=n;j++)
  printf("%14.8e  ",a[i][j]);
  printf("\n");
 }
 err=sl_decomp_qr(a,q,n);

  if (err==1)
 {
  printf("                              Q\n");
  for(i=1;i<=n;i++)
  {
   for(j=1;j<=n;j++)
   printf("%14.8e  ",q[i][j]);
   printf("\n");
  }
  sl_resol_qr(q,a,b,x,n);
  printf("Solution:\n");
  for(j=1;j<=n;j++)
  printf("x%d = %22.16e\n",j,x[j]);
  printf("\n");
 }
 else printf("Erreur, matrice singulière\n");
} 

نتائج البرنامج

نزيل التعليقات من الدالة sl_decomp_qr وفي كل مرحلة يظهر البرنامج :

  •  المصفوفة  A(k) 

  • العنصر الثاني b(k)

و يظهر بعد ذلك R, Q وحلول النظمة.

     A(1)
 8.0000000e+00  -4.0000000e+00   3.0000000e+00   7.0000000e+00  
 4.0000000e+00   2.0000000e+00  -6.0000000e+00   4.0000000e+00  
-1.6000000e+01   6.0000000e+00  -2.0000000e+00  -1.5000000e+01  
 6.0000000e+00   1.0000000e+01  -1.5000000e+01   1.0000000e+01  
     H(1)
-4.1478068e-01  -2.0739034e-01   8.2956136e-01  -3.1108551e-01  
-2.0739034e-01   9.6959900e-01   1.2160401e-01  -4.5601506e-02  
 8.2956136e-01   1.2160401e-01   5.1358394e-01   1.8240602e-01  
-3.1108551e-01  -4.5601506e-02   1.8240602e-01   9.3159774e-01  
     A(2)
-1.9287302e+01   3.1108551e+00   3.0071599e+00  -1.9287302e+01  
 0.0000000e+00   3.0423684e+00  -5.9989504e+00   1.4658833e-01  
 0.0000000e+00   1.8305265e+00  -2.0041982e+00   4.1364667e-01  
 0.0000000e+00   1.1563553e+01  -1.4998426e+01   4.2198825e+00  
     H(2)
 1.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00  
 0.0000000e+00  -2.5151050e-01  -1.5132836e-01  -9.5595094e-01  
 0.0000000e+00  -1.5132836e-01   9.8170189e-01  -1.1559031e-01  
 0.0000000e+00  -9.5595094e-01  -1.1559031e-01   2.6980861e-01  
     A(3)
-1.9287302e+01   3.1108551e+00   3.0071599e+00  -1.9287302e+01  
 0.0000000e+00  -1.2096387e+01   1.6149850e+01  -4.1334656e+00  
 0.0000000e+00   0.0000000e+00   6.7395888e-01  -1.0388279e-01  
 0.0000000e+00   0.0000000e+00   1.9196638e+00   9.5061582e-01  
     H(3)
 1.0000000e+00   0.0000000e+00   0.0000000e+00   0.0000000e+00  
 0.0000000e+00   1.0000000e+00   0.0000000e+00   0.0000000e+00  
 0.0000000e+00   0.0000000e+00  -3.3125953e-01  -9.4353968e-01  
 0.0000000e+00   0.0000000e+00  -9.4353968e-01   3.3125953e-01  
   R = A(4)
-1.9287302e+01   3.1108551e+00   3.0071599e+00  -1.9287302e+01  
 0.0000000e+00  -1.2096387e+01   1.6149850e+01  -4.1334656e+00  
 0.0000000e+00   0.0000000e+00  -2.0345343e+00  -8.6253158e-01  
 0.0000000e+00   0.0000000e+00   0.0000000e+00   4.1291809e-01  
      Q
-4.1478068e-01   2.2400717e-01  -3.0947113e-01  -8.2583618e-01  
-2.0739034e-01  -2.1867366e-01   9.0674174e-01  -2.9494149e-01  
 8.2956136e-01  -2.8267571e-01  -3.4674637e-02  -4.8033329e-01  
-3.1108551e-01  -9.0669568e-01  -2.8433202e-01   1.6853800e-02  
Solution:
x1 =  4.571428571428571e+00
x2 =  3.357142857142856e+00
x3 =  2.285714285714285e+00
x4 = -2.571428571428571e+00

 

void al_prod_mat(double a[NMAX][NMAX],double b[NMAX][NMAX],double c[NMAX][NMAX],int n) 
{
 double s;
 int i,j,k;

 for(i=1;i<=n;i++)
   for(j=1;j<=n;j++)
    {
      s=0;
      for(k=1;k<=n;k++)
        s += a[i][k]*b[k][j]; c[i][j]=s; 
    }
}

 الدالة الرئيسية sl_decomp_qr


تستدعي هذه الدالة الدالتين التاليتين:

الدالة sl_decomp_qr 

a : المصفوفة الرئيسية
q : جدول
n : درجة المصفوفة
وهي تستدعي الدالة 

 void al_prod_mat(double a[NMAX][NMAX],double b[NMAX][NMAX],double c[NMAX][NMAX],int n)

 التي يجب أن تعرف في البرنامج الرئيسي أيضا.

ترجع القيمة 1 عندما تكون المصفوفة مقلوبة رقميا والقيمة 0 إذا كان عكس ذلك.

سيكون الحل في الجدول q الذي سيعطي التفكيك QR للمصفوفة الرئيسية.

الدالة sl_resol_qr

  • q و r : يمثلان المصفوفتين Q و R
  •  b : جدول العنصر الثاني
  • x : جدول أحادي البعد
  • n : درجة المصفوفة
ترجع بالجدول x الذي يمثل الحل.

 

 

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

int sl_decomp_qr(double a[NMAX][NMAX],double q[NMAX][NMAX],int n)
{
 int i,j,k,sgn,err;
 double s,nu;
 double aa[NMAX][NMAX],qq[NMAX][NMAX],h[NMAX][NMAX];
 double bb[NMAX],v[NMAX],gam[NMAX];
 for(i=1;i<=n;i++)
 for(j=1;j<=n;j++)
 if (i==j) q[i][j]=1; else q[i][j]=0;
 err=1;
 k=1;
 while (err==1 && k<n)
 {
 for(i=1;i<=n;i++)
 {
 for(j=1;j<=n;j++)
 aa[i][j]=0;
 v[i]=0;
 }
 s=0;
 for(i=k;i<=n;i++)
 s+=a[i][k]*a[i][k];
 if (a[k][k]<0) sgn=-1; else sgn=1;
 aa[k][k]=-sgn*sqrt(s);
 v[k]=a[k][k]-aa[k][k];
 nu=-aa[k][k]*v[k];
 if (fabs(nu)<1e-15) err=0;
 else
 {
 for(i=k+1;i<=n;i++)
 v[i]=a[i][k];
 for(j=k+1;j<=n;j++)
 {
 s=0;
 for(i=k;i<=n;i++)
 s+=v[i]*a[i][j];
 gam[j]=s/nu;
 }
 for(i=k;i<=n;i++)
 for(j=k+1;j<=n;j++)
 aa[i][j]=a[i][j]-gam[j]*v[i];
 for(i=1;i<=k-1;i++)
 for(j=1;j<=n;j++)
 if (j==i) h[i][j]=1; else h[i][j]=0;
 for(i=k;i<=n;i++)
 for(j=1;j<=k-1;j++)
 h[i][j]=0;
 for(i=k;i<=n;i++)
 for(j=k;j<=n;j++)
 if (j==i) h[i][j]=1-v[i]*v[j]/nu; else h[i][j]=-v[i]*v[j]/nu;
 al_prod_mat(q,h,qq,n);
 for(i=k;i<=n;i++)
 for(j=k;j<=n;j++)
 a[i][j]=aa[i][j];
 for(i=1;i<=n;i++)
 for(j=1;j<=n;j++)
 q[i][j]=qq[i][j];
 /*
 printf(" H(%d)\n",k);
 for(i=1;i<=n;i++)
 {
 for(j=1;j<=n;j++)
 printf("%14.8e ",h[i][j]);
 printf("\n");
 }
 */
 k++;
 /*
 if(k==n)
 printf(" R = A(%d)\n",n);
 else printf(" A(%d)\n",k);
 for(i=1;i<=n;i++)
 {
 for(j=1;j<=n;j++)
 printf("%14.8e ",a[i][j]);
 printf("\n");
 }
 */
 }
 }
 if(fabs(a[n][n])<1e-15) err=0;
 return(err);
}

void sl_resol_qr(double q[NMAX][NMAX],double r[NMAX][NMAX],double b[NMAX],double x[NMAX],int n)
{
 int i,j;
 double s;
 double y[NMAX];
 for(i=1;i<=n;i++)
 {
 y[i]=0;
 for(j=1;j<=n;j++)
 y[i]+=q[j][i]*b[j];
 }
 x[n]=y[n]/r[n][n];
 for(i=n-1;i>0;i--)
 {
 s=y[i];
 for(j=i+1;j<=n;j++)
 s-=r[i][j]*x[j];
 x[i]=s/r[i][i];
 }
}
Remarque : si on veut que la fonction affiche les matrices intermédiaires, il suffit de supprimer les marques de commentaires /* et */.
مقالات أخرى من نفس الفئة « طريقة هوزهولدر طريقة يعقوبي »

أضف تعليقا

إصنعها يريد أن يتأكد أنك لست روبوتا، لذلك أحسب ما يلي:

كود امني
تحديث

عد إلى الأعلى

Main Menu-burger