math_util.c

multipleshooting to find the roots of the non-linear functions

     int dim1;
     {
        int *array;
  
        array = (int *)malloc(dim1 * sizeof(int));
  
        return( array );
  
     }
  /* ----------------------------------------------------------------------------
  
  A1D_FREE_INT
  
      This function frees the memory allocated by the function a1d_int.
  
  -----------------------------------------------------------------------------*/
  void a1d_free_int( array )
     int *array;
     {
  
        if( array == NULL )
           return;
        free( array );
     }
  /* ----------------------------------------------------------------------------
  
  A1D_ALLO_UINT
  
     This function allocates storage for a one dimensional unsigned int
     array of length dim1. The function returns a pointer to an array of ints.
  
  -----------------------------------------------------------------------------*/
  unsigned int *a1d_allo_uint( dim1 )
     int dim1;
     {
        unsigned int *array;
  
        array = (unsigned int *)malloc(dim1 * sizeof(unsigned int));
  
        return( array );
  
     }
  /* ----------------------------------------------------------------------------
  
  A1D_FREE_UINT
  
      This function frees the memory allocated by the function a1d_uint.
  
  -----------------------------------------------------------------------------*/
  void a1d_free_uint( array )
     unsigned int *array;
     {
  
        if( array == NULL )
           return;
        free( array );
     }
  /* ----------------------------------------------------------------------------
  
  A2D_ALLO_INT
  
      This function allocates storage for a two dimensional ints
      array. The array deminsions are (dim1,dim2). The function returns an
      array of pointers to ints.
  
  -----------------------------------------------------------------------------*/
  int **a2d_allo_int( dim1, dim2 )
     int dim1;
     int dim2;
     {
        int **array;
        int i;
        int j;
  
        array = (int **)malloc(dim1 * sizeof(int *));
  
        if( array == NULL )
           return( array );
  
        for( i = 0; i < dim1; i++ ) {
           array[i] = (int *)malloc(dim2 * sizeof(int));
           if( array[i] == NULL ) {
              for( j = i-1; j >= 0; j-- )
                 free( array[j] );
              free( array );
              array = NULL;
              return( array );
           }
        }
  
        return( array );
     }
  /* ----------------------------------------------------------------------------
  
  A2D_FREE_INT
  
      This function frees the memory allocated by the function a2d_allo.
  
  -----------------------------------------------------------------------------*/
  void a2d_free_int( array, dim1 )
     int **array;
     int dim1;
     {
        int i;
  
        if( array == NULL )
           return;
           
        for( i = dim1-1; i >= 0; i-- )
           if( array[i] != NULL )
              free( array[i] );
  
        if( array != NULL)
           free( array );
     }
  /* ----------------------------------------------------------------------------
  
  A2D_ALLO_UINT
  
      This function allocates storage for a two dimensional ints
      array. The array deminsions are (dim1,dim2). The function returns an
      array of pointers to ints.
  
  -----------------------------------------------------------------------------*/
  unsigned int **a2d_allo_uint( dim1, dim2 )
     int dim1;
     int dim2;
     {
        unsigned int **array;
        int i;
        int j;
  
        array = (unsigned int **)malloc(dim1 * sizeof(unsigned int *));
  
        if( array == NULL )
           return( array );
  
        for( i = 0; i < dim1; i++ ) {
           array[i] = (unsigned int *)malloc(dim2 * sizeof(unsigned int));
           if( array[i] == NULL ) {
              for( j = i-1; j >= 0; j-- )
                 free( array[j] );
              free( array );
              array = NULL;
              return( array );
           }
        }
  
        return( array );
     }
  /* ----------------------------------------------------------------------------
  
  A2D_FREE_UINT
  
      This function frees the memory allocated by the function a2d_allo.
  
  -----------------------------------------------------------------------------*/
  void a2d_free_uint( array, dim1 )
     unsigned int **array;
     int dim1;
     {
        int i;
  
        if( array == NULL )
           return;
           
        for( i = dim1-1; i >= 0; i-- )
           if( array[i] != NULL )
              free( array[i] );
  
        if( array != NULL)
           free( array );
     }
  /* ----------------------------------------------------------------------------
  
  VEC_NULL
  
      This function zeros a one dimensional array of long doubles. The length of the
      array is n.
                           {a} = {0}
  
  -----------------------------------------------------------------------------*/
  void vec_null( a, n )
     long double *a;
     int n;
     {
        int i;
  
        for( i = 0; i < n; i++ )
           *(a+i) = 0.0;
  
     }
  /* ----------------------------------------------------------------------------
  
  MAT_NULL
  
      This function zeros an m by n matrix of long doubles.
  
                       [a] = [0]
  
  -----------------------------------------------------------------------------*/
  void mat_null( a, m, n )
     long double **a;
     int m;
     int n;
     {
        int i;
        int j;
  
        for( i = 0; i < m; i++ )
           for( j = 0; j < n; j++ )
              a[i][j] = 0.0;
  
     }
  
  /* ----------------------------------------------------------------------------
  
  RK4
  
      This function uses the 4th order Runge-Kutta method to integrate a system
      of n ordinary differential equations. The integration goes from xstart to
      xend in nsteps. The result is put in y. On input, y = y(xstart).
  
      Return code -> 0 : Normal completion.
      Return code -> 1 : Error, number of equations <= 0.
      Return code -> 2 : Error, number of steps <= 0.
      Return code -> 3 : Error, xstart == xend.
      Return code -> 4 : Error, unable to allocate storage for coefficients.
  
  -----------------------------------------------------------------------------*/
  int rk4( fcn, n, y, xstart, xend, nsteps )
     void (*fcn)( int, long double, long double *, long double * );
     int n;
     long double *y;
     long double xstart;
     long double xend;
     int nsteps;
     {
        long double *a1d_allo_dbl( int );
        void a1d_free_dbl( long double * );
        long double *a;
        long double *b;
        long double *c;
        long double *d;
        long double *yprime;
        long double *yx;
        long double h;
        long double x;
        long double xstep;
        int i;
        int kstep;
  
  /*
     Do some error checks
  */
  
        if ( n <= 0 )
           return(1);
  
        if( nsteps <= 0 )
           return(2);
  
        if( xstart == xend )
           return(3);
  
  /*
     Compute step size
  */
  
        h = (xend - xstart) / (long double) nsteps;
  
  /*
     Allocate storage for a, b, c, d, yx and yprime
  */
  
        a = a1d_allo_dbl( n );
        b = a1d_allo_dbl( n );
        c = a1d_allo_dbl( n );
        d = a1d_allo_dbl( n );
        yx = a1d_allo_dbl( n );
        yprime = a1d_allo_dbl( n );
  
        if( a == NULL || b == NULL || c == NULL || d == NULL ||
            yx == NULL || yprime == NULL ) {
           a1d_free_dbl( a );
           a1d_free_dbl( b );
           a1d_free_dbl( c );
           a1d_free_dbl( d );
           a1d_free_dbl( yx );
           a1d_free_dbl( yprime );
           return(4);
        }
  
  /*
     Begin Runga-Kutta iterations
  */
  
        xstep = xstart;
  
        for( kstep = 1; kstep <= nsteps; kstep++ ) {
  /* printf(" STEP = %d\n",kstep); */
  /*
     Compute the a coefficients
  */
  
           x = xstep;
           for( i = 0; i < n; i++ )
              yx[i] = y[i];
  
           (*fcn)( n, x, yx, yprime );
  
           for( i = 0; i < n; i++ )
              a[i] = h * yprime[i];
  
  /*
     Compute the b coefficients
  */
  
           x = xstep + 0.5 * h;
           for( i = 0; i < n; i++ )
              yx[i] = y[i] + 0.5 * a[i];
  
           (*fcn)( n, x, yx, yprime );
  
           for( i = 0; i < n; i++ )
              b[i] = h * yprime[i];
  
  /*
     Compute the c coefficients
  */
  
           x = xstep + 0.5 * h;
           for( i = 0; i < n; i++ )
              yx[i] = y[i] + 0.5 * b[i];
  
           (*fcn)( n, x, yx, yprime );
  
           for( i = 0; i < n; i++ )
              c[i] = h * yprime[i];
  
  /*
     Compute the d coefficients
  */
  
           x = xstep + h;
           for( i = 0; i < n; i++ )
              yx[i] = y[i] + c[i];
  
           (*fcn)( n, x, yx, yprime );
  
           for( i = 0; i < n; i++ )
              d[i] = h * yprime[i];
  
  /*
     Compute the new value of y
  */
  
           for( i = 0; i < n; i++ )
              y[i] = y[i] + (a[i] + 2.0 * b[i] + 2.0 * c[i] + d[i]) / 6.0;           
  
  /*
     Increment xstep, and continue loop
  */
  
           xstep += h;
        }
  
  /*
     Free all variable allocated
  */
  
        a1d_free_dbl( a );
        a1d_free_dbl( b );
        a1d_free_dbl( c );
        a1d_free_dbl( d );
        a1d_free_dbl( yx );
        a1d_free_dbl( yprime );
        return(0);
     }
  /* ----------------------------------------------------------------------------
  
  VEC_COPY
  
      This function copies array a into array b. Both arrays are of demension n.
  
                               {a} = {b}
  
  -----------------------------------------------------------------------------*/
  void vec_copy( a, b, n )
     long double *a;
     long double *b;
     int n;
     {
        int i;
  
        for( i = 0; i < n; i++ )
           *(a+i) = *(b+i);
  
    }
  /* ----------------------------------------------------------------------------
  
  VEC_SUB
  
      This function substracts array b from c and puts the result in a. Both
      arrays are long double of length n.
  
                               {a} = {b} - {c}
  
  -----------------------------------------------------------------------------*/
  void vec_sub( a, b, c, n )
     long double *a;
     long double *b;
     long double *c;
     int n;
     {
        int i;
  
        for( i = 0; i < n; i++ )
           *(a+i) = *(b+i) - *(c+i);
  
     }
  /* ----------------------------------------------------------------------------
  
  VEC_MUL_BCT
  
      This function multiplies an n by 1 vector by a 1 by n vector to produce a
      n by n matrix. All arrays are long double.
      
                                           T    
                               [a] = {b}{c}
  
  -----------------------------------------------------------------------------*/
  void vec_mul_bct( a, b, c, n )
     long double **a;
     long double *b;
     long double *c;
     int n;
     {
        int i, j;
  
        for( i = 0; i < n; i++ )
           for( j = 0; j < n; j++ )
              a[i][j] = b[i]*c[j];
  
     }
  /* ----------------------------------------------------------------------------
  
  VEC_MUL_BTC
  
      This function multiplies an 1 by n vector by a n by 1 vector to produce a
      scalar. Both arrays are long double.
      
                                        T       
                             return( {b}{c} )
  
  -----------------------------------------------------------------------------*/