المبدأ

• توجد مصفوفة L^(k)  في كل مرحلة  من طريقة جوص بحيث ⁽A^(k+1) = L^(k)A^(k.

المصفوفة U = A⁽ⁿ⁾  مثلثية علوية و بالتالي  A = L^{(1) -1}L^{(2) -1} ... L^{(n-1) -1}U = LU .

المصفوفة L مثلثية سفلية للقطر الواحدي.

•  طريقة التفكيك LU تتجلى في إنشاء المصفوفتين L و U، ثم حل النظمتين المثلثيتين  Ly = b  و Ux = y  على التوالي.

الخوارزم

المرحلة 1: حساب L و U بحيث A = LU، مع L مثلثية سفلية للقطر الواحدي و U مثلثية علوية.

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

.U = A⁽ⁿ⁾

المرحلة  2: حل النظمة المثلثية   L y = b

المرحلة  3: حل النظمة المثلثية  U x = y

مثال

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

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

المرحلة 1: حساب L وU

المرحلة 2: لنحل النظمة المثلثية L y = b

المرحلة 3: لنحل النظمة المثلثية U x = y

البرمجة

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

int sl_decomp_lu(double a[NMAX][NMAX],double l[NMAX][NMAX],double u[NMAX][NMAX],int n);
void sl_resol_lu(double l[NMAX][NMAX],double u[NMAX][NMAX],double b[NMAX],int n);

main()
{
 double a[NMAX][NMAX],l[NMAX][NMAX],u[NMAX][NMAX],b[NMAX];
 int i,j,n,err;
 clrscr();
 n=4;
 printf("Méthode de la Décomposition LU\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 b\n");
 for(i=1;i<=n;i++)
 {
 for(j=1;j<=n;j++)
 printf("%12.6e ",a[i][j]);
 printf("I %12.6e \n",b[i]);
 }
 err=sl_decomp_lu(a,l,u,n);
 if(err==1)
 {
 printf(" L \n");
 for(i=1;i<=n;i++)
 {
 for(j=1;j<i;j++)
 printf("%12.6e ",l[i][j]);
 printf("%12.6e ",1.0);
 for(j=i+1;j<=n;j++)
 printf("%12.6e ",0.0);
 printf("\n");
 }
 printf(" U \n");
 for(i=1;i<=n;i++)
 {
 for(j=1;j<i;j++)
 printf("%12.6e ",0.0);
 for(j=i;j<=n;j++)
 printf("%12.6e ",u[i][j]);
 printf("\n");
 }
 sl_resol_lu(l,u,b,n);
 printf("Solution:\n");
 for(i=1;i<=n;i++)
 printf("x%d =%23.16e \n",i,b[i]);
 }
 else printf("Erreur, pivot nul");
}

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

يظهر البرنامج المصفوفة A والعنصر الثاني b للنظمة، وأيضا المصفوفتين L وU:

    A                                    b
 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  
   L            
 1.00000e+00   0.00000e+00   0.00000e+00   0.00000e+00  
 5.00000e-01   1.00000e+00   0.00000e+00   0.00000e+00  
-2.00000e+00  -5.00000e-01   1.00000e+00   0.00000e+00  
 7.50000e-01   3.25000e+00   2.85000e+01   1.00000e+00  
   U            
 8.00000e+00  -4.00000e+00   3.00000e+00   7.00000e+00  
 0.00000e+00   4.00000e+00  -7.50000e+00   5.00000e-01  
 0.00000e+00   0.00000e+00   2.50000e-01  -7.50000e-01  
 0.00000e+00   0.00000e+00   0.00000e+00   2.45000e+01  
Solution:
x1 =  4.571428571428571e+00  
x2 =  3.357142857142856e+00  
x3 =  2.285714285714285e+00  
x4 = -2.571428571428572e+00  

• ---

  برمجة الدالة الرئيسية sl_decomp_lu

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

الدالة sl_decomp_lu
a : المصفوفة
l و u : مصفوفتين (يمثلان L و U) 
n : درجة المصفوفة

ترجع القيمة 1 عندما لا يكون pivot غير منعدم والقيمة 0 إذا كان عكس ذلك.
إذا كانت قيمة الرجوع هي 1 فسيكون الحل في الجدولين l و u.

ملاحظة: في حالة لم تستطع الدالة حساب التفكيك LU . فيجب في هذه الحالة تغيير pivot وبالتالي فسيتم استخدام الدالتين ( sl_decomp_lu_pivmax و sl_resol_lu_pivmax).

• l و u : المصفوفتين L و U
• b: العنصر الثاني
• n : درجة المصفوفة

ترجع بالحل ممثلا في الجدول b.

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

int sl_decomp_lu(double a[NMAX][NMAX],double l[NMAX][NMAX],double u[NMAX][NMAX],int n)
{
 int i,j,k,err;
 double pivot,coef,s;
 for(i=1;i<=n;i++)
 {
 for(j=1;j<=n;j++)
 u[i][j]=a[i][j];
 }
 err=1;
 k=1;
 while (err==1 && k<n)
 {
 pivot=u[k][k];
 if(pivot!=0)
 for(i=k+1;i<=n;i++)
 {
 coef=u[i][k]/pivot;
 for(j=k+1;j<=n;j++)
 u[i][j] -= coef*u[k][j];
 l[i][k]=coef;
 }
 else err=0;
 k++;
 }
 if(u[n][n]==0) err=0;
 return(err);
}

void sl_resol_lu(double l[NMAX][NMAX],double u[NMAX][NMAX],double b[NMAX],int n)
{
 int i,j;
 double s;
 double y[NMAX];
 y[1]=b[1];
 for(i=2;i<=n;i++)
 {
 s=b[i];
 for(j=1;j<i;j++)
 s-=l[i][j]*y[j];
 y[i]=s;
 }
 b[n]=y[n]/u[n][n];
 for(i=n-1;i>=1;i--)
 {
 s=y[i];
 for(j=i+1;j<=n;j++)
 s-=u[i][j]*b[j];
 b[i]=s/u[i][i];
 }
}