طريقة تفكيك QR

طريقة تفكيك QR 

 QR
 

مصطلحات
العربية: طريقة التفكيك 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 */.

التعليقات   

 
# Guest 2014-11-10 09:24
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-10-04 06:16
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-09-09 22:07
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-26 23:37
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-18 04:58
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-15 06:56
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-08 16:46
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-08 11:50
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-06 23:42
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-06 23:07
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-06 17:47
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-08-03 09:13
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-07-24 22:01
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-07-13 07:10
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-07-01 05:22
قام المدير بحذف هذا التعليق
 
 
# Guest 2014-06-22 03:38
قام المدير بحذف هذا التعليق
 

أضف تعليق


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

المزيد في هذه الفئة :


Go to top