/* specfunc/mathieu_charv.c * * Copyright (C) 2002 Lowell Johnson * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at * your option) any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ /* Author: L. Johnson */ #include #include #include #include #include #include #include "mathieu.h" /* prototypes */ double solve_cubic(double c2, double c1, double c0); static double ceer(int order, double qq, double aa, int nterms) { double term, term1; int ii, n1; if (order == 0) term = 0.0; else { term = 2.0*qq*qq/aa; if (order != 2) { n1 = order/2 - 1; for (ii=0; ii= 0) { double t1 = rr + sqrt(ww); ss = fabs(t1)/t1*pow(fabs(t1), 1/3.); t1 = rr - sqrt(ww); tt = fabs(t1)/t1*pow(fabs(t1), 1/3.); } else { double theta = acos(rr/sqrt(-qq*qq*qq)); ss = 2*sqrt(-qq)*cos((theta + 4*M_PI)/3.); tt = 0.0; } return (ss + tt - c2/3); } /* Compute an initial approximation for the characteristic value. */ static double approx_c(int order, double qq) { double approx; double c0, c1, c2; if (order < 0) { printf("Error: Undefined order for Mathieu function.\n"); return 0.0; } switch (order) { case 0: if (qq <= 4) return (2 - sqrt(4 + 2*qq*qq)); /* Eqn. 31 */ else return asymptotic(order, qq); break; case 1: if (qq <= 4) return (5 + 0.5*(qq - sqrt(5*qq*qq - 16*qq + 64))); /* Eqn. 32 */ else return asymptotic(order, qq); break; case 2: if (qq <= 3) { c2 = -8.0; /* Eqn. 33 */ c1 = -48 - 3*qq*qq; c0 = 20*qq*qq; } else return asymptotic(order, qq); break; case 3: if (qq <= 6.25) { c2 = -qq - 8; /* Eqn. 34 */ c1 = 16*qq - 128 - 2*qq*qq; c0 = qq*qq*(qq + 8); } else return asymptotic(order, qq); break; default: if (order < 70) { if (1.7*order > 2*sqrt(qq)) { /* Eqn. 30 */ double n2 = (double)(order*order); double n22 = (double)((n2 - 1)*(n2 - 1)); double q2 = qq*qq; double q4 = q2*q2; approx = n2 + 0.5*q2/(n2 - 1); approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4)); approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/ (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9)); if (1.4*order < 2*sqrt(qq)) { approx += asymptotic(order, qq); approx *= 0.5; } } else approx = asymptotic(order, qq); return approx; } else return order*order; } /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */ approx = solve_cubic(c2, c1, c0); if ( approx < 0 && sqrt(qq) > 0.1*order ) return asymptotic(order-1, qq); else return (order*order + fabs(approx)); } static double approx_s(int order, double qq) { double approx; double c0, c1, c2; if (order < 1) { printf("Error: Undefined order for Mathieu function.\n"); return 0.0; } switch (order) { case 1: if (qq <= 4) return (5 - 0.5*(qq + sqrt(5*qq*qq + 16*qq + 64))); /* Eqn. 35 */ else return asymptotic(order-1, qq); break; case 2: if (qq <= 5) return (10 - sqrt(36 + qq*qq)); /* Eqn. 36 */ else return asymptotic(order-1, qq); break; case 3: if (qq <= 6.25) { c2 = qq - 8; /* Eqn. 37 */ c1 = -128 - 16*qq - 2*qq*qq; c0 = qq*qq*(8 - qq); } else return asymptotic(order-1, qq); break; default: if (order < 70) { if (1.7*order > 2*sqrt(qq)) { /* Eqn. 30 */ double n2 = (double)(order*order); double n22 = (double)((n2 - 1)*(n2 - 1)); double q2 = qq*qq; double q4 = q2*q2; approx = n2 + 0.5*q2/(n2 - 1); approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4)); approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/ (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9)); if (1.4*order < 2*sqrt(qq)) { approx += asymptotic(order-1, qq); approx *= 0.5; } } else approx = asymptotic(order-1, qq); return approx; } else return order*order; } /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */ approx = solve_cubic(c2, c1, c0); if ( approx < 0 && sqrt(qq) > 0.1*order ) return asymptotic(order-1, qq); else return (order*order + fabs(approx)); } int gsl_sf_mathieu_c_charv(int order, double qq, double *aa) { int even_odd, nterms = 12, ii, counter = 0, maxcount = 200; double a1, a2, fa, fa1, dela, aa_orig, da = 0.025; /* If the argument is 0, then the coefficient is simply the square of the order. */ if (qq == 0) { *aa = order*order; return GSL_SUCCESS; } /* Compute an initial approximation for the characteristic value. */ *aa = approx_c(order, qq); /* Save the original approximation for later comparison. */ aa_orig = *aa; even_odd = 0; if (order % 2 != 0) even_odd = 1; /* 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 = ceer(order, qq, a1, nterms); else fa1 = ceor(order, qq, a1, nterms); for (;;) { if (even_odd == 0) fa = ceer(order, qq, *aa, nterms); else fa = ceor(order, qq, *aa, nterms); a2 = a1; a1 = *aa; if (fa == fa1) break; *aa -= (*aa - a2)/(fa - fa1)*fa; dela = fabs(*aa - a2); if (dela < 1e-14) break; if (ii > 20) 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) break; if (*aa > aa_orig) *aa = aa_orig - da*counter; else *aa = aa_orig + da*counter; continue; } else break; } /* 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_s_charv(int order, double qq, double *aa) { int even_odd, nterms = 12, ii, counter = 0, maxcount = 200; double a1, a2, fa, fa1, dela, aa_orig, da = 0.025; /* If the argument is 0, then the coefficient is simply the square of the order. */ if (qq == 0) { *aa = order*order; return GSL_SUCCESS; } /* Compute an initial approximation for the characteristic value. */ *aa = approx_s(order, qq); /* Save the original approximation for later comparison. */ aa_orig = *aa; even_odd = 0; if (order % 2 != 0) even_odd = 1; /* 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) break; *aa -= (*aa - a2)/(fa - fa1)*fa; dela = fabs(*aa - a2); if (dela < 1e-14) break; if (ii > 20) 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) break; if (*aa > aa_orig) *aa = aa_orig - da*counter; else *aa = aa_orig + da*counter; continue; } else break; } /* 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 0.0) { ee[ii] = sqrt(e2[ii]); dd[ii] = tt[3*ii+1]; continue; } if (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)); } } dd[ii] = tt[3*ii+1]; } return 0; } double *gsl_sf_mathieu_c_charv_array(double qq, gsl_sf_mathieu_workspace *work) { 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; gsl_matrix_view mat, evec; gsl_vector_view eval; gsl_eigen_symmv_workspace *wmat = work->wmat; /* Allocate required storage. */ aa = (double *)malloc((work->size+1)*sizeof(double)); if (aa == NULL) { GSL_ERROR_NULL("Error allocating memory for the eigenvalues", GSL_ENOMEM); return NULL; } /* Convert the nonsymmetric tridiagonal matrix to a symmetric tridiagonal form. */ tt[0] = 0.0; tt[1] = 0.0; tt[2] = qq; for (ii=1; iieval, 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; iieval, 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; iieven_order-1, odd_order = work->odd_order, extra_values = work->extra_values, ii, jj; double *zz = work->zz, *bb; gsl_matrix_view mat, evec; gsl_vector_view eval; gsl_eigen_symmv_workspace *wmat = work->wmat; /* Allocate required storage. */ bb = (double *)malloc((work->size+1)*sizeof(double)); if (bb == NULL) { GSL_ERROR_NULL("Error allocating memory for the eigenvalues", GSL_ENOMEM); return NULL; } /* Fill the period \pi matrix. */ for (ii=0; iieval, 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; iieval, 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