高斯消去法過程很直覺,先將左下三角消成 0,在將中間對角線消成 1,最後將右上三角消成 0,此時方程組之解就在矩陣最右邊,詳細說明可以參考 高斯消去法 | 線代啟示錄。
目前沒有加入 nonsingular matrix 的判斷,當方程組有多個解或無解時,此程式會有問題,所以正確的方法應該要在運算前判斷矩陣是否非奇異。
/* 高斯消去法 Gaussian Elimination */
int8_t GaussianElimination( uint8_t col, uint8_t row, float *arr )
{
float tmp;
int8_t i = 0, j = 0, k = 0;
/* check nonsingular */
// if matrix is nonsingular
// return ERROR;
/* left-down to zero */
for(i = 0; i < col; i++) {
for(j = i + 1; j < row; j++) {
for(k = col; k > i - 1; k--) {
arr[j*row+k] = arr[j*row+k] - arr[j*row+i] / arr[i*row+i] * arr[i*row+k];
}
}
}
/* diagonal to one */
for(i = 0; i < col; i++) {
tmp = arr[i*row+i];
for(j = i; j < row; j++) {
arr[i*row+j] = arr[i*row+j] / tmp;
}
}
/* right-up to zero */
for(j = row - 2; j > 0; j--) {
for(i = 0; i < j; i++) {
arr[i*row+row-1] = arr[i*row+row-1] - arr[i*row+j] * arr[j*row+row-1];
}
}
return SUCCESS;
}
/* 主程式 */
#define col_n 4
#define row_n 5
float arr[col_n][row_n] = {
{1.0f,3.0f,-5.0f,2.0f, 94.0f}, // 1*x1 + 3*x2 - 5*x3 + 2*x4 = 94
{4.0f,-3.0f,1.0f,5.0f, 58.0f}, // 4*x1 - 3*x2 + 1*x3 + 5*x4 = 58
{6.0f,-2.0f,2.0f,4.0f, 72.0f}, // 6*x1 - 2*x2 + 2*x3 + 4*x4 = 72
{0.0f,2.0f,3.0f,-7.0f, -69.0f} // 0*x1 + 2*x2 + 3*x3 - 7*x4 = -69
};
int main()
{
uint8_t i = 0;
float *ptr = arr;
GaussianElimination(col_n, row_n, arr);
for(i = 0; i < col_n; i++) {
printf("%.1f\n", ptr[i*row_n+(row_n-1)]);
}
}
沒有留言:
張貼留言