Где я? Кто я?
Профиль
Группа: Экс. модератор
Сообщений: 3094
Регистрация: 25.3.2002
Где: СПб
Репутация: 18 Всего: 62
|
Решение систем линейных уравненийМетод квадратного корня Код | ////////////////////////////////////////////////////////////////////////////// // // 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("=%f\n",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=%f\n",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("=%f\n",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=%f\n",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=%f\n",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("=%f\n",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=%f\n",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("=%f\n",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=%f\n",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("=%f\n",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=%f\n",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("%+f\n",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=%f\n",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=0\n",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=%f\n",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("=%f\n",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=%f\n",i+1,x[i]); printf("%d iterations were made",iterations); }
|
|