مصطلحات

العربية: طريقة التفكيك LU
الإنجليزية: Householder Method
الفرنسية: Méthode de Householder

المبدأ

• يساعد الخوارزم على الحصول على نظمة مثلثية مكافئة والى حل للنظمة.

• لتكن  A⁽¹⁾ = A  و b⁽¹⁾ = b.

• نضرب (في اليسار) الجدول(A⁽¹⁾ | b⁽¹⁾)  في مصفوفة متعامدة ومتماثلة على الشكل التالي :

H⁽¹⁾ = I - 2u.u^T نختار u بحيث تكون المصفوفةH⁽¹⁾A⁽¹⁾   كما يلي :

   و   

النظمة الأصلية مكافئة للنظمة الجديدة A⁽²⁾x = b⁽²⁾

• نضرب بعد ذلك(A⁽²⁾ | b⁽²⁾) في المصفوفةH⁽²⁾  (متعامدة و متماثلة، تم إنشاءها كالمصفوفة H⁽¹⁾) بطريقة تنعدم فيها العناصر ما فوق القطر للعمود الثاني، وهكذا حتى تصبح المصفوفة مثلثية علوية.

• ونحل بعد ذلك النظمة المثلثية.

الخوارزم

 المرحلة 1:  التحويل المثلثي للمصفوفة A، أي تحويل النظمة إلى نظمة مكافئة ذات مصفوفة مثلثية علوية.

A⁽¹⁾ = A, b⁽¹⁾ = b.

المرحلة 2: حل النظمة المثلثية A⁽ⁿ⁾ x = b⁽ⁿ⁾

مثال

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

  و 

الحلول الصحيحة هي :  x₁ = 32/7,  x₂ = 47/14,  x₃ = 16/7,  x₄ = -18/7.

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

 المرحلة 1: التحويل المثلثي لـ A، تحويل العمود الأول

= -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

         = (8 + 2)(12) + (4)(1) + (-16)(-19) + (6)(1) = 410 + 24

البرمجة

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

int sl_householder_aff(double a[NMAX][NMAX],double b[NMAX],int n);

main()
{
 double a[NMAX][NMAX];
 double b[NMAX];
 int i,j,n;
 int err;
 n=4;
 clrscr();
 printf("Méthode de Householder\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)                               b(1)\n");
 for(i=1;i<=n;i++)
 {
  for(j=1;j<=n;j++)
  printf("%12.6e ",a[i][j]);
  printf(" I  %12.6e",b[i]);
  printf("\n");
 }
 err=sl_householder_aff(a,b,n);
 if(err==1)
 {
  printf("Solution :\n");
  for(j=1;j<=n;j++)
  printf("x%d = %22.16e\n",j,b[j]);
 }
 else printf("Erreur, matrice singulière\n");
}

			 A(1)                               b(1)
 8.00000e+00 -4.00000e+00  3.00000e+00  7.00000e+00  I   1.20000e+01
 4.00000e+00  2.00000e+00 -6.00000e+00  4.00000e+00  I   1.00000e+00
-1.60000e+01  6.00000e+00 -2.00000e+00 -1.50000e+01  I  -1.90000e+01
 6.00000e+00  1.00000e+01 -1.50000e+01  1.00000e+01  I   1.00000e+00
			 A(2)                               b(2)
-1.92873e+01  3.11086e+00  3.00716e+00 -1.92873e+01  I  -2.12575e+01
 0.00000e+00  3.04237e+00 -5.99895e+00  1.46588e-01  I  -3.87516e+00
 0.00000e+00  1.83053e+00 -2.00420e+00  4.13647e-01  I   5.00651e-01
 0.00000e+00  1.15636e+01 -1.49984e+01  4.21988e+00  I  -6.31274e+00
			 A(3)                               b(3)
-1.92873e+01  3.11086e+00  3.00716e+00 -1.92873e+01  I  -2.12575e+01
 0.00000e+00 -1.20964e+01  1.61499e+01 -4.13347e+00  I   6.93356e+00
 0.00000e+00  0.00000e+00  6.73959e-01 -1.03883e-01  I   1.80760e+00
 0.00000e+00  0.00000e+00  1.91966e+00  9.50616e-01  I   1.94336e+00
			 A(4)                               b(4)
-1.92873e+01  3.11086e+00  3.00716e+00 -1.92873e+01  I  -2.12575e+01
 0.00000e+00 -1.20964e+01  1.61499e+01 -4.13347e+00  I   6.93356e+00
 0.00000e+00  0.00000e+00 -2.03453e+00 -8.62532e-01  I  -2.43243e+00
 0.00000e+00  0.00000e+00  0.00000e+00  4.12918e-01  I  -1.06179e+00
Solution :
x1 =  4.571428571428572e+00
x2 =  3.357142857142857e+00
x3 =  2.285714285714285e+00
x4 = -2.571428571428572e+00

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

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

a : المصفوفة
b : العنصر الثاني
n : درجة المصفوفة

القيمة الرجعية

ترجع القيمة 1 عندما تكون المصفوفة مقلوبة رقميا والقيمة 0 إذا كان عكس ذلك.
إذا كانت قيمة الرجوع هي 1 فسيكون الحل في الجدول b.
خلال كل كرة يتم إظهار المصفوفة Ak والعنصر الثاني bk.
تساوي الثابتة NMAX البعد القصوي للمصفوفة زائد 1، وتساوي الثابتة N2MAX القيمة NMAX+1

int sl_householder_aff(double a[NMAX][NMAX],double b[NMAX],int n)
{
 int i,j,k,err,sgn;
 double s,nu;
 double aa[NMAX][NMAX];
 double bb[NMAX],v[NMAX];
 double gam[NMAX1];
 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;
 s=0;
 for(i=k;i<=n;i++)
 s+=v[i]*b[i];
 gam[n+1]=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];
 bb[i]=b[i]-gam[n+1]*v[i];
 }
 for(i=k;i<=n;i++)
 {
 for(j=k;j<=n;j++)
 a[i][j]=aa[i][j];
 b[i]=bb[i];
 }
 k++;
 printf(" A(%d) b(%d)\n",k,k);
 for(i=1;i<=n;i++)
 {
 for(j=1;j<=n;j++)
 printf("%12.6e ",a[i][j]);
 printf(" I %12.6e",b[i]);
 printf("\n");
 }
 }
 }
 if(fabs(a[n][n])<1e-15)err=0;
 if(err==1)
 {
 b[n]=b[n]/a[n][n];
 for(i=n-1;i>0;i--)
 {
 s=b[i];
 for(j=i+1;j<=n;j++)
 s-=a[i][j]*b[j];
 b[i]=s/a[i][i];
 }
 }
 return(err);
}