コードは次のようになる。
#include <stdio.h> #include <stdlib.h> #include <memory.h> /************************************************************ * ガウスの消去法による線形方程式の数値計算 * * (C) 2012- Gofer. ************************************************************/ // 方程式の表示 void display_equtation(int N, double *A) { int i, j; for(i=0; i<N; ++i) { for(j=0; j<N; ++j) { printf("%f x[%d]", A[i*(N+1) + j], j); if(j < N-1) printf(" + "); } printf(" = %f\n", A[i*(N+1) + N]); } putchar('\n'); return; } // n行目とm行目を入れ替える double* replace_row(int n, int m, int N, double* A) { double* line_temp = (double*)malloc(sizeof(double)*(N+1)); memmove(line_temp, &A[n*(N+1)], sizeof(double)*(N+1)); memmove(&A[n*(N+1)], &A[m*(N+1)], sizeof(double)*(N+1)); memmove(&A[m*(N+1)], line_temp, sizeof(double)*(N+1)); free(line_temp); return A; } // i行目の前進消去を行う double* forward_elimination(int i, int N, double *A) { int j, k, J = -1; double divide; // 対角要素が0ならピボット選択を行う if(A[i*(N+1) + i] == 0) { for(j=i+1; j<N; ++j) { if(A[j*(N+1) + i] != 0) { J = j; break; } } if(J == -1) return NULL; replace_row(i, J, N, A); } // 割算する係数を決定 divide = A[i*(N+1) + i]; // 標準化 for(j=0; j<(N+1); ++j) A[i*(N+1) + j] /= divide; // i行目以外の行から引く for(k=i+1; k<N; ++k) { for(j=(N+1)-1; j>=0; --j) A[k*(N+1) + j] -= A[k*(N+1) + i] * A[i*(N+1) + j]; } return A; } // i行目の後退代入 double* backward_substitution(int i, int N, double *A, double *x) { // (N-1)-i個の解は決定済み int j, k; // 1ステップ前に求めた解を代入 k = i+1; if(k < N) { for(j=k-1; j>=0; --j) { A[j*(N+1) + N] -= x[k] * A[j*(N+1) + k]; A[j*(N+1) + k] = 0; } } // 第i行を用いて1文字消去 if(A[i*(N+1) + i] != 0) { A[i*(N+1) + N] /= A[i*(N+1) + i]; A[i*(N+1) + i] = 1; } // xベクタに代入 x[i] = A[i*(N+1) + N]; return A; } int main(void) { //int N = 2; /* double A[2][3] = { {1, 2, 5}, {2, 3, 8}}; */ //int N = 3; /* double A[3][4] = { {2, 4, 2, 8}, {4, 10, 3, 17}, {3, 7, 1, 11}}; */ //int N = 4; /* double A[4][5] = { {1, 1, 2, 1, 1}, {1, 1, 3, 2, -2}, {2, -2, 2, -1, 1}, {-1, 1, 0, 1, -1}}; */ int N = 4; double A[4][5] = { {1, -1, 2, -1, -8}, {2, -2, 3, -3, -20}, {1, 1, 1, 0, -2}, {1, -1, 4, 3, 4}}; int i; double* x = (double*)malloc(sizeof(double) * N); display_equtation(N, (double*)A); // 前進消去 for(i=0; i<N; ++i) { if(forward_elimination(i, N, (double*)A) == NULL) { fprintf(stderr, "Error!\n"); free(x); return -1; } } // 後退代入 for(i=N-1; i>=0; --i) backward_substitution(i, N, (double*)A, x); for(i=0; i<N; ++i) printf("x[%d] = %f\n", i, x[i]); free(x); return 0; }
0 件のコメント:
コメントを投稿