/* input matrix
-x1 + x2 + 2x3 = 2
3x1 - x2 + x3 = 6
-x1 + 3x2 + 4x3 = 4
*/
/* ex6-1.c based on Gaussian Elimination method
* for solving the n x n linear algebra system
* a11 x1+a12 x2+...+a1n xn=b1
* a21 x1+a22 x2+...+a2n xn=b2
* . . . .
* . . . .
* an1 x1+an2 x2+...+ann xn=bn
* Input number of unknowns and equations n
* with coefficent a11,a12,...,ann and b1,b2,
* ...bn. Output solution x1,x2,x3,...,xn.
*/
#include <stdio.h>
#include <math.h>
#define MAX 3
void gaussh(int n,double a[MAX][MAX+1],double x[]);
void printmatrix(int n , double a[MAX][MAX+1]) ;
void backSub(double a[MAX][MAX+1]) ;
int main()
{
int i,j,k,m,n;
n=MAX;
double x[MAX];
double a[MAX][MAX+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");
printf("================================\n");
printmatrix(n,a);
gaussh(n,a,x); /* call the function gaussh() */
return 0;
}
void gaussh(int n,double a[MAX][MAX+1],double x[])
{
int i,j,k,m;
double temp,bb,cc;
for(k=0;k<=n-1;k++)
{
/* check if a[k][k]=0 is true then interchange */
/* E(k) and E(k+1).............................*/
if(a[k][k]==0)
{
for(m=0;m<=n;m++)
{
temp=a[k][m];
a[k][m]=a[k+1][m];
a[k+1][m]=temp;
}
}
printf("\n");
printf("================================\n");
printmatrix(n,a);
/* To reduce the matrix to triangular form */
for(i=k;i<=n-1;i++)
{
bb=a[i+1][k]/a[k][k];
for(j=k;j<=n+1;j++)
a[i+1][j]=a[i+1][j]-bb*a[k][j];
}
}
if(fabs(a[n-1][n-1])==0.0)
{
printf("NO UNIQUE SOLUTION!!!\n");
return;
}
printf("\n\n上三角矩陣\n");
printf("================================\n");
printmatrix(n,a);
backSub(a);
/* To start backward substitution */
/*
x[n]=a[n-1][n]/a[n-1][n-1];
for(i=n;i>=0;i--)
{
cc=0.0;
for(j=i+1;j<=n;j++)
cc=cc+a[i][j]*x[j];
x[i]=(a[i][n+1]-cc)/a[i][i];
}
*/
return;
}
void printmatrix(int n , double a[MAX][MAX+1])
{
for (int i=0; i<n; i++)
{
for (int j=0; j<=n; j++)
printf("%4.2lf ", a[i][j]);
printf("\n");
}
}
// function to calculate the values of the unknowns
void backSub(double a[MAX][MAX+1])
{
double x[MAX]; // An array to store solution
/* Start calculating from last equation up to the
first */
for (int i = MAX-1; i >= 0; i--)
{
/* start with the RHS of the equation */
x[i] = a[i][MAX];
/* Initialize j to i+1 since matrix is upper
triangular*/
for (int j=i+1; j<MAX; j++)
{
/* subtract all the lhs values
* except the coefficient of the variable
* whose value is being calculated */
x[i] -= a[i][j]*x[j];
}
/* divide the RHS by the coefficient of the
unknown being calculated */
x[i] = x[i]/a[i][i];
}
printf("\n行列式 解 \n");
printf("================================\n");
for (int i=0; i<MAX; i++)
printf("x%1d = %3.3lf\n", i, x[i]);
}
輸出畫面
原始行列式
================================
-1.00 1.00 2.00 2.00
3.00 -1.00 1.00 6.00
-1.00 3.00 4.00 4.00
================================
-1.00 1.00 2.00 2.00
3.00 -1.00 1.00 6.00
-1.00 3.00 4.00 4.00
================================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 2.00 2.00 2.00
================================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 0.00 -5.00 -10.00
上三角矩陣
================================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 0.00 -5.00 -10.00
行列式 解
================================
x0 = 1.000
x1 = -1.000
x2 = 2.000
沒有留言:
張貼留言