高斯消元法

高斯消元法

 

方法1：

 

#include "stdio.h"

#include "stdlib.h"

 

void RKT ( t , y , n , h , k , z )

int n ;                                      /* 微分方程组中方程的个数，也是未知函数的个数 */

int k ;                                      /* 积分的步数 ( 包括起始点这一步 )*/

double t ;                     /* 积分的起始点 t0*/

double h ;                    /* 积分的步长 */

double y [];                           /* 存放 n 个未知函数在起始点 t 处的函数值 , 返回时 , 其初值在二维数组 z 的第零列中 */

double z [];                           /* 二维数组 , 体积为 n x k. 返回 k 个积分点上的 n 个未知函数值 */

{

         extern void Func ();                                  /* 声明要求解的微分方程组 */

    int i , j , l ;

    double a [4],* b ,* d ;

    b = malloc ( n * sizeof ( double ));                  /* 分配存储空间 */

         if ( b == NULL )

         {

                   printf ( " 内存分配失败 /n" );

                   exit (1);

         }

    d = malloc ( n * sizeof ( double ));                  /* 分配存储空间 */

         if ( d == NULL )

         {

                   printf ( " 内存分配失败 /n" );

                   exit (1);

         }

         /* 后面应用 RK4 公式中用到的系数 */

    a [0]= h /2.0;                                               

         a [1]= h /2.0;

    a [2]= h ;

         a [3]= h ;

    for ( i =0; i <= n -1; i ++)

                  z [ i * k ]= y [ i ];                                     /* 将初值赋给数组 z 的相应位置 */

    for ( l =1; l <= k -1; l ++)

    {

                  Func ( y , d );

        for ( i =0; i <= n -1; i ++)

                            b [ i ]= y [ i ];

        for ( j =0; j <=2; j ++)

        {

                            for ( i =0; i <= n -1; i ++)

            {

                                     y [ i ]= z [ i * k + l -1]+ a [ j ]* d [ i ];

                b [ i ]= b [ i ]+ a [ j +1]* d [ i ]/3.0;

            }

            Func ( y , d );

        }

        for ( i =0; i <= n -1; i ++)

          y [ i ]= b [ i ]+ h * d [ i ]/6.0;

        for ( i =0; i <= n -1; i ++)

          z [ i * k + l ]= y [ i ];

        t = t + h ;

         }

    free ( b );                          /* 释放存储空间 */

         free ( d );                         /* 释放存储空间 */

    return ;

}

main ()

{

         int i , j ;

    double t , h , y [3], z [3][11];

    y [0]=-1.0;

         y [1]=0.0;

         y [2]=1.0;

    t =0.0;

         h =0.01;

    RKT ( t , y ,3, h ,11, z );

    printf ( "/n" );

    for ( i =0; i <=10; i ++)                            /* 打印输出结果 */

    {

                   t = i * h ;

        printf ( "t=%5.2f/t   " , t );

        for ( j =0; j <=2; j ++)

          printf ( "y(%d)=%e " , j , z [ j ][ i ]);

        printf ( "/n" );

         }

}

 

void Func ( y , d )

double y [], d [];

{

         d [0]= y [1];          /*y0'=y1*/

         d [1]=- y [0];                  /*y1'=y0*/

         d [2]=- y [2];                  /*y2'=y2*/

         return ;

}

高斯消元法：

#include < stdio .h>

#include <stdlib.h>

#include < malloc .h>

#include <math.h>

 

int GS ( int , double **, double *, double );

double ** TwoArrayAlloc ( int , int );

void TwoArrayFree ( double **);

 

void main ()

{

         int i , n ;

         double ep ,** a ,* b ;

         n = 3;

         ep = 1e-4;

         a = TwoArrayAlloc ( n , n );

         b = ( double *) calloc ( n , sizeof ( double ));

         if ( b == NULL )

         {

                   printf ( " 内存分配失败 /n" );

                   exit (1);

         }

         a [0][0]= 2; a [0][1]= 6; a [0][2]=-1;

         a [1][0]= 5; a [1][1]=-1; a [1][2]= 2;

         a [2][0]=-3; a [2][1]=-4; a [2][2]= 1;

         b [0] = -12; b [1] = 29;  b [2] = 5;

         if (! GS ( n , a , b , ep ))

         {

                   printf ( " 不可以用高斯消去法求解 /n" );

                   exit (0);

         }

         printf ( " 该方程组的解为： /n" );

         for ( i =0; i <3; i ++)

                  printf ( "x%d = %.2f/n" , i , b [ i ]);

         TwoArrayFree ( a );

         free ( b );

}

 

int GS ( n , a , b , ep )

int n ;

double ** a ;

double * b ;

double ep ;

{

         int i , j , k , l ;

         double t ;

         for ( k =1; k <= n ; k ++)

         {

                  for ( l = k ; l <= n ; l ++)

                            if ( fabs ( a [ l -1][ k -1])> ep )

                                     break ;

                            else if ( l == n )

                                     return (0);

                   if ( l != k )

                   {

                            for ( j = k ; j <= n ; j ++)

                            {

                                     t = a [ k -1][ j -1];

                                     a [ k -1][ j -1]= a [ l -1][ j -1];

                                     a [ l -1][ j -1]= t ;

                            }

                            t = b [ k -1];

                            b [ k -1]= b [ l -1];

                            b [ l -1]= t ;

                   }

                  t =1/ a [ k -1][ k -1];

                  for ( j = k +1; j <= n ; j ++)

                            a [ k -1][ j -1]= t * a [ k -1][ j -1];

                  b [ k -1]*= t ;

                  for ( i = k +1; i <= n ; i ++)

                   {

                            for ( j = k +1; j <= n ; j ++)

                                     a [ i -1][ j -1]-= a [ i -1][ k -1]* a [ k -1][ j -1];

                            b [ i -1]-= a [ i -1][ k -1]* b [ k -1];

                   }

         }

         for ( i = n -1; i >=1; i --)

                  for ( j = i +1; j <= n ; j ++)

                            b [ i -1]-= a [ i -1][ j -1]* b [ j -1];

return (1);

}

 

double ** TwoArrayAlloc ( int r , int c )

{

         double * x ,** y ;

         int n ;

         x =( double *) calloc ( r * c , sizeof ( double ));

         y =( double **) calloc ( r , sizeof ( double *));

         for ( n =0; n <= r -1;++ n )

                  y [ n ]=& x [ c * n ];

         return ( y );

}

 

void TwoArrayFree ( double ** x )

{

         free ( x [0]);

         free ( x );

}