2019年5月3日 星期五

C語言 例題6-1 (a) 利用Guassian Elimination 解行列式

C語言 例題6-1 (a) 利用Guassian Elimination 解行列式

/* 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.

 3
-1.0   1.0  2.0  2.0
 3.0  -1.0  1.0  6.0
-1.0   3.0  4.0  4.0


 3
 2.0  -1.5  3.0  1.0
-1.0   0.0  2.0  3.0
 4.0  -4.5  5.0  1.0


 4
2.0   0.0  0.0  0.0   3.0
1.0   1.5  0.0  0.0   4.5
0.0  -3.0  0.5  0.0  -6.6
2.0  -2.0  1.0  1.0   0.8


 4
 4.01  1.23   1.43  -0.73   5.94
 1.23  7.41   2.41   3.02  14.07
 1.43  2.41   5.79  -1.11   8.52
-0.73  3.02  -1.11   6.41   7.59


  4
 1.0   1.0   0.0  1.0   2.0
 2.0   1.0   -1.0  1.0   1.0
 4.0   -1.0   -2.0  2.0  0.0
 3.0  -1.0  -1.0  2.0   -3.0



 */
#include <stdio.h>
#include <math.h>
#define MAX 20
void gaussh( int m, int n ,double arr[][n], double dist[]);
void print(double *arr, int m, int n);
void Copy2Darray(double *arr, int m, int n , double *dist ) ;

int N[]={3,3,4,4,4};

double a1[][4]= { {-1.0 ,   1.0  , 2.0 ,  2.0},
                    { 3.0 ,  -1.0  , 1.0 ,  6.0},
                    {-1.0 ,   3.0  , 4.0 ,  4.0}
                };

double a2[][4]= { {  2.0 , -1.5 , 3.0 , 1.0 },
                  { -1.0 ,  0.0 , 2.0 , 3.00},
                  {  4.0 , -4.5 , 5.0 , 1.00}
                };



int  main()
{
    int i,j,k,m1,n1,n;
    double a[MAX][MAX],x[MAX];
    n=3;
    m1 = 3;
    n1 = 4;
    Copy2Darray((double *)a1, m1, n1, (double *)a ) ;   //copy a1 --> a array

    m1 = 3;
    n1 = 4;
 
    // We can also use "print(&arr[0][0], m, n);"
    print((double *)a, m1, n1);
   
    gaussh(m1, n1, a,x ); /* call the function gaussh() */
   
    for(i=0;i<m1;i++)
       printf("x%d=%6.3lf\n",i,x[i]);
   
    return 0;
}

void print(double *arr, int m, int n)
{
    int i, j;
   
    for (i = 0; i < m; i++)
    {
      for (j = 0; j < n; j++)
        printf("%4.2lf ", *((arr+i*n) + j));
      printf("\n");   
    }
}

void Copy2Darray(double *arr, int m, int n , double *dist )
{
    int i, j;
   
    for (i = 0; i < m; i++)
    {
      for (j = 0; j < n; j++)
        {
            *((dist+i*n) + j)=*((arr+i*n) + j);
            //printf("%4.2lf ", *((arr+i*n) + j));
        }
      //printf("\n");   
    }
}

void gaussh( int m, int n , double arr[][n], double dist[])
{
    int i,j,k,l,cn=0;
    double temp,bb,cc;
    for(k=0;k<m;k++)
    {
       /* check if a[k][k]=0 is true then interchange */
       /* E(k) and E(k+1).............................*/
       if(arr[k][k]==0)
       {
        for(l=1;l<=m+1;l++)
        {
            temp=arr[k][l];
            arr[k][l]=arr[k+1][l];
            arr[k+1][l]=temp;
        }
       }
       /* To reduce the matrix to triangular form */
       for(i=k;i<=m-1;i++)
       {
        bb=arr[i+1][k]/arr[k][k];
        for(j=k;j<=m+1;j++)
            arr[i+1][j]=arr[i+1][j]-bb*arr[k][j];
       }
     
       printf("\n");
       printf("===========%1d===================\n",cn);     
       print((double *)arr, m, n); 
       cn++;
   
    }
   
    printf("\n");
    printf("===========上三角矩陣================\n");
    print((double *)arr, m, n); 

    if(fabs(arr[n][n])==0.0)
    {
        printf("NO UNIQUE SOLUTION!!!\n");
        return ;
    }
   
    printf("\n行列式 解 \n");
    printf("================================\n");

/* To start backward substitution */
/* To start backward substitution */
/* Start calculating from last equation up to the
       first */
    for (int i = m-1; i >= 0; i--)
    {
        /* start with the RHS of the equation */
        dist[i] = arr[i][m];
          /* Initialize j to i+1 since matrix is upper
           triangular*/
        for (int j=i+1; j<m; j++)
        {
            /* subtract all the lhs values
             * except the coefficient of the variable
             * whose value is being calculated */
            dist[i] -= arr[i][j]*dist[j];
        }
 
        /* divide the RHS by the coefficient of the
           unknown being calculated */
        dist[i] = dist[i]/arr[i][i];
    }


    return;
}


輸出畫面
-1.00 1.00 2.00 2.00
3.00 -1.00 1.00 6.00
-1.00 3.00 4.00 4.00

===========0===================
-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===================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 0.00 -5.00 -10.00

===========2===================
-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

沒有留言:

張貼留言

Messaging API作為替代方案

  LINE超好用功能要沒了!LINE Notify明年3月底終止服務,有什麼替代方案? LINE Notify將於2025年3月31日結束服務,官方建議改用Messaging API作為替代方案。 //CHANNEL_ACCESS_TOKEN = 'Messaging ...