1 /* specfunc/mathieu_charv.c
3 * Copyright (C) 2002 Lowell Johnson
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
20 /* Author: L. Johnson */
26 #include <gsl/gsl_math.h>
27 #include <gsl/gsl_eigen.h>
28 #include <gsl/gsl_errno.h>
29 #include <gsl/gsl_sf_mathieu.h>
33 static double solve_cubic(double c2, double c1, double c0);
36 static double ceer(int order, double qq, double aa, int nterms)
53 for (ii=0; ii<n1; ii++)
54 term = qq*qq/(aa - 4.0*(ii+1)*(ii+1) - term);
62 for (ii=0; ii<nterms; ii++)
64 (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
69 return (term + term1 - aa);
73 static double ceor(int order, double qq, double aa, int nterms)
79 n1 = (int)((float)order/2.0 - 0.5);
81 for (ii=0; ii<n1; ii++)
82 term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);
86 for (ii=0; ii<nterms; ii++)
88 (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
90 return (term + term1 - aa);
94 static double seer(int order, double qq, double aa, int nterms)
102 for (ii=0; ii<n1; ii++)
103 term = qq*qq/(aa - 4.0*(ii + 1)*(ii + 1) - term);
107 for (ii=0; ii<nterms; ii++)
109 (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
111 return (term + term1 - aa);
115 static double seor(int order, double qq, double aa, int nterms)
122 n1 = (int)((float)order/2.0 - 0.5);
123 for (ii=0; ii<n1; ii++)
124 term = qq*qq/(aa - (2.0*ii + 1.0)*(2.0*ii + 1.0) - term);
128 for (ii=0; ii<nterms; ii++)
130 (aa - (order + 2.0*(nterms - ii))*(order + 2.0*(nterms - ii)) - term1);
132 return (term + term1 - aa);
136 /*----------------------------------------------------------------------------
137 * Asymptotic and approximation routines for the characteristic value.
139 * Adapted from F.A. Alhargan's paper,
140 * "Algorithms for the Computation of All Mathieu Functions of Integer
141 * Orders," ACM Transactions on Mathematical Software, Vol. 26, No. 3,
142 * September 2000, pp. 390-407.
143 *--------------------------------------------------------------------------*/
144 static double asymptotic(int order, double qq)
147 double nn, n2, n4, n6;
148 double hh, ah, ah2, ah3, ah4, ah5;
151 /* Set up temporary variables to simplify the readability. */
164 /* Equation 38, p. 397 of Alhargan's paper. */
165 asymp = -2*qq + nn*hh - 0.125*(n2 + 1);
166 asymp -= 0.25*nn*( n2 + 3)/ah;
167 asymp -= 0.25* ( 5*n4 + 34*n2 + 9)/ah2;
168 asymp -= 0.25*nn*( 33*n4 + 410*n2 + 405)/ah3;
169 asymp -= ( 63*n6 + 1260*n4 + 2943*n2 + 486)/ah4;
170 asymp -= nn*(527*n6 + 15617*n4 + 69001*n2 + 41607)/ah5;
176 /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
177 static double solve_cubic(double c2, double c1, double c0)
179 double qq, rr, ww, ss, tt;
182 qq = (3*c1 - c2*c2)/9;
183 rr = (9*c2*c1 - 27*c0 - 2*c2*c2*c2)/54;
184 ww = qq*qq*qq + rr*rr;
188 double t1 = rr + sqrt(ww);
189 ss = fabs(t1)/t1*pow(fabs(t1), 1/3.);
191 tt = fabs(t1)/t1*pow(fabs(t1), 1/3.);
195 double theta = acos(rr/sqrt(-qq*qq*qq));
196 ss = 2*sqrt(-qq)*cos((theta + 4*M_PI)/3.);
200 return (ss + tt - c2/3);
204 /* Compute an initial approximation for the characteristic value. */
205 static double approx_c(int order, double qq)
213 GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);
220 return (2 - sqrt(4 + 2*qq*qq)); /* Eqn. 31 */
222 return asymptotic(order, qq);
227 return (5 + 0.5*(qq - sqrt(5*qq*qq - 16*qq + 64))); /* Eqn. 32 */
229 return asymptotic(order, qq);
235 c2 = -8.0; /* Eqn. 33 */
240 return asymptotic(order, qq);
246 c2 = -qq - 8; /* Eqn. 34 */
247 c1 = 16*qq - 128 - 2*qq*qq;
251 return asymptotic(order, qq);
257 if (1.7*order > 2*sqrt(qq))
260 double n2 = (double)(order*order);
261 double n22 = (double)((n2 - 1)*(n2 - 1));
264 approx = n2 + 0.5*q2/(n2 - 1);
265 approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));
266 approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/
267 (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));
268 if (1.4*order < 2*sqrt(qq))
270 approx += asymptotic(order, qq);
275 approx = asymptotic(order, qq);
283 /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
284 approx = solve_cubic(c2, c1, c0);
286 if ( approx < 0 && sqrt(qq) > 0.1*order )
287 return asymptotic(order-1, qq);
289 return (order*order + fabs(approx));
293 static double approx_s(int order, double qq)
301 GSL_ERROR_VAL("Undefined order for Mathieu function", GSL_EINVAL, 0.0);
308 return (5 - 0.5*(qq + sqrt(5*qq*qq + 16*qq + 64))); /* Eqn. 35 */
310 return asymptotic(order-1, qq);
315 return (10 - sqrt(36 + qq*qq)); /* Eqn. 36 */
317 return asymptotic(order-1, qq);
323 c2 = qq - 8; /* Eqn. 37 */
324 c1 = -128 - 16*qq - 2*qq*qq;
328 return asymptotic(order-1, qq);
334 if (1.7*order > 2*sqrt(qq))
337 double n2 = (double)(order*order);
338 double n22 = (double)((n2 - 1)*(n2 - 1));
341 approx = n2 + 0.5*q2/(n2 - 1);
342 approx += (5*n2 + 7)*q4/(32*n22*(n2 - 1)*(n2 - 4));
343 approx += (9*n2*n2 + 58*n2 + 29)*q4*q2/
344 (64*n22*n22*(n2 - 1)*(n2 - 4)*(n2 - 9));
345 if (1.4*order < 2*sqrt(qq))
347 approx += asymptotic(order-1, qq);
352 approx = asymptotic(order-1, qq);
360 /* Solve the cubic x^3 + c2*x^2 + c1*x + c0 = 0 */
361 approx = solve_cubic(c2, c1, c0);
363 if ( approx < 0 && sqrt(qq) > 0.1*order )
364 return asymptotic(order-1, qq);
366 return (order*order + fabs(approx));
370 int gsl_sf_mathieu_a(int order, double qq, gsl_sf_result *result)
372 int even_odd, nterms = 50, ii, counter = 0, maxcount = 200;
373 double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;
380 /* If the argument is 0, then the coefficient is simply the square of
384 result->val = order*order;
389 /* Use symmetry characteristics of the functions to handle cases with
390 negative order and/or argument q. See Abramowitz & Stegun, 20.8.3. */
396 return gsl_sf_mathieu_a(order, -qq, result);
398 return gsl_sf_mathieu_b(order, -qq, result);
401 /* Compute an initial approximation for the characteristic value. */
402 aa = approx_c(order, qq);
404 /* Save the original approximation for later comparison. */
407 /* Loop as long as the final value is not near the approximate value
408 (with a max limit to avoid potential infinite loop). */
409 while (counter < maxcount)
414 fa1 = ceer(order, qq, a1, nterms);
416 fa1 = ceor(order, qq, a1, nterms);
421 fa = ceer(order, qq, aa, nterms);
423 fa = ceor(order, qq, aa, nterms);
430 result->err = GSL_DBL_EPSILON;
433 aa -= (aa - a2)/(fa - fa1)*fa;
434 dela = fabs(aa - a2);
435 if (dela < GSL_DBL_EPSILON)
437 result->err = GSL_DBL_EPSILON;
449 /* If the solution found is not near the original approximation,
450 tweak the approximate value, and try again. */
451 if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)))
454 if (counter == maxcount)
456 result->err = fabs(aa - aa_orig);
460 aa = aa_orig - da*counter;
462 aa = aa_orig + da*counter;
472 /* If we went through the maximum number of retries and still didn't
473 find the solution, let us know. */
474 if (counter == maxcount)
476 GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);
483 int gsl_sf_mathieu_b(int order, double qq, gsl_sf_result *result)
485 int even_odd, nterms = 50, ii, counter = 0, maxcount = 200;
486 double a1, a2, fa, fa1, dela, aa_orig, da = 0.025, aa;
493 /* The order cannot be 0. */
496 GSL_ERROR("Characteristic value undefined for order 0", GSL_EFAILED);
499 /* If the argument is 0, then the coefficient is simply the square of
503 result->val = order*order;
508 /* Use symmetry characteristics of the functions to handle cases with
509 negative order and/or argument q. See Abramowitz & Stegun, 20.8.3. */
515 return gsl_sf_mathieu_b(order, -qq, result);
517 return gsl_sf_mathieu_a(order, -qq, result);
520 /* Compute an initial approximation for the characteristic value. */
521 aa = approx_s(order, qq);
523 /* Save the original approximation for later comparison. */
526 /* Loop as long as the final value is not near the approximate value
527 (with a max limit to avoid potential infinite loop). */
528 while (counter < maxcount)
533 fa1 = seer(order, qq, a1, nterms);
535 fa1 = seor(order, qq, a1, nterms);
540 fa = seer(order, qq, aa, nterms);
542 fa = seor(order, qq, aa, nterms);
549 result->err = GSL_DBL_EPSILON;
552 aa -= (aa - a2)/(fa - fa1)*fa;
553 dela = fabs(aa - a2);
556 result->err = GSL_DBL_EPSILON;
568 /* If the solution found is not near the original approximation,
569 tweak the approximate value, and try again. */
570 if (fabs(aa - aa_orig) > (3 + 0.01*order*fabs(aa_orig)))
573 if (counter == maxcount)
575 result->err = fabs(aa - aa_orig);
579 aa = aa_orig - da*counter;
581 aa = aa_orig + da*counter;
591 /* If we went through the maximum number of retries and still didn't
592 find the solution, let us know. */
593 if (counter == maxcount)
595 GSL_ERROR("Wrong characteristic Mathieu value", GSL_EFAILED);
602 /* Eigenvalue solutions for characteristic values below. */
605 /* figi.c converted from EISPACK Fortran FIGI.F.
607 * given a nonsymmetric tridiagonal matrix such that the products
608 * of corresponding pairs of off-diagonal elements are all
609 * non-negative, this subroutine reduces it to a symmetric
610 * tridiagonal matrix with the same eigenvalues. if, further,
611 * a zero product only occurs when both factors are zero,
612 * the reduced matrix is similar to the original matrix.
616 * n is the order of the matrix.
618 * t contains the input matrix. its subdiagonal is
619 * stored in the last n-1 positions of the first column,
620 * its diagonal in the n positions of the second column,
621 * and its superdiagonal in the first n-1 positions of
622 * the third column. t(1,1) and t(n,3) are arbitrary.
628 * d contains the diagonal elements of the symmetric matrix.
630 * e contains the subdiagonal elements of the symmetric
631 * matrix in its last n-1 positions. e(1) is not set.
633 * e2 contains the squares of the corresponding elements of e.
634 * e2 may coincide with e if the squares are not needed.
637 * zero for normal return,
638 * n+i if t(i,1)*t(i-1,3) is negative,
639 * -(3*n+i) if t(i,1)*t(i-1,3) is zero with one factor
640 * non-zero. in this case, the eigenvectors of
641 * the symmetric matrix are not simply related
642 * to those of t and should not be sought.
644 * questions and comments should be directed to burton s. garbow,
645 * mathematics and computer science div, argonne national laboratory
647 * this version dated august 1983.
649 static int figi(int nn, double *tt, double *dd, double *ee,
654 for (ii=0; ii<nn; ii++)
658 e2[ii] = tt[3*ii]*tt[3*(ii-1)+2];
662 /* set error -- product of some pair of off-diagonal
663 elements is negative */
667 if (e2[ii] == 0.0 && (tt[3*ii] != 0.0 || tt[3*(ii-1)+2] != 0.0))
669 /* set error -- product of some pair of off-diagonal
670 elements is zero with one member non-zero */
671 return (-1*(3*nn + ii));
674 ee[ii] = sqrt(e2[ii]);
684 int gsl_sf_mathieu_a_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])
686 unsigned int even_order = work->even_order, odd_order = work->odd_order,
687 extra_values = work->extra_values, ii, jj;
689 double *tt = work->tt, *dd = work->dd, *ee = work->ee, *e2 = work->e2,
690 *zz = work->zz, *aa = work->aa;
691 gsl_matrix_view mat, evec;
692 gsl_vector_view eval;
693 gsl_eigen_symmv_workspace *wmat = work->wmat;
695 if (order_max > work->size || order_max <= order_min || order_min < 0)
697 GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);
700 /* Convert the nonsymmetric tridiagonal matrix to a symmetric tridiagonal
706 for (ii=1; ii<even_order-1; ii++)
709 tt[3*ii+1] = 4*ii*ii;
712 tt[3*even_order-3] = qq;
713 tt[3*even_order-2] = 4*(even_order - 1)*(even_order - 1);
714 tt[3*even_order-1] = 0.0;
718 status = figi((signed int)even_order, tt, dd, ee, e2);
722 GSL_ERROR("Internal error in tridiagonal Mathieu matrix", GSL_EFAILED);
725 /* Fill the period \pi matrix. */
726 for (ii=0; ii<even_order*even_order; ii++)
731 for (ii=1; ii<even_order-1; ii++)
733 zz[ii*even_order+ii-1] = ee[ii];
734 zz[ii*even_order+ii] = dd[ii];
735 zz[ii*even_order+ii+1] = ee[ii+1];
737 zz[even_order*(even_order-1)+even_order-2] = ee[even_order-1];
738 zz[even_order*even_order-1] = dd[even_order-1];
740 /* Compute (and sort) the eigenvalues of the matrix. */
741 mat = gsl_matrix_view_array(zz, even_order, even_order);
742 eval = gsl_vector_subvector(work->eval, 0, even_order);
743 evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);
744 gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
745 gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
747 for (ii=0; ii<even_order-extra_values; ii++)
748 aa[2*ii] = gsl_vector_get(&eval.vector, ii);
750 /* Fill the period 2\pi matrix. */
751 for (ii=0; ii<odd_order*odd_order; ii++)
753 for (ii=0; ii<odd_order; ii++)
754 for (jj=0; jj<odd_order; jj++)
757 zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);
758 else if (ii == jj + 1 || ii + 1 == jj)
759 zz[ii*odd_order+jj] = qq;
763 /* Compute (and sort) the eigenvalues of the matrix. */
764 mat = gsl_matrix_view_array(zz, odd_order, odd_order);
765 eval = gsl_vector_subvector(work->eval, 0, odd_order);
766 evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);
767 gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
768 gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
770 for (ii=0; ii<odd_order-extra_values; ii++)
771 aa[2*ii+1] = gsl_vector_get(&eval.vector, ii);
773 for (ii = order_min ; ii <= order_max ; ii++)
775 result_array[ii - order_min] = aa[ii];
782 int gsl_sf_mathieu_b_array(int order_min, int order_max, double qq, gsl_sf_mathieu_workspace *work, double result_array[])
784 unsigned int even_order = work->even_order-1, odd_order = work->odd_order,
785 extra_values = work->extra_values, ii, jj;
786 double *zz = work->zz, *bb = work->bb;
787 gsl_matrix_view mat, evec;
788 gsl_vector_view eval;
789 gsl_eigen_symmv_workspace *wmat = work->wmat;
791 if (order_max > work->size || order_max <= order_min || order_min < 0)
793 GSL_ERROR ("invalid range [order_min,order_max]", GSL_EINVAL);
796 /* Fill the period \pi matrix. */
797 for (ii=0; ii<even_order*even_order; ii++)
799 for (ii=0; ii<even_order; ii++)
800 for (jj=0; jj<even_order; jj++)
803 zz[ii*even_order+jj] = 4*(ii + 1)*(ii + 1);
804 else if (ii == jj + 1 || ii + 1 == jj)
805 zz[ii*even_order+jj] = qq;
808 /* Compute (and sort) the eigenvalues of the matrix. */
809 mat = gsl_matrix_view_array(zz, even_order, even_order);
810 eval = gsl_vector_subvector(work->eval, 0, even_order);
811 evec = gsl_matrix_submatrix(work->evec, 0, 0, even_order, even_order);
812 gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
813 gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
816 for (ii=0; ii<even_order-extra_values; ii++)
817 bb[2*(ii+1)] = gsl_vector_get(&eval.vector, ii);
819 /* Fill the period 2\pi matrix. */
820 for (ii=0; ii<odd_order*odd_order; ii++)
822 for (ii=0; ii<odd_order; ii++)
823 for (jj=0; jj<odd_order; jj++)
826 zz[ii*odd_order+jj] = (2*ii + 1)*(2*ii + 1);
827 else if (ii == jj + 1 || ii + 1 == jj)
828 zz[ii*odd_order+jj] = qq;
833 /* Compute (and sort) the eigenvalues of the matrix. */
834 mat = gsl_matrix_view_array(zz, odd_order, odd_order);
835 eval = gsl_vector_subvector(work->eval, 0, odd_order);
836 evec = gsl_matrix_submatrix(work->evec, 0, 0, odd_order, odd_order);
837 gsl_eigen_symmv(&mat.matrix, &eval.vector, &evec.matrix, wmat);
838 gsl_eigen_symmv_sort(&eval.vector, &evec.matrix, GSL_EIGEN_SORT_VAL_ASC);
840 for (ii=0; ii<odd_order-extra_values; ii++)
841 bb[2*ii+1] = gsl_vector_get(&eval.vector, ii);
843 for (ii = order_min ; ii <= order_max ; ii++)
845 result_array[ii - order_min] = bb[ii];