Оригинальный DVD-ROM: eXeL@B DVD !
eXeL@B ВИДЕОКУРС !

ВИДЕОКУРС ВЗЛОМ
выпущен 2 июня!


УЗНАТЬ БОЛЬШЕ >>
Домой | Статьи | RAR-cтатьи | FAQ | Форум | Скачать | Видеокурс
Новичку | Ссылки | Программирование | Интервью | Архив | Связь

ПРОГРАММИРОВАНИЕ НА C и С++



Слушай, дружище, зачем так мучиться с этим языком С++, ты ведь не Билл Гейтс. Возьми тот же Python и программируй, он кроссплатформенный, под Windows тоже работает. Я сам давно заметил: то что на Си пишешь в страницу кода, на питоне решается в одну-две строки. При том, питон намного проще, я его сам недавно изучил по видеокурсу вот этому. Кстати, автор отлично там объясняет. Буквально день-два и уже будешь писать на нём, чего не скажешь про сложный С++.

Нижне-верхняя (LU) декомпозиция матрицы (алгоритм Краута)

Первым из алгоритмов, посвященным большому разделу решения систем линейных уравнений, представляем алгоритм Краута. Это фактически метод решения систем общего вида, конкурирующий по быстродействию с общеизвестным методом Гаусса-Жордана, но позволяющий более эффективно использовать решение.

Если мы можем разложить матрицу линейной системы A в произведение A=L*U, где L - нижняя, а U - верхняя треугольные матрицы, то решение системы уравнений с произвольной правой частью производится весьма просто, применением двух обратных подстановок. Более того, в отличие от известного метода Гаусса-Жордана, разложенная матрица позволяет быстро решать серии линейных уравнений с различными правыми частями при одной и той же матрице.

Метод Краута позволяет провести LU-декомпозицию матрицы примерно за то же число операций, что и "прямая" часть метода Гаусса-Жордана. Итоговые коэффициенты двух треугольных матриц упаковываются в матрицу того же размера, что и A, и на том же месте в памяти; при этом верхняя матрица U размещается в наддиагональной части и на диагонали; нижняя L в поддиагональной части, а диагональные элементы L считаются все равными 1 (без ограничения общности) и не выводятся.

Алгоритм Краута представляет из себя следующее:

  1. Положить Lii=1 i=0,...,N-1 (здесь ничего не производится, а только имеется в виду для дальнейших шагов;
  2. Для каждого j=0,1,...,N-1 проделать:
    1. Для всех i=0,...,j вычислить Uij:
      Uij=Aij - SUM(0<=k<i)(Lik*Uik)
      (при i=0 сумму полагать нулем);
    2. Для всех i=j+1,...,N-1 вычислить:
      Lij = Aij - SUM(0<=k<j)(Lik*Uik))/Uii
      Обе процедуры выполняются до перехода к новому j.

Устойчивая работа алгоритма Краута возможно только при условии выбора ведущего элемента для каждого обращения к j из п.2 алгоритма. Выбор производится перед выполнением п.2 для очередного j перестановки необработанных строк матрицы так, чтобы строка, подставляемая на место j, содержала наибольший по модулю элемент в колонке j. Порядок перестановки необходимо запомнить, для чего достаточно использовать дополнительный вектор целых величин. Еще один вектор используется внутри алгоритма для масштабирования.

Алгоритм Краута имеет несколько приложений, одно из которых - решение системы линейных уравнений (обратная подстановка с учетом порядка перестановки строк), другие - вычисление обратной матрицы и вычисление детерминанта. Подробнее вызовы функций, связанных с алгоритмом Краута и реализация деталей алгоритма находятся в комментариях к программе и в ней самой.

Программа самой декомпозиции и некоторых ее приложений приводится ниже



 
 /*
 
   LU-decomposition according to Crout's algorithm with pivoting.
 
 	Description:
 
    int LU_decompos(double **a,int n,int *indx,int *d,double *vv);
 
 	Parameters:
 
    a - source matrix (n x n) on input, destination on output;
 
    n - the matrix size;
 
    indx - integer array (size n) to remember permutations;
 
    d - on output, contains +1 or -1 for even or odd permutations number.
 
    vv - temporary array (size n).
 
 	   Returns:
 
    0 - the source matrix is singular (invalid for decomposition),
 
    1 - if OK.
 
 
 
   Back substitution, using LU decomposed matrix.
 
 	Description:
 
   void LU_backsub(double **a,int n,int *indx,double *b);
 
 	Parameters:
 
   a - the matrix decomposed by Crout;
 
   n - the matrix size;
 
   indx - permutation order obtained by decomposition algorithm;
 
   b - the vector (size n) to be substituted on input, the result
 
       of the substitution on output.
 
   Note: a and indx are not modified by this routine and could be
 
   used in multiple calls.
 
 
 
   Invertation of matrix, using LU decomposed matrix.
 
 	Description:
 
   void LU_invert(double **a,int n,int *indx,double **inv,double *col);
 
 	Parameters:
 
   a - the matrix decomposed by Crout;
 
   n - the matrix size;
 
   indx - permutation order obtained by decomposition algorithm;
 
   inv - the destination matrix;
 
   col - temporary array (size n).
 
   Note: test for singularity has been already obtained on the
 
   matrix decomposition, a and indx are not modified by this routine,
 
   the routine uses multiple backsubstitutions (previous algorithm).
 
 
 
   Obtaining the matrix determinant, using LU-decomposed matrix
 
 	Description:
 
   double LU_determ(double **a,int n,int *indx,int *d);
 
 	Parameters:
 
   a - the matrix decomposed by Crout;
 
   n - the matrix size;
 
   indx - permutation order obtained by decomposition algorithm;
 
   d - the parity sign (+1 or -1) obtained at decomposition.
 
 	Returns:
 
   the determinant value. Note: non-zero (the matrix cannot be
 
   singular, if decomposed properly); a, indx and d are not modified
 
   by this routine.
 
 
 
 */
 
 
 
 /* for fabs(); inclusion could be removed if fabs() is implemented inline */
 
 #include <math.h>
 
 /* "small number" to avoid overflow in some cases */
 
 #define TINY 1.e-30
 
 
 
 /* the decomposition itself */
 
 int LU_decompos(double **a, int n, int *indx, int *d, double *vv) {
 
   register int i,imax,j,k;
 
   double big,sum,temp;
 
   /* set 1 for initially even (0) transmutations number */
 
   *d=1;
 
   /* search for the largest element in each row; save the scaling in the
 
      temporary array vv and return zero if the matrix is singular */
 
   for(i=0;i<n;i++) {
 
     big=0.;
 
     for(j=0;j<n;j++) if((temp=fabs(a[i][j]))>big) big=temp;
 
     if(big==0.) return(0);
 
     vv[i]=big;
 
   }
 
   /* the main loop for the Crout's algorithm */
 
   for(j=0;j<n;j++) {
 
     /* this is the part a) of the algorithm except for i==j */
 
     for(i=0;i<j;i++) {
 
       sum=a[i][j];
 
       for(k=0;k<i;k++) sum-=a[i][k]*a[k][j];
 
       a[i][j]=sum;
 
     }
 
     /* initialize for the search for the largest pivot element */
 
     big=0.; imax=j;
 
     /* this is the part a) for i==j and part b) for i>j
 
        plus pivot search */
 
     for(i=j;i<n;i++) {
 
       sum=a[i][j];
 
       for(k=0;k<j;k++) sum-=a[i][k]*a[k][j];
 
       a[i][j]=sum;
 
       /* is the figure of merit for the pivot better than
 
          the best so far? */
 
       if((temp=vv[i]*fabs(sum))>=big) {big=temp;imax=i;}
 
     }
 
     /* interchange rows, if needed, change parity and the scale factor */
 
     if(imax!=j) {
 
       for(k=0;k<n;k++) {
 
         temp=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=temp;
 
       }
 
       *d=-(*d); vv[imax]=vv[j];
 
     }
 
     /* store the index */
 
     indx[j]=imax;
 
     /* if the pivot element is zero, the matrix is singular but for some
 
        applications a tiny number is desirable instead */
 
     if(a[j][j]==0.) a[j][j]=TINY;
 
     /* finally, divide by the pivot element */
 
     if(j<n-1) {
 
       temp=1./a[j][j];
 
       for(i=j+1;i<n;i++) a[i][j]*=temp;
 
     }
 
   }
 
   return(1);
 
 }
 
 
 
 /* the back substitution */
 
 void LU_backsub(double **a,int n,int *indx,double *b) {
 
   register int i,j,ip,ii=-1;
 
   double sum;
 
   /* First step of backsubstitution; the only wrinkle is to unscramble
 
      the permutation order. Note: the algorithm is optimized for a
 
      possibility of large amount of zeroes in b */
 
   for(i=0;i<n;i++) {
 
     ip=indx[i];
 
     sum=b[ip];b[ip]=b[i];
 
     if(ii>=0) for(j=ii;j<i;j++) sum-=a[i][j]*b[j];
 
     else if(sum) ii=i;	/* a nonzero element encounted */
 
     b[i]=sum;
 
   }
 
   /* the second step */
 
   for(i=n-1;i>=0;i--) {
 
     sum=b[i];
 
     for(j=i+1;j<n;j++) sum-=a[i][j]*b[j];
 
     b[i]=sum/a[i][i];
 
   }
 
 }
 
 
 
 /* the matrix invertation */
 
 void LU_invert(double **a,int n,int *indx,double **inv,double *col) {
 
   register int i,j;
 
   for(j=0;j<n;j++) {
 
     for(i=0;i<n;i++) col[i]=0.;
 
     col[j]=1.;
 
     LU_backsub(a,n,indx,col);
 
     for(i=0;i<n;i++) inv[i][j]=col[i];
 
   }
 
 }
 
 
 
 /* calculating determinant */
 
 double LU_determ(double **a,int n,int *indx,int *d) {
 
   register int j;
 
   double res=(double)(*d);
 
   for(j=0;j<n;j++) res*=a[j][j];
 
   return(res);
 
 }
 
 

Здесь Вы можете скачать текст программы crout.c

<< ВЕРНУТЬСЯ В ПОДРАЗДЕЛ

<< ВЕРНУТЬСЯ В ОГЛАВЛЕНИЕ




Материалы находятся на сайте https://exelab.ru/pro/



Оригинальный DVD-ROM: eXeL@B DVD !


Вы находитесь на EXELAB.rU
Проект ReactOS