-x1 + x2 + 2x3 = 2
3x1 - x2 + x3 = 6
-x1 + 3x2 + 4x3 = 4
參考 https://www.geeksforgeeks.org/gaussian-elimination/
// C program to demostrate working of Guassian Elimination
// method
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 3 // Number of unknowns
// function to reduce matrix to r.e.f. Returns a value to
// indicate whether matrix is singular or not
int forwardElim(double mat[N][N+1]);
// function to calculate the values of the unknowns
void backSub(double mat[N][N+1]);
// function to get matrix content
void gaussianElimination(double mat[N][N+1])
{
/* reduction into r.e.f. */
int singular_flag = forwardElim(mat);
/* if matrix is singular */
if (singular_flag != -1)
{
printf("Singular Matrix.\n");
/* if the RHS of equation corresponding to
zero row is 0, * system has infinitely
many solutions, else inconsistent*/
if (mat[singular_flag][N])
printf("Inconsistent System.");
else
printf("May have infinitely many "
"solutions.");
return;
}
/* get solution to system and print it using
backward substitution */
backSub(mat);
}
// function for elemntary operation of swapping two rows
void swap_row(double mat[N][N+1], int i, int j)
{
//printf("Swapped rows %d and %d\n", i, j);
for (int k=0; k<=N; k++)
{
double temp = mat[i][k];
mat[i][k] = mat[j][k];
mat[j][k] = temp;
}
}
// function to print matrix content at any stage
void print(double mat[N][N+1])
{
for (int i=0; i<N; i++, printf("\n"))
for (int j=0; j<=N; j++)
printf("%lf ", mat[i][j]);
printf("\n");
}
// function to reduce matrix to r.e.f.
int forwardElim(double mat[N][N+1])
{
for (int k=0; k<N; k++)
{
// Initialize maximum value and index for pivot
int i_max = k;
int v_max = mat[i_max][k];
/* find greater amplitude for pivot if any */
for (int i = k+1; i < N; i++)
if (abs(mat[i][k]) > v_max)
v_max = mat[i][k], i_max = i;
/* if a prinicipal diagonal element is zero,
* it denotes that matrix is singular, and
* will lead to a division-by-zero later. */
if (!mat[k][i_max])
return k; // Matrix is singular
/* Swap the greatest value row with current row */
if (i_max != k)
swap_row(mat, k, i_max);
printf("第 %1d 次的行列式 \n",k);
print(mat);
for (int i=k+1; i<N; i++)
{
/* factor f to set current row kth elemnt to 0,
* and subsequently remaining kth column to 0 */
double f = mat[i][k]/mat[k][k];
/* subtract fth multiple of corresponding kth
row element*/
for (int j=k+1; j<=N; j++)
mat[i][j] -= mat[k][j]*f;
/* filling lower triangular matrix with zeros*/
mat[i][k] = 0;
}
print(mat); //for matrix state
}
printf("上三角矩陣\n");
print(mat); //for matrix state
return -1;
}
// function to calculate the values of the unknowns
void backSub(double mat[N][N+1])
{
double x[N]; // An array to store solution
/* Start calculating from last equation up to the
first */
for (int i = N-1; i >= 0; i--)
{
/* start with the RHS of the equation */
x[i] = mat[i][N];
/* Initialize j to i+1 since matrix is upper
triangular*/
for (int j=i+1; j<N; j++)
{
/* subtract all the lhs values
* except the coefficient of the variable
* whose value is being calculated */
x[i] -= mat[i][j]*x[j];
}
/* divide the RHS by the coefficient of the
unknown being calculated */
x[i] = x[i]/mat[i][i];
}
printf("\nSolution for the system:\n");
for (int i=0; i<N; i++)
printf("%lf\n", x[i]);
}
// Driver program
int main()
{
/* input matrix
-x1 + x2 + 2x3 = 2
3x1 - x2 + x3 = 6
-x1 + 3x2 + 4x3 = 4
*/
double mat[N][N+1] = {{-1.0, 1.0,2.0, 2.0},
{3.0, -1.0, 1.0, 6.0},
{-1.0, 3.0 , 4.0, 4.0}
};
printf("原始行列式\n");
print(mat); //for matrix state
gaussianElimination(mat);
return 0;
}
輸出畫面
原始行列式
-1.000000 1.000000 2.000000 2.000000
3.000000 -1.000000 1.000000 6.000000
-1.000000 3.000000 4.000000 4.000000
第 0 次的行列式
3.000000 -1.000000 1.000000 6.000000
-1.000000 1.000000 2.000000 2.000000
-1.000000 3.000000 4.000000 4.000000
3.000000 -1.000000 1.000000 6.000000
0.000000 0.666667 2.333333 4.000000
0.000000 2.666667 4.333333 6.000000
第 1 次的行列式
3.000000 -1.000000 1.000000 6.000000
0.000000 2.666667 4.333333 6.000000
0.000000 0.666667 2.333333 4.000000
3.000000 -1.000000 1.000000 6.000000
0.000000 2.666667 4.333333 6.000000
0.000000 0.000000 1.250000 2.500000
第 2 次的行列式
3.000000 -1.000000 1.000000 6.000000
0.000000 2.666667 4.333333 6.000000
0.000000 0.000000 1.250000 2.500000
3.000000 -1.000000 1.000000 6.000000
0.000000 2.666667 4.333333 6.000000
0.000000 0.000000 1.250000 2.500000
上三角矩陣
3.000000 -1.000000 1.000000 6.000000
0.000000 2.666667 4.333333 6.000000
0.000000 0.000000 1.250000 2.500000
Solution for the system:
1.000000
-1.000000
2.000000
沒有留言:
張貼留言