Полезное для программистов:

Фриланс
Новости
Статьи
   
Рубрики:


Решение систем линейных уравнений

Поиск:
Метод квадратного корня

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (square root method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    From given matrix A we build two auxulary triangular matrixes S & D
//    Then solving equations Gy=b (where g[I][j]=s[j][I]*d[j]) and Sx=y
//    we obtain vector x.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>

#define N    4     // size of matrix
#define N1   N+1

float matrix[N][N1]=
    {{10.9,1.2,2.1,0.9,23.2},
  {1.2,11.2,1.5,2.5,38.1},
  {2.7,1.5,9.8,1.3,40.3},
  {0.9,2.5,1.3,12.1,58.2}
    };

void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%fn",matrix[i][N]);
  }
}

void main(void)
{
  // Variables declaration
  register char i,j,k;
  float s[N][N],x[N],y[N],d[N],sum;

  // Printing given matrix
  ShowMatrix();

  // Building matrixes S and D
  for (j=0;j<N;j++)
    for (i=0;i<=j;i++)
    {
      sum=matrix[i][j];
      for (k=0;k<i;k++)
        sum-=s[k][i]*d[k]*s[k][j];
      if (i==j)
      {
        d[i]=(sum>0?1:-1);
        s[i][j]=sqrt(fabs(sum));
      }
      else s[i][j]=sum/(d[i]*s[i][i]);
    }

  // Solving equation Gy=b (G: g[I][j]=s[j][I]*d[j])
  for (i=0;i<N;i++)
  {
    y[i]=matrix[i][N];
    for (j=0;j<i;j++) y[i]-=s[j][i]*d[j]*y[j];
    y[i]/=s[i][i]*d[i];
  }

  // Solving equation Sx=y
  for (i=N-1;i>=0;i--)
  {
    x[i]=y[i];
    for (j=i+1;j<N;j++) x[i]-=s[i][j]*x[j];
    x[i]/=s[i][i];
  }

  // Printing solution
  printf("nSolution:n");
  for (i=0;i<N;i++)
  printf("x%d=%fn",i+1,x[i]);
}



Метод Некрасова

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Nekrasov method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This is an iterative method of solving systems of linear equations.
//    Criteria for end of iterations is ||xn||-||x(n-1)||<epsilon, where
//    ||x|| is norm of vector x.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>

#define N     4     // size of matrix
#define N1    N+1

float matrix[N][N1]=
      {{10.9,1.2,2.1,0.9,23.2},
       {1.2,11.2,1.5,2.5,38.1},
       {2.7,1.5,9.8,1.3,40.3},
       {0.9,2.5,1.3,12.1,58.2}
     };

void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%fn",matrix[i][N]);
  }
}

// this function calculates ||x||
float norma(float *x)
{
  float sum=0;

  for (unsigned i=0;i<N;i++) sum+=fabs(*(x+i));
  return sum;
}

void main(void)
{
  // Variables declaration
  register char i,j;
  float x_current[N],x_next[N],sum1,sum2,epsilon;

  epsilon=1e-7;
  // Printing given matrix
  ShowMatrix();

  // Solving
  do
  {
    for (i=0;i<N;i++) x_current[i]=x_next[i];
    for (i=0;i<N;i++)
    {
      sum1=0;sum2=0;
      for(j=0;j<i;j++) sum1+=matrix[i][j]*x_next[j];
      for(j=i+1;j<N;j++) sum2+=matrix[i][j]*x_current[j];
      x_next[i]=(matrix[i][N]-sum1-sum2)/matrix[i][i];
    }
  } while (fabs(norma(x_current)-norma(x_next))>epsilon);

  // Printing solution
  printf("nSolution:n");
  for (i=0;i<N;i++)
    printf("x%d=%fn",i+1,x_current[i]);
}



Метод 3-диагональных матриц

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (3-diagonal matrix)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This method is developed for 3-diagonal matrixes:
//    (a11 a12 0 ..............)
//    (a21 a22 a23 0 ..........)
//    (0   a32 a33 a34 0 ......) = M
//    (........................)
//    (.......... 0 an(n-1) ann)
//    Cause all elements except three diagonals is 0 we'll store
//    only non zero elements in the A,B,C arrays. Array D stores
//    vector B (M*X=B). We build two auxulary vectors P and Q and
//    then calculate vector X using the following recurrent formula:
//      X(k-1)=P(k)*X(k)+Q(k)
//
//////////////////////////////////////////////////////////////////////////////


#include <stdio.h>

#define N    20
#define N1   N+1

void main(void)
{
  // Variables declaration
  register int i,i1,k;  // using processor registers for better performance
  float a[N]={0,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1};
  float b[N]={2,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,4,2};
  float c[N]={-1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0};
  float d[N]={0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0};
  float x[N1],p[N1],q[N1],tmp;

  // Some initializations

  p[0]=0;
  q[0]=0;
  x[N]=0;
  c[N-1]=0; // to avoid overflow

  // Building vectors P and Q
  for (i=0;i<N;i++)
  {
    i1=i+1;
    if ((tmp=b[i]+a[i]*p[i])==0)
    {
      printf("Wrong method for such system...n");
      return;
    };
    p[i1]=-c[i]/tmp;
    q[i1]=(d[i]-a[i]*q[i])/tmp;
  }

  // Building vector X
  for (k=N;k>0;k--) x[k-1]=p[k]*x[k]+q[k];

  // Printing solution
  printf("Solution:n");
  for (i=0;i<N;i++)
    printf("x%d=%fn",i+1,x[i]);
}



Метод Гаусса

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Gauss method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This methos based on excluding variables from equations.
//    Using first equation we exclude x1 from all other equations
//    Using second equation equation we exclude x2 from all equations
//   (excluding first) and so on. So we have triangular matrix from which
//    we can easily get vector X. (AX=B)
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>

#define N     4     // size of matrix
#define N1    N+1

float matrix[N][N1]=
                  {{2,5,4,1,28},
                   {1,3,2,1,17},
                   {2,10,9,7,77},
                   {3,8,9,2,54}
                  };

void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%fn",matrix[i][N]);
  }
}

void main(void)
{
  // Variables declaration
  float tmp,x[N1];
  register short int i,j,k;

  // Printing given matrix
  ShowMatrix();
  // Main loop
  for (i=0;i<N;i++)
  {
    // Excluding variable x[i] from equations
    tmp=matrix[i][i];
    for (j=N;j>=i;j--) matrix[i][j]/=tmp;
    for (j=i+1;j<N;j++)
    {
      tmp=matrix[j][i];
      for (k=N;k>=i;k--)
        matrix[j][k]-=tmp*matrix[i][k];
    }
  }

  // Calculating vector x
  x[N-1]=matrix[N-1][N];
  for (i=N-2;i>=0;i--)
  {
    x[i]=matrix[i][N];
    for (j=i+1;j<N;j++) x[i]-=matrix[i][j]*x[j];
  }

  // Printing solution
  printf("nSolution:n");
  for (i=0;i<N;i++)
    printf("x%d=%fn",i+1,x[i]);
}



Метод Гаусса с выбором максимального элемента

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Gauss method with max element selecting)
//  (c) Johna Smith, 1996
//
//  Method description:
//    This methos based on excluding variables from equations.
//    Using first equation we exclude x1 from all other equations
//    Using second equation equation we exclude x2 from all equations
//    (excluding first) and so on. So we have triangular matrix from which
//    we can easily get vector X. (AX=B)
//    'Selecting max element' means that on each step we select equation
//    with max element in current column and use it as next equation to exclude
//    variable.
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>

#define N      4     // size of matrix
#define N1     N+1

float matrix[N][N1]=
    {{2,5,4,1,28},
  {1,3,2,1,17},
   {2,10,9,7,77},
  {3,8,9,2,54}
        };

void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%fn",matrix[i][N]);
  }
}

void main(void)
{
  // Variables declaration
  float max,swp,tmp,x[N1];
  register short int row_with_max_element,i,j,k;

  // Printing given matrix
  ShowMatrix();
  // Main loop
  for (i=0;i<N;i++)
  {
    // Searching for max element in the current column (i)
    max=fabs(matrix[i][i]);
    row_with_max_element=i;
    for (j=i+1;j<N;j++)
      if (fabs(matrix[j][i])>max)
      {
        row_with_max_element=j;
        max=fabs(matrix[j][i]);
      }
    // Swapping 2 lines of matrix - row_with_max_element & i
    if (row_with_max_element!=i)
    for (j=0;j<N1;j++)
    {
      swp=matrix[row_with_max_element][j];
      matrix[row_with_max_element][j]=matrix[i][j];
      matrix[i][j]=swp;
    }
    // Excluding variable x[i] from equations
    tmp=matrix[i][i];
    for (j=N;j>=i;j--) matrix[i][j]/=tmp;
    for (j=i+1;j<N;j++)
    {
      tmp=matrix[j][i];
      for (k=N;k>=i;k--)
        matrix[j][k]-=tmp*matrix[i][k];
    }
  }

  // Calculating vector x
  x[N-1]=matrix[N-1][N];
  for (i=N-2;i>=0;i--)
  {
    x[i]=matrix[i][N];
    for (j=i+1;j<N;j++) x[i]-=matrix[i][j]*x[j];
  }

  // Printing solution
  printf("nSolution:n");
  for (i=0;i<N;i++)
    printf("x%d=%fn",i+1,x[i]);
}



Метод итераций

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (iterations method)
//  (c) Johna Smith, 1996
//
//  Method description:
//  This program automatically transforms system of equations to
//  iterative form and solves it.
//  Iterative form:
//     x1=a11*x1+a12*x2+... + a1n*xn+b1
//     x2=a21*x1+a22*x2+... + a2n*xn+b2
//     ...
//     xn=an1*x1+an2*x2+... + ann*xn+bn
//  Setting first iteration of vector X we can get next iteration from
//  such form of system of equations. Iterations will be finished when
//  difference between two consequitive iterations will be less than epsilon
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>

#define N      3     // size of matrix
#define N1     N+1

float matrix[N][N1]=
    {{14.38,-2.41,1.39,5.86},
  {1.84,25.36,-3.31,-2.28},
  {2.46,-3.49,16.37,4.47}
    };
int maxiterations=10;  // maximum number of iterations
float epsilon=0.0001;  // required accuracy

void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%fn",matrix[i][N]);
  }
}

void main(void)
{
  // Variables declaration
  float x[N],y[N],t;
  register short int i,j,k;
  int iterations=0;

  // Printing given matrix
  ShowMatrix();

  // setting first iteration of vector X
  for (i=0;i<N;i++) x[i]=matrix[i][N];
  do
  {
    for (i=0;i<N;i++)
    {
      t=-matrix[i][N];
      for (j=0;j<N;j++)
        t+=matrix[i][j]*x[j];
        y[i]=(-t+matrix[i][i]*x[i])/matrix[i][i];
    }
    iterations++;
    k=0;
    // checking solution
    while (fabs(x[k]-y[k])<epsilon && k<N) k++;
    // new iteration becomes old
    for(i=0;i<N;i++) x[i]=y[i];
  } while (k!=N && iterations<maxiterations);

  if (iterations==maxiterations)
  {
    printf("Iterations are very slow...");
  } else
  {
    // Printing solution
    printf("nSolution:n");
    for (i=0;i<N;i++)
      printf("x%d=%fn",i+1,x[i]);
    printf("%d iterations were made",iterations);
  }
}



Метод Монте-Карло

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Monte-Karlo method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    We should translate given system of linear equations to special form:
//      x1=a(11)x1+a(12)x2+...+a(1n)xn+a(1n+1)
//      x2=a(21)x2+a(22)x2+...+a(2n)xn+a(2n+1)
//      ....
//      xn=a(n1)x1+a(n2)x2+...+a(nn)xn+a(nn+1)
//            n
//     where SUM |aij|<1 (i=1, 2, ... ,n)
//           j=1
//    Lets divide interval [0;1] into N+1 smaller intervals.
//    Imagine a particle that is moving randomly at [0;1] interval.
//    Remember its moves until it reachs one of the bounds of the interval.
//    So we'll get a trajectory of the particle s(i),s(j),s(k),s(l),s(m),...
//    If particle moves from interval s(i) to s(j) we'll write v(ij)
//    y[i]=v(ij)*v(jk)*...*v(tm)*w(m), v(ij)=sign(aij), v(jk)=sign(ajk)
//                    n
//   w(m)=a(mn+1)/(1-SUM a(mj)),   x[i]=y/M, where M is number of particle runs
//                   j=1
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define N     3     // size of matrix
#define N1    N+1
#define M     1000 // number of tries

float matrix[N][N1]=
    {{3,1,-1,-2},
  {1,4,-2,-3},
  {2,-1,5,-15}
    };

void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    printf("x%d=",i+1);
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],j+1);
    printf("%+fn",matrix[i][N]);
  }
}

void main(void)
{
  // Variables declaration
  float b[N][N],w[N],y,c,x[N],l;
  register short int i,j,v;
  unsigned char finish;
  int row,   // currently analysed row
      tries; // number of tries to calculate true root by generating randoms

  // Transforming matrix
  for (i=0;i<N;i++)
  {
    l=0;
    for (j=0;j<N1;j++) l+=fabs(matrix[i][j]);
    v=(matrix[i][i]>0?-1:1);
    for (j=0;j<N1;j++) matrix[i][j]=v*(matrix[i][j])/l+(i==j?1:0);
  }
  // Printing transformed matrix
  ShowMatrix();

  // filling matrix B
  for (i=0;i<N;i++)
  {
    b[i][0]=fabs(matrix[i][0]);
    for(j=1;j<N;j++)
      b[i][j]=b[i][j-1]+fabs(matrix[i][j]);
  }
  // filling matrix w
  for(i=0;i<N;i++) w[i]=matrix[i][N]/(1-b[i][N-1]);

  // for each equation in the system
  for(i=0;i<N;i++)
  {
    tries=0;
    row=i;
    y=0;
    v=1;
    while (tries<M)
    {
      // generating random number
      c=random(32767)/32767.0;
      j=N-1;
      finish=0;
      while(j>=0 && !finish)
      {
        if (c>b[row][j])
        {
          if (j==N-1)
          {
            tries++;
            y+=v*w[row];
            row=i;v=1;
          }
          else
          {
            v*=(matrix[row][j+1]>=0?1:-1);
            row=j+1;
          }
          finish=1;
        }
        j--;
      }
      if (!finish)
      {
        v*=(matrix[row][j+1]>=0?1:-1);
        row=0;
      }
    }
    // calculating root
    x[i]=y/M;
  }

  // Printing solution
  printf("nSolution:n");
  for (i=0;i<N;i++)
    printf("x%d=%fn",i+1,x[i]);
}



Метод ортогонализации

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (orthogonalization method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given system of linear equations:
//      a11 x1 + a12 x2 + ... + a1n xn + a(1n+1) = 0
//      a21 x1 + a22 x2 + ... + a2n xn + a(2n+1) = 0
//
//      an1 x1 + an2 x2 + ... + ann xn + a(nn+1) = 0
//
//    Left part of each equation is a result of scalar multiplication of
//    two vectors: ai=(ai1,ai2,...,ain,a(in+1)) and x=(x1,x2,..,xn,1)
//    So to solve this system we need only to build vector x which will be
//    orthogonal to each of ai vectors.
//                            n+1
//    Let u1=a1,   b1=u1/sqrt(SUM u(1j)^2)
//                            j=1
//
//    Other ui we can get from the following iterative formula:
//                      k
//       u(k+1)=u(k+1)-SUM {u(k+1), bj}bj, where {..} means scalar multiplication
//                     j=1
//                          n+1
//       b(k+1)=u(k+1)/sqrt(SUM u(k+1,j)^2)
//                          j=1
//    And finally we can obtain roots of this system from this formula:
//
//    xi=b(n+1,i)/b(n+1n+1)
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define N     3     // size of matrix
#define N1    N+1

double matrix[N1][N1]=
    {{3,1,-1,-2},
  {1,4,-2,-3},
  {2,-1,5,-15}
    };

void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("%+f=0n",matrix[i][N]);
  }
}

void main(void)
{
  // Variables declaration
  double c[N],x[N],s;
  int i,j,k,l,m;

  ShowMatrix();

  // expanding system by vector (0,0,0,...,0,1)
  for (i=0;i<N;i++) matrix[N][i]=0;
  matrix[N][N]=1;

  // for all equations
  for (i=0;i<N1;i++)
  {
    if (i>0)
    {
      k=0;
      // make some iterations to increase accuracy of calculations
      while (k<=3)
      {
        for (l=0;l<i;l++)
        {
          c[l]=0;
          for(m=0;m<N1;m++) c[l]+=matrix[l][m]*matrix[i][m];
        }
        for (j=0;j<N1;j++)
        {
          s=0;
          for(l=0;l<i;l++) s+=c[l]*matrix[l][j];
          matrix[i][j]-=s;
        }
        k++;
      }
    }
     // normalizing vector
     s=0;
     for (k=0;k<N1;k++) s+=matrix[i][k]*matrix[i][k];
     s=sqrt(s);
     for (j=0;j<N1;j++) matrix[i][j]/=s;
  }

  s=matrix[N][N];
  for (j=0;j<N;j++) x[j]=matrix[N][j]/s;

  // Printing solution
  printf("nSolution:n");
  for (i=0;i<N;i++)
    printf("x%d=%fn",i+1,x[i]);
}



Метод Зайделя

Код

//////////////////////////////////////////////////////////////////////////////
//
//  Solving system of linear equations (Zaidel method)
//  (c) Johna Smith, 1996
//
//  Method description:
//    Given: system of linear equations
//       a11 x1 + a12 x2 + ... + a1n xn = b1
//       a21 x1 + a22 x2 + ... + a2n xn = b2
//       ...
//       an1 x1 + an2 x2 + ... + ann xn = bn
//
//    We use following iterative formula to solve this system:
//
//                     1    i-1                N
//     xi(j+1)=xi(j)- --- ( SUM aik xk(j+1) + SUM aik xk(j) -bi)
//                    aii   k=1               k=i
//
//     where xi(j+1) is j+1-th iteration, xi(j) is j-th iteration
//
//////////////////////////////////////////////////////////////////////////////

#include <stdio.h>
#include <math.h>

#define N    3     // size of matrix
#define N1   N+1

float matrix[N][N1]=
    {{10,1,1,12},
  {2,10,1,13},
  {2,2,10,14}
    };
float z[N]={0,0,0};  // first iteration

float epsilon=0.0001;  // required accuracy

void ShowMatrix(void)
{
  for (int i=0;i<N;i++)
  {
    for (int j=0;j<N;j++)
      printf("%+f*x%d",matrix[i][j],i+1);
    printf("=%fn",matrix[i][N]);
  }
}

void main(void)
{
  // Variables declaration
  float x[N];
  short int i,j;
  int iterations=0,finish=0;

  // Printing given matrix
  ShowMatrix();

  while (!finish)
  {
    finish=1;
    for (i=0;i<N;i++)
    {
      x[i]=-matrix[i][N];
      for (j=0;j<N;j++) x[i]+=matrix[i][j]*z[j];
      // don't stop iterations until required accuracy is reached
      if (fabs(x[i]/matrix[i][i])>=epsilon) finish=0;
      x[i]=z[i]-x[i]/matrix[i][i];
      // next iteration
      z[i]=x[i];
    }
    iterations++;
  }

  // Printing solution
  printf("nSolution:n");
  for (i=0;i<N;i++)
    printf("x%d=%fn",i+1,x[i]);
  printf("%d iterations were made",iterations);
}
Сайт: forum.vingrad.ru






Просмотров: 2969

 

 

Новые статьи:


Популярные:
  1. Как сделать цикличным проигрывание MIDI-файла?
  2. Создание AVI файла из рисунков
  3. Как устройство "отключить в данной конфигурации"?
  4. Kто в данный момент присоединен через Сеть?
  5. Как узнать количество доступной памяти?
  6. Как реализовать в RichEdit разноцветный текст?
  7. Как скрыть свое приложение от ProcessViewer
  8. Как программно нажать/скрыть/показ кнопку "Start"?
  9. Модуль работы с ресурсами в PE файлах
10. Функции вызова диалоговых окон выбора
11. Проверка граматики средствами Word'а из Delphi.
12. Модуль для упрощенного вызова сообщений
13. Функции для записи и чтение своих данных в, ЕХЕ- файле
14. Рекурсивный просмотр директорий
15. Network Traffic Monitor
16. Разные модули
17. Универсальная функция для обращения к любым экспортируем функциям DLL
18. Библиотека от VladS
19. Протектор для UPX'а
20. Еще об ICQ, сообщения по контакт листу?
21. Использование открытых интерфейсов
22. Теория и практика использования RTTI
23. Работа с TApplication
24. Примеры использования Drag and Drop для различных визуальных компонентов
25. Что такое порт? Правила для работы с портами
26. Симфония на клавиатуре
27. Загрузка DLL
28. Исправление автоинкремента
29. Взаимодействие с чужими окнами
30. Проверить дубляжи в столбце


 

 

 
 
На главную