連立方程式

#include <stdio.h>
#include <math.h>

#define DIM  4

int pivot( double a[ DIM ][ DIM + 1 ], int i ){

      double p = fabs( a[ i ][ i ] ), t;
      int ii, j, k = i;

      for ( ii = i + 1; ii < DIM; ii++ ){
          if ( fabs( a[ ii ][ i ] ) > p ){
              k = ii;
              p = fabs( a[ ii ][ i ] );
          }
      }
      if ( k != i ){
         for ( j = i; j < DIM + 1; j++ ){
             t = a[ k ][ j ];
             a[ k ][ j ] = a[ i ][ j ];
             a[ i ][ j ] = t;
         }
      }
      return k;
}

int main(){

    double a[ DIM ][ DIM + 1 ] = { { 2, -3, 1, -1, -3 }, 
                                   { 1, 2, -3, 1, 6 }, 
                                   { 3, -1, 1, 2, -4 },
                                   { 2, 1,  2, 0, 2 } };
    double x[ DIM ], s, t;
    int i, ii, j, imax, jmax;

    imax = DIM;
    jmax = DIM + 1;

/*  forward elimination */
    for ( i = 0; i < imax; i++ ){
       pivot( a, i );
       for ( ii = i + 1; ii < imax; ii++ ){
           t = a[ ii ][ i ] / a[ i ][ i ];
           for ( j = 0; j < jmax; j++ ) 
               a[ ii ][ j ] -= t * a[ i ][ j ];			
       }
    }

/*  backward substitution */
    for ( i = imax - 1; i >= 0; i-- ){
       for ( j = i + 1, s = 0; j < imax; j++ ) 
           s += a[ i ][ j ] * x[ j ];
       x[ i ] = ( a[ i ][ jmax - 1 ] - s ) / a[ i ][ i ];
    }

/*  output solutions */
    for ( i = 0; i < imax; i++ )
       printf( "x( %d ) = %10.5f\n", i, x[ i ] );

	return 0;

}

戻る

SEO		[PR] 爆速!無料ブログ 無料ホームページ開設 無料ライブ放送