mathieu_charv.c

来自「math library from gnu」· C语言 代码 · 共 850 行 · 第 1/2 页

          a1 = aa;

          if (fa == fa1)
          {
              result->err = GSL_DBL_EPSILON;
              break;
          }
          aa -= (aa - a2)/(fa - fa1)*fa;
          dela = fabs(aa - a2);
          if (dela < GSL_DBL_EPSILON)
          {
              result->err = GSL_DBL_EPSILON;
              break;
          }
          if (ii > 20)
          {
              result->err = dela;
              break;
          }
          fa1 = fa;
          ii++;
      }

      /* If the solution found is not near the original approximation,
         tweak the approximate value, and try again. */
      if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)))
      {
          counter++;
          if (counter == maxcount)
          {
              result->err = fabs(aa - aa_orig);
              break;
          }
          if (aa > aa_orig)
              aa = aa_orig - da*counter;
          else
              aa = aa_orig + da*counter;

          continue;
      }
      else
          break;
  }

  result->val = aa;
      
  /* If we went through the maximum number of retries and still didn't
     find the solution, let us know. */
  if (counter == maxcount)
  {
      GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);
  }
  
  return GSL_SUCCESS;
}


int gsl_sf_mathieu_b(int order, double qq, gsl_sf_result *result)
{
  int even_odd, nterms = 50, ii, counter = 0, maxcount = 200;
  double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;


  even_odd = 0;
  if (order % 2 != 0)
      even_odd = 1;

  /* The order cannot be 0. */
  if (order == 0)
  {
      GSL_ERROR("Characteristic value undefined for order 0", GSL_EFAILED);
  }

  /* If the argument is 0, then the coefficient is simply the square of
     the order. */
  if (qq == 0)
  {
      result->val = order*order;
      result->err = 0.0;
      return GSL_SUCCESS;
  }

  /* Use symmetry characteristics of the functions to handle cases with
     negative order and/or argument q.  See Abramowitz & Stegun, 20.8.3. */
  if (order < 0)
      order *= -1;
  if (qq < 0.0)
  {
      if (even_odd == 0)
          return gsl_sf_mathieu_b(order, -qq, result);
      else
          return gsl_sf_mathieu_a(order, -qq, result);
  }
  
  /* Compute an initial approximation for the characteristic value. */
  aa = approx_s(order, qq);
  
  /* Save the original approximation for later comparison. */
  aa_orig = aa;
  
  /* Loop as long as the final value is not near the approximate value
     (with a max limit to avoid potential infinite loop). */
  while (counter < maxcount)
  {
      a1 = aa + 0.001;
      ii = 0;
      if (even_odd == 0)
          fa1 = seer(order, qq, a1, nterms);
      else
          fa1 = seor(order, qq, a1, nterms);

      for (;;)
      {
          if (even_odd == 0)
              fa = seer(order, qq, aa, nterms);
          else
              fa = seor(order, qq, aa, nterms);
      
          a2 = a1;
          a1 = aa;

          if (fa == fa1)
          {
              result->err = GSL_DBL_EPSILON;
              break;
          }
          aa -= (aa - a2)/(fa - fa1)*fa;
          dela = fabs(aa - a2);
          if (dela < 1e-18)
          {
              result->err = GSL_DBL_EPSILON;
              break;
          }
          if (ii > 20)
          {
              result->err = dela;
              break;
          }
          fa1 = fa;
          ii++;
      }
      
      /* If the solution found is not near the original approximation,
         tweak the approximate value, and try again. */
      if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)))
      {
          counter++;
          if (counter == maxcount)
          {
              result->err = fabs(aa - aa_orig);
              break;
          }
          if (aa > aa_orig)
              aa = aa_orig - da*counter;
          else
              aa = aa_orig + da*counter;
          
          continue;
      }
      else
          break;
  }
  
  result->val = aa;
      
  /* If we went through the maximum number of retries and still didn't
     find the solution, let us know. */
  if (counter == maxcount)
  {
      GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);
  }
  
  return GSL_SUCCESS;
}


/* Eigenvalue solutions for characteristic values below. */


/*  figi.c converted from EISPACK Fortran FIGI.F.
 *
 *   given a nonsymmetric tridiagonal matrix such that the products
 *    of corresponding pairs of off-diagonal elements are all
 *    non-negative, this subroutine reduces it to a symmetric
 *    tridiagonal matrix with the same eigenvalues.  if, further,
 *    a zero product only occurs when both factors are zero,
 *    the reduced matrix is similar to the original matrix.
 *
 *    on input
 *
 *       n is the order of the matrix.
 *
 *       t contains the input matrix.  its subdiagonal is
 *         stored in the last n-1 positions of the first column,
 *         its diagonal in the n positions of the second column,
 *         and its superdiagonal in the first n-1 positions of
 *         the third column.  t(1,1) and t(n,3) are arbitrary.
 *
 *    on output
 *
 *       t is unaltered.
 *
 *       d contains the diagonal elements of the symmetric matrix.
 *
 *       e contains the subdiagonal elements of the symmetric
 *         matrix in its last n-1 positions.  e(1) is not set.
 *
 *       e2 contains the squares of the corresponding elements of e.
 *         e2 may coincide with e if the squares are not needed.
 *
 *       ierr is set to
 *         zero       for normal return,
 *         n+i        if t(i,1)*t(i-1,3) is negative,
 *         -(3*n+i)   if t(i,1)*t(i-1,3) is zero with one factor
 *                    non-zero.  in this case, the eigenvectors of
 *                    the symmetric matrix are not simply related
 *                    to those of  t  and should not be sought.
 *
 *    questions and comments should be directed to burton s. garbow,
 *    mathematics and computer science div, argonne national laboratory
 *
 *    this version dated august 1983.
 */
static int figi(int nn, double *tt, double *dd, double *ee,
                double *e2)
{
  int ii;

  for (ii=0; ii<nn; ii++)
  {
      if (ii != 0)
      {
          e2[ii] = tt[3*ii]*tt[3*(ii-1)+2];

          if (e2[ii] < 0.0)
          {
              /* set error -- product of some pair of off-diagonal
                 elements is negative */
              return (nn + ii);
          }

          if (e2[ii] == 0.0 && (tt[3*ii] != 0.0 || tt[3*(ii-1)+2] != 0.0))
          {
              /* set error -- product of some pair of off-diagonal
                 elements is zero with one member non-zero */
              return (-1*(3*nn + ii));
          }

          ee[ii] = sqrt(e2[ii]);
      }

      dd[ii] = tt[3*ii+1];
  }

  return 0;
}


int gsl_sf_mathieu_a_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])
{
  unsigned int even_order = work->even_order, odd_order = work->odd_order,
      extra_values = work->extra_values, ii, jj;
  int status;
  double *tt = work->tt, *dd = work->dd, *ee = work->ee, *e2 = work->e2,
         *zz = work->zz, *aa = work->aa;
  gsl_matrix_view mat, evec;
  gsl_vector_view eval;
  gsl_eigen_symmv_workspace *wmat = work->wmat;
  
  if (order_max > work->size || order_max <= order_min || order_min < 0)
    {
      GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);
    }
  
  /* Convert the nonsymmetric tridiagonal matrix to a symmetric tridiagonal
     form. */

  tt[0] = 0.0;
  tt[1] = 0.0;
  tt[2] = qq;
  for (ii=1; ii<even_order-1; ii++)
  {
      tt[3*ii] = qq;
      tt[3*ii+1] = 4*ii*ii;
      tt[3*ii+2] = qq;
  }
  tt[3*even_order-3] = qq;
  tt[3*even_order-2] = 4*(even_order - 1)*(even_order - 1);
  tt[3*even_order-1] = 0.0;

  tt[3] *= 2;
  
  status = figi((signed int)even_order, tt, dd, ee, e2);

  if (status) 
    {
      GSL_ERROR("Internal error in tridiagonal Mathieu matrix", GSL_EFAILED);
    }

  /* Fill the period \pi matrix. */
  for (ii=0; ii<even_order*even_order; ii++)
      zz[ii] = 0.0;

  zz[0] = dd[0];
  zz[1] = ee[1];
  for (ii=1; ii<even_order-1; ii++)
  {
      zz[ii*even_order+ii-1] = ee[ii];
      zz[ii*even_order+ii] = dd[ii];
      zz[ii*even_order+ii+1] = ee[ii+1];
  }
  zz[even_order*(even_order-1)+even_order-2] = ee[even_order-1];
  zz[even_order*even_order-1] = dd[even_order-1];
  
  /* Compute (and sort) the eigenvalues of the matrix. */
  mat = gsl_matrix_view_array(zz, even_order, even_order);
  eval = gsl_vector_subvector(work->eval, 0, even_order);
  evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);
  gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
  gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
  
  for (ii=0; ii<even_order-extra_values; ii++)
      aa[2*ii] = gsl_vector_get(&eval.vector, ii);
  
  /* Fill the period 2\pi matrix. */
  for (ii=0; ii<odd_order*odd_order; ii++)
      zz[ii] = 0.0;
  for (ii=0; ii<odd_order; ii++)
      for (jj=0; jj<odd_order; jj++)
      {
          if (ii == jj)
              zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);
          else if (ii == jj + 1 || ii + 1 == jj)
              zz[ii*odd_order+jj] = qq;
      }
  zz[0] += qq;

  /* Compute (and sort) the eigenvalues of the matrix. */
  mat = gsl_matrix_view_array(zz, odd_order, odd_order);
  eval = gsl_vector_subvector(work->eval, 0, odd_order);
  evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);
  gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
  gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);

  for (ii=0; ii<odd_order-extra_values; ii++)
      aa[2*ii+1] = gsl_vector_get(&eval.vector, ii);

  for (ii = order_min ; ii <= order_max ; ii++)
    {
      result_array[ii - order_min] = aa[ii];
    }
  
  return GSL_SUCCESS;
}


int gsl_sf_mathieu_b_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])
{
  unsigned int even_order = work->even_order-1, odd_order = work->odd_order,
      extra_values = work->extra_values, ii, jj;
  double *zz = work->zz, *bb = work->bb;
  gsl_matrix_view mat, evec;
  gsl_vector_view eval;
  gsl_eigen_symmv_workspace *wmat = work->wmat;

  if (order_max > work->size || order_max <= order_min || order_min < 0)
    {
      GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);
    }

  /* Fill the period \pi matrix. */
  for (ii=0; ii<even_order*even_order; ii++)
      zz[ii] = 0.0;
  for (ii=0; ii<even_order; ii++)
      for (jj=0; jj<even_order; jj++)
      {
          if (ii == jj)
              zz[ii*even_order+jj] = 4*(ii + 1)*(ii + 1);
          else if (ii == jj + 1 || ii + 1 == jj)
              zz[ii*even_order+jj] = qq;
      }

  /* Compute (and sort) the eigenvalues of the matrix. */
  mat = gsl_matrix_view_array(zz, even_order, even_order);
  eval = gsl_vector_subvector(work->eval, 0, even_order);
  evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);
  gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
  gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);

  bb[0] = 0.0;
  for (ii=0; ii<even_order-extra_values; ii++)
      bb[2*(ii+1)] = gsl_vector_get(&eval.vector, ii);
  
  /* Fill the period 2\pi matrix. */
  for (ii=0; ii<odd_order*odd_order; ii++)
      zz[ii] = 0.0;
  for (ii=0; ii<odd_order; ii++)
      for (jj=0; jj<odd_order; jj++)
      {
          if (ii == jj)
              zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);
          else if (ii == jj + 1 || ii + 1 == jj)
              zz[ii*odd_order+jj] = qq;
      }

  zz[0] -= qq;

  /* Compute (and sort) the eigenvalues of the matrix. */
  mat = gsl_matrix_view_array(zz, odd_order, odd_order);
  eval = gsl_vector_subvector(work->eval, 0, odd_order);
  evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);
  gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
  gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
  
  for (ii=0; ii<odd_order-extra_values; ii++)
      bb[2*ii+1] = gsl_vector_get(&eval.vector, ii);  

  for (ii = order_min ; ii <= order_max ; ii++)
    {
      result_array[ii - order_min] = bb[ii];
    }

  return GSL_SUCCESS;
}