【計算方法】線形方程式群を解く4つの方法

1.ガウス列主消元法を用いて解く

#include
#include
#include
#include
/* 2018/4/18                  。                ，     。           ，         ，         ，       。 */ 
/*************************************** *              * 1.       * 2.    * 3.       * 4.  --       * *****************************************/ 
#define N 5
#define M 4
void dep(double *p);
void display(double array[][N]);



static double A[M][N] = { {8,7,0,0,0}, {6,12,5,0,-2}, {0,4,9,3,8}, {0,0,1,2,6} };
double B[N] = {0,-2,8,6};

//   
void swap(double *p, double *t) { 
    double s;
    s = *p;
    *p = *t;
    *t = s; 
}

//      max     
void f(double *p, int row) {

    double max = fabs(*p);

    int r1 = row; 
    int r2 = row;       //       
    int i = 1;

    for(;row < M - 1; row++, i++) { 

        if (max < fabs( *(p + i * N) ) ) {
            max = *(p + i * N);
            r1 = row;     //        
        }
    }   

    //       
    for(int i = 0; i < N; i++) {     
        swap(&A[r2][i], &A[r1][i]);  
    }

    dep(p);

    return ;
}

void display(double array[][N]) {

    for(int i = 0; i < M; i++) {    
        for(int j = 0; j < N; j++) {
            printf("%10f",array[i][j]);
        }
        printf("
");
    } 
    return ;
}

//   ----      
void dep(double* p) {
    double result;
    double *temp; 
    double *tail;
    //  row+1 ----N  row    
    for (int i = 1; i < M ; i++) {
        tail = p;
        result =  *(p + i * N) / *p;
        temp = p + i * N;
        for(int j = 0; j < N; j++) {    
            *temp -= *tail * result;    //       ，   0 
            temp++;
            tail++;
        }
    } 

    return ; 
}
/* //     double* solve() { double x[M]; int i; memset(x, 0, sizeof(x)); //      x[M - 1] = A[M - 1][N -1] / A[M - 1][N - 2]; for(i = M -2; i >= 0; i--){ double s = A[i][N - 1]; for (int j = i + 1; j < N -1; j++) { s -= (A[i][j]*x[i + 1]); } x[i] = s / A[i][i]; } return x; } void merge() { double newArray[N][N]; //         for(int i = 0; i < N ; i++ ) { for(int j = 0; j < N + 1; j++) { newArray[i][j] = A[i][j]; if(j == N) { newArray[i][j] = B[i]; } } } } */
int main() {     
/********************************************** *  1.       * 1）    * 2）   * 3）     ***********************************************/    
    double x[M];
    int i;
    memset(x, 0, sizeof(x));        //      

    display(A);
    for(i = 0; i < 4; i++) {
        f(&A[i][i],i);
    }
    printf("      
");
    system("pause");
    system("cls");
    display(A);
    //     
    x[M - 1] = A[M - 1][N -1] / A[M - 1][N - 2];

    for(i = M -2; i >= 0; i--){

        double s = A[i][N - 1];
        for (int j = i + 1; j < N -1; j++) {
            s -= (A[i][j]*x[i + 1]);
        }
        x[i] = s / A[i][i];
    }
    for (i = 0; i < M; i++) {
        printf("%-10lf",x[i]);
    } 
    return 0; 
}

2.追っ手で解く

#include
#include
#include
#include
/**
*           
*/

#define N 5 

void display(double *temp) {
    int i  = 0;
    for(; i < N; i++, temp++) {
        printf("%-8lf
",*temp);
    } 
}

int main() {

    double a[N] = {0, 0, 6, 4, 1};
    double b[N] = {0, 8, 12, 9, 2}; 
    double c[N] = {0, 7, 5, 3, 0};

    double f[N] = {0, 0, -2, 8, 6}; //      

    double Y[N];
    double U[N];
    double X[N];        //     

    double L;
    int i,j;
    memset(Y,0,sizeof(Y));
    memset(U,0,sizeof(U));
    memset(X,0,sizeof(X));
    /**
    *1.  Y     
    *    (  ) 
    */ 
    Y[1] = f[1] / b[1];
    U[1] = c[1] / b[1];

    for (i = 2; i < N - 1; i++) {
        L = b[i] - a[i] * U[i - 1]; //  L   
        U[i] = c[i] / L;    
        Y[i] = (f[i] - a[i] * Y[i - 1]) / L; 
    }

    /**
    *    (  ) 
    */ 
    Y[N - 1] = (f[N-1] - Y[N - 2] * a[N - 1]) / (b[N -1] - U[N - 2] * a[N - 1]); 
    X[N-1] = Y[N - 1];
    for (i = N - 2; i > 0; i--) {       
        X[i] = Y[i] - U[i] * X[i + 1]; 
    }
    printf("
");
    display(X);

}

ヤクビーとガウスセデル反復法