# include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include # include int main ( ); void bessel_j0_values ( int *n_data, double *x, double *fx ); void bessel_j1_values ( int *n_data, double *x, double *fx ); void dawson_values ( int *n_data, double *x, double *fx ); void gsl_eigen_nonsymm_test ( ); void gsl_eigen_symm_test ( ); void gsl_fft_complex_test ( ); double gsl_integration_f ( double x, void *params ); void gsl_integration_qng_test ( ); void gsl_linalg_LU_test ( ); void gsl_midpoint_test ( ); void gsl_min_fminimizer_test ( ); double gsl_min_func ( double x, void *params ); double gsl_monte_f ( double x[], size_t dim_num, void *params ); void gsl_multiroot_fsolver_test ( ); void gsl_multiroot_print_state ( size_t iter, gsl_multiroot_fsolver *s ); void gsl_poly_complex_solve_test ( ); void gsl_poly_eval_test ( ); void gsl_qrng_niederreiter_2_test ( ); void gsl_qrng_sobol_test ( ); void gsl_sf_bessel_J0_test ( ); void gsl_sf_bessel_J1_test ( ); void gsl_sf_coupling_3j_test ( ); void gsl_sf_coupling_6j_test ( ); void gsl_sf_coupling_9j_test ( ); void gsl_sf_dawson_test ( ); void gsl_sf_hyperg_2F1_test ( ); void gsl_spmatrix_test ( ); void hyper_2f1_values ( int *n_data, double *a, double *b, double *c, double *x, double *fx ); void nine_j_values ( int *n_data, double *j1, double *j2, double *j3, double *j4, double *j5, double *j6, double *j7, double *j8, double *j9, double *fx ); int ode_func ( double t, const double y[], double f[], void *params ); int ode_jac ( double t, const double y[], double *dfdy, double dfdt[], void *params ); int rosenbrock_f ( const gsl_vector *x, void *params, gsl_vector *f ); void six_j_values ( int *n_data, double *j1, double *j2, double *j3, double *j4, double *j5, double *j6, double *fx ); void three_j_values ( int *n_data, double *j1, double *j2, double *j3, double *m1, double *m2, double *m3, double *fx ); void timestamp ( ); /* Rosenbrock global parameter structure. */ struct rparams { double a; double b; }; /******************************************************************************/ int main ( ) /******************************************************************************/ /* Purpose: gsl_test() tests gsl(). Licensing: This code is distributed under the MIT license. Modified: 28 September 2024 Author: John Burkardt */ { timestamp ( ); printf ( "\n" ); printf ( "gsl_test():\n" ); printf ( " C version\n" ); printf ( " Test gsl(), the GNU Scientific Library (GSL).\n" ); gsl_eigen_nonsymm_test ( ); gsl_eigen_symm_test ( ); gsl_fft_complex_test ( ); gsl_integration_qng_test ( ); gsl_linalg_LU_test ( ); gsl_midpoint_test ( ); gsl_min_fminimizer_test ( ); gsl_multiroot_fsolver_test ( ); gsl_poly_complex_solve_test ( ); gsl_poly_eval_test ( ); gsl_qrng_niederreiter_2_test ( ); gsl_qrng_sobol_test ( ); gsl_sf_bessel_J0_test ( ); gsl_sf_bessel_J1_test ( ); gsl_sf_coupling_3j_test ( ); gsl_sf_coupling_6j_test ( ); gsl_sf_coupling_9j_test ( ); gsl_sf_dawson_test ( ); gsl_sf_hyperg_2F1_test ( ); gsl_spmatrix_test ( ); /* Terminate. */ printf ( "\n" ); printf ( "gsl_test():\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ void bessel_j0_values ( int *n_data, double *x, double *fx ) /******************************************************************************/ /* Purpose: bessel_j0_values() returns some values of the J0 Bessel function. Discussion: In Mathematica, the function can be evaluated by: BesselJ[0,x] Licensing: This code is distributed under the MIT license. Modified: 10 August 2004 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. The routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 21 static double fx_vec[N_MAX] = { -0.1775967713143383E+00, -0.3971498098638474E+00, -0.2600519549019334E+00, 0.2238907791412357E+00, 0.7651976865579666E+00, 0.1000000000000000E+01, 0.7651976865579666E+00, 0.2238907791412357E+00, -0.2600519549019334E+00, -0.3971498098638474E+00, -0.1775967713143383E+00, 0.1506452572509969E+00, 0.3000792705195556E+00, 0.1716508071375539E+00, -0.9033361118287613E-01, -0.2459357644513483E+00, -0.1711903004071961E+00, 0.4768931079683354E-01, 0.2069261023770678E+00, 0.1710734761104587E+00, -0.1422447282678077E-01 }; static double x_vec[N_MAX] = { -5.0E+00, -4.0E+00, -3.0E+00, -2.0E+00, -1.0E+00, 0.0E+00, 1.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, 7.0E+00, 8.0E+00, 9.0E+00, 10.0E+00, 11.0E+00, 12.0E+00, 13.0E+00, 14.0E+00, 15.0E+00 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *x = 0.0; *fx = 0.0; } else { *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ void bessel_j1_values ( int *n_data, double *x, double *fx ) /******************************************************************************/ /* Purpose: bessel_j1_values() returns some values of the J1 Bessel function. Discussion: In Mathematica, the function can be evaluated by: BesselJ[1,x] Licensing: This code is distributed under the MIT license. Modified: 12 August 2004 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. The routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 21 static double fx_vec[N_MAX] = { 0.3275791375914652E+00, 0.6604332802354914E-01, -0.3390589585259365E+00, -0.5767248077568734E+00, -0.4400505857449335E+00, 0.0000000000000000E+00, 0.4400505857449335E+00, 0.5767248077568734E+00, 0.3390589585259365E+00, -0.6604332802354914E-01, -0.3275791375914652E+00, -0.2766838581275656E+00, -0.4682823482345833E-02, 0.2346363468539146E+00, 0.2453117865733253E+00, 0.4347274616886144E-01, -0.1767852989567215E+00, -0.2234471044906276E+00, -0.7031805212177837E-01, 0.1333751546987933E+00, 0.2051040386135228E+00 }; static double x_vec[N_MAX] = { -5.0E+00, -4.0E+00, -3.0E+00, -2.0E+00, -1.0E+00, 0.0E+00, 1.0E+00, 2.0E+00, 3.0E+00, 4.0E+00, 5.0E+00, 6.0E+00, 7.0E+00, 8.0E+00, 9.0E+00, 10.0E+00, 11.0E+00, 12.0E+00, 13.0E+00, 14.0E+00, 15.0E+00 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *x = 0.0; *fx = 0.0; } else { *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ void dawson_values ( int *n_data, double *x, double *fx ) /******************************************************************************/ /* Purpose: dawson_values() returns some values of Dawson's integral. Discussion: The definition of Dawson's integral is D(X) = exp ( -X * X ) * Integral ( 0 <= Y <= X ) exp ( Y * Y ) dY Dawson's integral has a maximum at roughly X = 0.9241388730 In Mathematica, the function can be evaluated by: Sqrt[Pi] * Exp[-x^2] * I * Erf[I*x] / 2 Licensing: This code is distributed under the MIT license. Modified: 18 August 2004 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Eric Weisstein, CRC Concise Encyclopedia of Mathematics, CRC Press, 1998. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. The routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 21 static double fx_vec[N_MAX] = { 0.0000000000000000E+00, 0.9933599239785286E-01, 0.1947510333680280E+00, 0.2826316650213119E+00, 0.3599434819348881E+00, 0.4244363835020223E+00, 0.4747632036629779E+00, 0.5105040575592318E+00, 0.5321017070563654E+00, 0.5407243187262987E+00, 0.5380795069127684E+00, 0.5262066799705525E+00, 0.5072734964077396E+00, 0.4833975173848241E+00, 0.4565072375268973E+00, 0.4282490710853986E+00, 0.3999398943230814E+00, 0.3725593489740788E+00, 0.3467727691148722E+00, 0.3229743193228178E+00, 0.3013403889237920E+00 }; static double x_vec[N_MAX] = { 0.0E+00, 0.1E+00, 0.2E+00, 0.3E+00, 0.4E+00, 0.5E+00, 0.6E+00, 0.7E+00, 0.8E+00, 0.9E+00, 1.0E+00, 1.1E+00, 1.2E+00, 1.3E+00, 1.4E+00, 1.5E+00, 1.6E+00, 1.7E+00, 1.8E+00, 1.9E+00, 2.0E+00 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *x = 0.0; *fx = 0.0; } else { *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ void gsl_eigen_nonsymm_test ( ) /******************************************************************************/ /* Purpose: gsl_eigen_nonsymm_test() tests gsl_eigen_nonsymm(). Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Reference: Mark Gelassi, Jim Davies, James Tyler, Bryan Gough, Reid Priedhorsky, Gerard Jungman, Michael Booth, Fabrice Rossi, GNU Scientific Library Reference Manual. */ { double data[] = { -1.0, 1.0, -1.0, 1.0, -8.0, 4.0, -2.0, 1.0, 27.0, 9.0, 3.0, 1.0, 64.0, 16.0, 4.0, 1.0 }; /* Convert C array to GSL matrix. */ gsl_matrix_view m = gsl_matrix_view_array ( data, 4, 4 ); /* Allocate GSL matrix and vector space for eigenvectors and eigenvalues. */ gsl_matrix_complex *evec = gsl_matrix_complex_alloc ( 4, 4 ); gsl_vector_complex *eval = gsl_vector_complex_alloc ( 4 ); /* Set up workspace W. */ gsl_eigen_nonsymmv_workspace *w = gsl_eigen_nonsymmv_alloc ( 4 ); printf ( "\n" ); printf ( "gsl_eigen_nonsymm_test():\n" ); printf ( " gsl_eigen_nonsymm() computes the eigenvalues and eigenvectors\n" ); printf ( " of a nonsymmetric matrix.\n" ); /* Compute the eigenvalues and eigenvectors. */ gsl_eigen_nonsymmv ( &m.matrix, eval, evec, w ); /* Free the workspace memory. */ gsl_eigen_nonsymmv_free ( w ); /* Sort the results. */ gsl_eigen_nonsymmv_sort ( eval, evec, GSL_EIGEN_SORT_ABS_DESC ); /* Print the results. */ for ( int i = 0; i < 4; i++ ) { gsl_complex eval_i = gsl_vector_complex_get ( eval, i ); gsl_vector_complex_view evec_i = gsl_matrix_complex_column ( evec, i ); printf ( "\n" ); printf ( " Eigenvalue(%d) = %g + %gi\n", i, GSL_REAL ( eval_i ), GSL_IMAG ( eval_i ) ); printf ( "\n" ); printf ( " Eigenvector(%d) = \n", i ); for ( int j = 0; j < 4; j++ ) { gsl_complex z = gsl_vector_complex_get ( &evec_i.vector, j ); printf ( " %g + %gi\n", GSL_REAL ( z ), GSL_IMAG ( z ) ); } } /* Free memory. */ gsl_vector_complex_free ( eval ); gsl_matrix_complex_free ( evec ); return; } /******************************************************************************/ void gsl_eigen_symm_test ( ) /******************************************************************************/ /* Purpose: gsl_eigen_symm_test() tests gsl_eigen_symm(). Licensing: This code is distributed under the MIT license. Modified: 07 June 2024 Reference: Mark Gelassi, Jim Davies, James Tyler, Bryan Gough, Reid Priedhorsky, Gerard Jungman, Michael Booth, Fabrice Rossi, GNU Scientific Library Reference Manual. */ { gsl_matrix *A; int i; int j; int n; double value; n = 5; A = gsl_matrix_alloc ( n, n ); for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { if ( j == i + 1 ) { value = sqrt ( ( double ) ( ( i + 1 ) * ( n - i - 1 ) ) ); } else if ( i == j + 1 ) { value = sqrt ( ( double ) ( ( j + 1 ) * ( n - j - 1 ) ) ); } else { value = 0.0; } gsl_matrix_set ( A, i, j, value ); } } /* Allocate space for eigenvectors and eigenvalues. */ gsl_matrix *evec = gsl_matrix_alloc ( n, n ); gsl_vector *eval = gsl_vector_alloc ( n ); /* Set up workspace W. */ gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc ( 4 * n ); printf ( "\n" ); printf ( "gsl_eigen_symm_test():\n" ); printf ( " gsl_eigen_symm() computes the eigenvalues and eigenvectors\n" ); printf ( " of a symmetric matrix.\n" ); /* Compute the eigenvalues and eigenvectors. */ gsl_eigen_symmv ( A, eval, evec, w ); /* Sort the results in ascending order. */ gsl_eigen_symmv_sort ( eval, evec, GSL_EIGEN_SORT_ABS_ASC ); /* Print the eigenvalues. */ printf ( "\n" ); printf ( " Eigenvalues:\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { printf ( " %d: %g\n", i, gsl_vector_get ( eval, i ) ); } /* Print the eigenvectors. */ printf ( "\n" ); printf ( " Eigenvectors:\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { for ( j = 0; j < n; j++ ) { printf ( " %12.4g", gsl_matrix_get ( evec, i, j ) ); } printf ( "\n" ); } /* Free memory. */ gsl_eigen_symmv_free ( w ); return; } /******************************************************************************/ void gsl_fft_complex_test ( ) /******************************************************************************/ /* Purpose: gsl_fft_complex_test() tests gsl_fft_complex(). Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Reference: Mark Gelassi, Jim Davies, James Tyler, Bryan Gough, Reid Priedhorsky, Gerard Jungman, Michael Booth, Fabrice Rossi, GNU Scientific Library Reference Manual. */ { double data[2*128]; int i; printf ( "\n" ); printf ( "gsl_fft_complex_test():\n" ); printf ( " gsl_fft_complex() computes the fast Fourier transform\n" ); printf ( " of a complex vector of data.\n" ); for ( i = 0; i < 128; i++ ) { data[0+i*2] = 0.0; data[1+i*2] = 0.0; } data[0+0*2] = 1.0; /* Insert a pulse. Make it symmetric so the transform will actually be real. */ for ( i = 1; i <= 10; i++ ) { data[0+i*2] = 1.0; data[0+(128-i)*2] = 1.0; } /* Print the input. */ printf ( "\n" ); printf ( " Input data (a symmetric pulse):\n" ); printf ( "\n" ); for ( i = 0; i < 128; i++ ) { printf ( "%3d %e %e\n", i, data[0+i*2], data[1+i*2] ); } /* Compute the FFT. */ gsl_fft_complex_radix2_forward ( data, 1, 128 ); /* Print the scaled output. */ printf ( "\n" ); printf ( " Output data:\n" ); printf ( "\n" ); for ( i = 0; i < 128; i++ ) { printf ( "%3d %e %e\n", i, data[0+i*2] / sqrt ( 128.0 ), data[1+i*2] / sqrt ( 128.0 ) ); } return; } /******************************************************************************/ double gsl_integration_f ( double x, void *params ) /******************************************************************************/ /* Purpose: gsl_integration_f() evaluates the integrand for quadrature test. Discussion: Integral ( 0 <= x <= 1 ) 2/(2+sin(10*pi*x)) = 2/sqrt(3). Licensing: This code is distributed under the MIT license. Modified: 14 February 2022 Author: Original C version by GSL reference manual. This version by John Burkardt. Reference: Mark Gelassi, Jim Davies, James Tyler, Bryan Gough, Reid Priedhorsky, Gerard Jungman, Michael Booth, Fabrice Rossi, GNU Scientific Library Reference Manual. Input: double *X, the evaluation point. void *PARAMS, optional parameters. Output: double GSL_INTEGRATION_F, the value of the integrand at X. */ { double value; value = 2.0 / ( 2.0 + sin ( 10.0 * M_PI * x ) ); return value; } /******************************************************************************/ void gsl_integration_qng_test ( ) /******************************************************************************/ /* Purpose: gsl_integration_qng_test() tests gsl_integration_qng(). Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: Original C version by GSL reference manual. This version by John Burkardt. Reference: Mark Gelassi, Jim Davies, James Tyler, Bryan Gough, Reid Priedhorsky, Gerard Jungman, Michael Booth, Fabrice Rossi, GNU Scientific Library Reference Manual. */ { double a; double abserr; double abserr_exact; double b; double epsabs; double epsrel; int error_code; gsl_function F; size_t neval; double result; double result_exact; printf ( "\n" ); printf ( "gsl_integration_qng_test():\n" ); printf ( " gsl_integration_qng() uses the Gauss-Kronrod rule to estimate\n" ); printf ( " the integral of f(x) over [a,b].\n" ); F.function = &gsl_integration_f; F.params = NULL; a = 0.0; b = 1.0; /* If we reduce EPSABS much more, then the code returns an error. */ epsabs = 1.0E-01; epsrel = 0.0E+00; error_code = gsl_integration_qng ( &F, a, b, epsabs, epsrel, &result, &abserr, &neval ); if ( error_code != 0 ) { printf ( "\n" ); printf ( " gsl_integration_qng() returned error code = %d\n", error_code ); } else { printf ( "\n" ); printf ( " Integral estimate = %g\n", result ); printf ( " Error estimate = %g\n", abserr ); printf ( " Number of function evaluations was %ld\n", neval ); result_exact = 2.0 / sqrt ( 3.0 ); abserr_exact = fabs ( result - result_exact ); printf ( " Exact integral = %g\n", result_exact ); printf ( " Exact error = %g\n", abserr_exact ); } return; } /******************************************************************************/ void gsl_linalg_LU_test ( ) /******************************************************************************/ /* Purpose: gsl_linalg_LU_test() tests gsl_linalg_LU(). Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Reference: Mark Gelassi, Jim Davies, James Tyler, Bryan Gough, Reid Priedhorsky, Gerard Jungman, Michael Booth, Fabrice Rossi, GNU Scientific Library Reference Manual. */ { double a_data[] = { 0.18, 0.60, 0.57, 0.96, 0.41, 0.24, 0.99, 0.58, 0.14, 0.30, 0.97, 0.66, 0.51, 0.13, 0.19, 0.85 }; double b_data[] = { 1.0, 2.0, 3.0, 4.0 }; /* Create GSL_MATRIX M. */ gsl_matrix_view m = gsl_matrix_view_array ( a_data, 4, 4 ); /* Create GSL_VECtOR B. */ gsl_vector_view b = gsl_vector_view_array ( b_data, 4 ); /* Allocate GSL_VECTOR X (solution). */ gsl_vector *x = gsl_vector_alloc ( 4 ); int s; /* Allocate GSL_PERMUTATION P (pivot vector). */ gsl_permutation *p = gsl_permutation_alloc ( 4 ); printf ( "\n" ); printf ( "gsl_linalg_LU_test():\n" ); printf ( " gsl_linalg_LU_decomp() computes LU factors of a matrix A.\n" ); printf ( " gsl_linalg_LU_solve() solves linear systems A*x=b.\n" ); /* LU factor the matrix M. */ gsl_linalg_LU_decomp ( &m.matrix, p, &s ); /* Solve M*x=b. */ gsl_linalg_LU_solve ( &m.matrix, p, &b.vector, x ); /* Display the solution. */ printf ( "\n" ); printf ( " Computed solution vector X: \n" ); printf ( "\n" ); gsl_vector_fprintf ( stdout, x, "%g" ); /* Free memory. */ gsl_permutation_free ( p ); gsl_vector_free ( x ); return; } /******************************************************************************/ void gsl_midpoint_test ( ) /******************************************************************************/ /* Purpose: gsl_midpoint_test() demonstrates the GSL midpoint ODE solver. Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: Original C version by GSL reference manual. This version by John Burkardt. */ { double dt; FILE *data; char *data_filename = "gsl_midpoint_data.txt"; double epsabs = 1e-6; double epsrel = 0.0; double h_start = 1e-6; int i; double mu = 10.0; int n = 100; double t = 0.0; double t0 = 0.0; double tstop = 100.0; double y[2] = { 1.0, 0.0 }; /* Define "sys", the definition of the ODE data. */ gsl_odeiv2_system sys = { ode_func, ode_jac, 2, &mu }; /* Define "d", a pointer to the ODE driver to be used. Here, "gsl_odeiv2_step_rk2imp" specifies we want to use the GSL midpoint driver. */ gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new ( &sys, gsl_odeiv2_step_rk2imp, h_start, epsabs, epsrel ); printf ( "\n" ); printf ( "gsl_midpoint_test():\n" ); printf ( " Test the GSL implicit midpoint method ODE solver.\n" ); printf ( "\n" ); printf ( " parameters:\n" ); printf ( " mu = %g\n", mu ); printf ( " t0 = %g\n", t0 ); printf ( " y0 = (%g,%g)\n", y[0], y[1] ); printf ( " tstop = %g\n", tstop ); printf ( " n = %d\n", n ); data = fopen ( data_filename, "wt" ); fprintf ( data, " %g %g %g\n", t0, y[0], y[1] ); dt = ( tstop - t0 ) / n; for ( i = 1; i <= n; i++ ) { double ti = i * dt; int status = gsl_odeiv2_driver_apply ( d, &t, ti, y ); if ( status != GSL_SUCCESS ) { printf ( "error, return value = %d\n", status ); break; } printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); fprintf ( data, " %g %g %g\n", t, y[0], y[1] ); } fclose ( data ); printf ( "\n" ); printf ( " Solution data stored in '%s'.\n", data_filename ); /* Free memory. */ gsl_odeiv2_driver_free ( d ); return; } /******************************************************************************/ void gsl_min_fminimizer_test ( ) /******************************************************************************/ /* Purpose: gsl_min_fminimizer_test() demonstrates the GSL function minimizer. Licensing: This code is distributed under the MIT license. Modified: 28 September 2024 Author: Original C version by GSL reference manual. This version by John Burkardt. */ { double a = 0.0; double b = 6.0; gsl_function F; int iter = 0; double m = 2.0; double m_expected = M_PI; int max_iter = 100; gsl_min_fminimizer *s; int status; const gsl_min_fminimizer_type *T; printf ( "\n" ); printf ( "gsl_min_fminimizer_test():\n" ); printf ( " Find minimizer of f(x) within interval [a,b].\n" ); F.function = &gsl_min_func; F.params = 0; T = gsl_min_fminimizer_brent; s = gsl_min_fminimizer_alloc (T); gsl_min_fminimizer_set (s, &F, m, a, b); printf ( "using %s method\n", gsl_min_fminimizer_name ( s ) ); printf ("%5s [%9s, %9s] %9s %10s %9s\n", "iter", "lower", "upper", "minimizer", "err", "err(est)"); printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n", iter, a, b, m, m - m_expected, b - a); do { iter++; status = gsl_min_fminimizer_iterate (s); m = gsl_min_fminimizer_x_minimum (s); a = gsl_min_fminimizer_x_lower (s); b = gsl_min_fminimizer_x_upper (s); status = gsl_min_test_interval (a, b, 0.001, 0.0); if ( status == GSL_SUCCESS ) { printf ("Converged:\n"); } printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n", iter, a, b, m, m - m_expected, b - a); } while ( status == GSL_CONTINUE && iter < max_iter ); /* Free memory. */ gsl_min_fminimizer_free ( s ); return; } /******************************************************************************/ double gsl_min_func ( double x, void *params ) /******************************************************************************/ /* Purpose: gsl_min_func() evaluates the function whose minimizer is desired. Licensing: This code is distributed under the MIT license. Modified: 28 September 2024 Author: Original C version by GSL reference manual. This version by John Burkardt. Input: double x: the argument. void *params: optional parameters, not used here. Output: double value: the value of f(x). */ { double value; /* avoid unused parameter warning. */ ( void ) ( params ); value = cos ( x ) + 1.0; return value; } /******************************************************************************/ double gsl_monte_f ( double x[], size_t dim_num, void *params ) /******************************************************************************/ /* Purpose: gsl_monte_f() evaluates the integrand for a Monte Carlo test. Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: Original C version by GSL reference manual. This version by John Burkardt. Reference: Mark Gelassi, Jim Davies, James Tyler, Bryan Gough, Reid Priedhorsky, Gerard Jungman, Michael Booth, Fabrice Rossi, GNU Scientific Library Reference Manual. Input: double x[], the evaluation point. size_t dim, the spatial dimension. void *params, optional parameters. Output: double gsl_monte_f, the value of the integrand at X. */ { double c; double den; size_t dim; double value; c = 1.0 / pow ( 2.0 * M_PI, dim_num ); den = 1.0; for ( dim = 0; dim < dim_num; dim++ ) { den = den * cos ( x[dim] ); } den = 1.0 - den; if ( den == 0.0 ) { value = 100000.0; } else { value = c / den; } return value; } /******************************************************************************/ void gsl_multiroot_fsolver_test ( ) /******************************************************************************/ /* Purpose: gsl_multiroot_fsolver_test() tests the N-dimensional root finder. Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: Original C version by GSL reference manual. This version by John Burkardt. */ { size_t iter = 0; const size_t n = 2; struct rparams p = { 1.0, 10.0 }; gsl_multiroot_fsolver *s; int status; const gsl_multiroot_fsolver_type *T; gsl_vector *x; double x_init[2] = { -10.0, -5.0 }; printf ( "\n" ); printf ( "gsl_multiroot_fsolver_test():\n" ); printf ( " gsl_multiroot_fsolver() finds a root of a\n" ); printf ( " set of nonlinear equations.\n" ); printf ( "\n" ); printf ( " In this case, we have two functions in two unknowns,\n" ); printf ( " and the only root is X = (1,1).\n" ); printf ( "\n" ); gsl_multiroot_function f = { &rosenbrock_f, n, &p }; x = gsl_vector_alloc ( n ); gsl_vector_set ( x, 0, x_init[0] ); gsl_vector_set ( x, 1, x_init[1] ); T = gsl_multiroot_fsolver_hybrids; s = gsl_multiroot_fsolver_alloc ( T, 2 ); gsl_multiroot_fsolver_set ( s, &f, x ); gsl_multiroot_print_state ( iter, s ); do { iter++; status = gsl_multiroot_fsolver_iterate ( s ); gsl_multiroot_print_state ( iter, s ); if ( status ) { break; } status = gsl_multiroot_test_residual ( s->f, 1.0E-07 ); } while ( status == GSL_CONTINUE && iter < 1000 ); printf ( " status = %s\n", gsl_strerror ( status ) ); /* Free memory. */ gsl_multiroot_fsolver_free ( s ); gsl_vector_free ( x ); return; } /******************************************************************************/ void gsl_multiroot_print_state ( size_t iter, gsl_multiroot_fsolver *s ) /******************************************************************************/ /* Purpose: gsl_multiroot_print_state() prints the state of the root finding iteration. Licensing: This code is distributed under the MIT license. Modified: 05 August 2005 Author: Original C version by GSL reference manual. This version by John Burkardt. Input: size_t iter: the iteration counter. gsl_multiroot_fsolver *s, the solver state. */ { printf ( " iter = %3lu x = %.3f %.3f f(x) = %.3e %.3e\n", iter, gsl_vector_get ( s->x, 0 ), gsl_vector_get ( s->x, 1 ), gsl_vector_get ( s->f, 0 ), gsl_vector_get ( s->f, 1 ) ); return; } /******************************************************************************/ void gsl_poly_complex_solve_test ( ) /******************************************************************************/ /* Purpose: gsl_poly_complex_solve() tests gsl_poly_complex_solve(). Licensing: This code is distributed under the MIT license. Modified: 07 June 2024 Author: Original C version by GSL reference manual. This version by John Burkardt. */ { int n = 5; int i; /* Space for n+1 coefficients. Here P(x) = -1 + x^5 */ double a[6] = { -1, 0, 0, 0, 0, 1 }; /* Space for real and imaginary parts of n roots. */ double z[2*n]; printf ( "\n" ); printf ( "gsl_poly_complex_solve_test():\n" ); printf ( " gsl_poly_complex_solve() finds all roots of a polynomial.\n" ); printf ( " Here, p(x) = x^5-1.\n" ); gsl_poly_complex_workspace *w = gsl_poly_complex_workspace_alloc ( n + 1 ); gsl_poly_complex_solve ( a, n + 1, w, z ); gsl_poly_complex_workspace_free ( w ); printf ( "\n" ); printf ( " Computed roots:\n" ); printf ( "\n" ); for ( i = 0; i < n; i++ ) { printf ( "z%d = %+.18f %+.18f\n", i, z[2*i], z[2*i+1]); } return; } /******************************************************************************/ void gsl_poly_eval_test ( ) /******************************************************************************/ /* Purpose: gsl_poly_eval_test() tests gsl_poly_eval(). Licensing: This code is distributed under the MIT license. Modified: 14 February 2022 Author: Original C version by GSL reference manual. This version by John Burkardt. */ { double c[4] = {-18.0, 27.0, -10.0, 1.0 }; int i; int len = 4; double p; double x; printf ( "\n" ); printf ( "gsl_poly_eval_test():\n" ); printf ( " gsl_poly_eval() evaluates a real polynomial.\n" ); printf ( " p(x) = (x-1)*(x-3)*(x-6)\n" ); printf ( " = x^3 - 10x^2 + 27x - 18.\n" ); printf ( "\n" ); printf ( " x p(x)\n" ); printf ( "\n" ); for ( i = 0; i <= 20; i++ ) { x = ( double ) i / 2.0; p = gsl_poly_eval ( c, len, x ); printf ( " %8.3g %8.3g\n", x, p ); } return; } /******************************************************************************/ void gsl_qrng_niederreiter_2_test ( ) /******************************************************************************/ /* Purpose: gsl_qrng_niederreiter_2_test() tests gsl_qrng_niederreiter_2(). Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: Original C version by GSL reference manual. This version by John Burkardt. */ { # define M 2 # define N 25 int i; int j; double v[M]; gsl_qrng * q = gsl_qrng_alloc ( gsl_qrng_niederreiter_2, M ); printf ( "\n" ); printf ( "gsl_qrng_niederreiter_2_test():\n" ); printf ( " gsl_qrng_alloc sets aside space for a sequence;\n" ); printf ( " gsl_qrng_niederreiter_2 requests the Niederreiter_2 sequence;\n" ); printf ( " gsl_qrng_get gets the next entry of the requested sequence;\n" ); printf ( "\n" ); printf ( " Determine the first %d points of the Niederreiter2\n", N ); printf ( " quasi-random sequence in %d dimensions.\n", M ); printf ( "\n" ); printf ( " I X(I)\n" ); printf ( "\n" ); for ( i = 0; i < N; i++ ) { gsl_qrng_get ( q, v ); printf ( "%6d ", i ); for ( j = 0; j < M; j++ ) { printf ( "%12.6g", v[j] ); } printf ( "\n" ); } return; # undef M # undef N } /******************************************************************************/ void gsl_qrng_sobol_test ( ) /******************************************************************************/ /* Purpose: gsl_qrng_sobol_test() tests the GSL Sobol sequence routine. Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: Original C version by GSL reference manual. This version by John Burkardt. */ { # define M 2 # define N 25 int i; int j; double v[M]; gsl_qrng * q = gsl_qrng_alloc ( gsl_qrng_sobol, M ); printf ( "\n" ); printf ( "gsl_qrng_sobol_test():\n" ); printf ( " gsl_qrng_alloc() sets aside space for a sequence;\n" ); printf ( " gsl_qrng_sobol() requests the Sobol sequence;\n" ); printf ( " gsl_qrng_get() gets the next entry of the requested sequence;\n" ); printf ( "\n" ); printf ( " Determine the first %d points of the Sobol\n", N ); printf ( " quasi-random sequence in %d dimensions.\n", M ); printf ( "\n" ); printf ( " I X(I)\n" ); printf ( "\n" ); for ( i = 0; i < N; i++ ) { gsl_qrng_get ( q, v ); printf ( "%6d ", i ); for ( j = 0; j < M; j++ ) { printf ( "%12.6g", v[j] ); } printf ( "\n" ); } return; # undef M # undef N } /******************************************************************************/ void gsl_sf_bessel_J0_test ( ) /******************************************************************************/ /* Purpose: gsl_sf_bessel_J0_test() tests gsl_sf_bessel_J0(). Licensing: This code is distributed under the MIT license. Modified: 14 February 2022 Author: John Burkardt */ { double fx1; double fx2; int n_data; double x; printf ( "\n" ); printf ( "gsl_sf_bessel_J0_test():\n" ); printf ( " gsl_sf_bessel_J0() evaluates the J0 Bessel function.\n" ); printf ( "\n" ); printf ( " X Exact Computed\n" ); printf ( "\n" ); n_data = 0; for ( ; ; ) { bessel_j0_values ( &n_data, &x, &fx1 ); if ( n_data <= 0 ) { break; } fx2 = gsl_sf_bessel_J0 ( x ); printf ( " %12.6g %12.6g %12.6g\n", x, fx1, fx2 ); } return; } /******************************************************************************/ void gsl_sf_bessel_J1_test ( ) /******************************************************************************/ /* Purpose: gsl_sf_bessel_J1_test() tests gsl_sf_bessel_J1(). Licensing: This code is distributed under the MIT license. Modified: 14 February 2022 Author: John Burkardt */ { double fx1; double fx2; int n_data; double x; printf ( "\n" ); printf ( "gsl_sf_bessel_J1_test():\n" ); printf ( " gsl_sf_bessel_J1() evaluates the J1 Bessel function.\n" ); printf ( "\n" ); printf ( " X Exact Computed\n" ); printf ( "\n" ); n_data = 0; for ( ; ; ) { bessel_j1_values ( &n_data, &x, &fx1 ); if ( n_data <= 0 ) { break; } fx2 = gsl_sf_bessel_J1 ( x ); printf ( " %12.6g %12.6g %12.6g\n", x, fx1, fx2 ); } return; } /******************************************************************************/ void gsl_sf_coupling_3j_test ( ) /******************************************************************************/ /* Purpose: gsl_sf_coupling_3j_test() tests gsl_sf_coupling_3j(). Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: John Burkardt */ { double fx; double fx2; double j1; double j2; double j3; double m1; double m2; double m3; int n_data; int twoj1; int twoj2; int twoj3; int twom1; int twom2; int twom3; printf ( "\n" ); printf ( "gsl_sf_coupling_3j_test():\n" ); printf ( " gsl_sf_coupling_3j() returns values of\n" ); printf ( " the Wigner 3J coefficient.\n" ); printf ( "\n" ); printf ( " J1 J2 J3 M1 M2 M3 THREE_J\n" ); printf ( "\n" ); n_data = 0; for ( ; ; ) { three_j_values ( &n_data, &j1, &j2, &j3, &m1, &m2, &m3, &fx ); if ( n_data == 0 ) { break; } twoj1 = round ( 2.0 * j1 ); twoj2 = round ( 2.0 * j2 ); twoj3 = round ( 2.0 * j3 ); twom1 = round ( 2.0 * m1 ); twom2 = round ( 2.0 * m2 ); twom3 = round ( 2.0 * m3 ); fx2 = gsl_sf_coupling_3j ( twoj1, twoj2, twoj3, twom1, twom2, twom3 ); printf ( " %6g %6g %6g %6g %6g %6g %24.16g\n", j1, j2, j3, m1, m2, m3, fx ); printf ( " %24.16g\n", fx2 ); } return; } /******************************************************************************/ void gsl_sf_coupling_6j_test ( ) /******************************************************************************/ /* Purpose: gsl_sf_coupling_6j_test() tests gsl_sf_coupling_6j(). Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: John Burkardt */ { double fx; double fx2; double j1; double j2; double j3; double j4; double j5; double j6; int n_data; int twoj1; int twoj2; int twoj3; int twoj4; int twoj5; int twoj6; printf ( "\n" ); printf ( "gsl_sf_coupling_6j_test():\n" ); printf ( " gsl_sf_coupling_6j() returns values of\n" ); printf ( " the Wigner 6J coefficient.\n" ); printf ( "\n" ); printf ( " J1 J2 J3 J4 J5 J6 SIX_J\n" ); printf ( "\n" ); n_data = 0; for ( ; ; ) { six_j_values ( &n_data, &j1, &j2, &j3, &j4, &j5, &j6, &fx ); if ( n_data == 0 ) { break; } twoj1 = round ( 2.0 * j1 ); twoj2 = round ( 2.0 * j2 ); twoj3 = round ( 2.0 * j3 ); twoj4 = round ( 2.0 * j4 ); twoj5 = round ( 2.0 * j5 ); twoj6 = round ( 2.0 * j6 ); fx2 = gsl_sf_coupling_6j ( twoj1, twoj2, twoj3, twoj4, twoj5, twoj6 ); printf ( " %6g %6g %6g %6g %6g %6g %24.16g\n", j1, j2, j3, j4, j5, j6, fx ); printf ( " %24.16g\n", fx2 ); } return; } /******************************************************************************/ void gsl_sf_coupling_9j_test ( ) /******************************************************************************/ /* Purpose: gsl_sf_coupling_9j_test() tests gsl_sf_coupling_9j(). Licensing: This code is distributed under the MIT license. Modified: 15 February 2022 Author: John Burkardt */ { double fx1; double fx2; double j1; double j2; double j3; double j4; double j5; double j6; double j7; double j8; double j9; int n_data; int twoj1; int twoj2; int twoj3; int twoj4; int twoj5; int twoj6; int twoj7; int twoj8; int twoj9; printf ( "\n" ); printf ( "gsl_sf_coupling_9j_test():\n" ); printf ( " gsl_sf_coupling_9j() returns values of\n" ); printf ( " the Wigner 9J coefficient.\n" ); printf ( "\n" ); printf ( " SOMETHING IS SERIOUSLY WRONG HERE.\n" ); printf ( " ALL RESULTS ARE 0.\n" ); printf ( " JVB, 15 February 2022\n" ); printf ( "\n" ); printf ( " J1 J2 J3 J4 J5 J6 J7 J8 J9 NINE_J\n" ); printf ( "\n" ); n_data = 0; for ( ; ; ) { nine_j_values ( &n_data, &j1, &j2, &j3, &j4, &j5, &j6, &j7, &j8, &j9, &fx1 ); if ( n_data == 0 ) { break; } twoj1 = ( int ) ( 2.0 * j1 ); twoj2 = ( int ) ( 2.0 * j2 ); twoj3 = ( int ) ( 2.0 * j3 ); twoj4 = ( int ) ( 2.0 * j4 ); twoj5 = ( int ) ( 2.0 * j5 ); twoj6 = ( int ) ( 2.0 * j6 ); twoj7 = ( int ) ( 2.0 * j7 ); twoj8 = ( int ) ( 2.0 * j8 ); twoj9 = ( int ) ( 2.0 * j9 ); fx2 = gsl_sf_coupling_9j ( twoj1, twoj2, twoj3, twoj4, twoj5, twoj6, twoj7, twoj8, twoj9 ); printf ( " %6g %6g %6g %6g %6g %6g %6g %6g %6g %24.16g\n", j1, j2, j3, j4, j5, j6, j7, j8, j9, fx1 ); printf ( " %6d %6d %6d %6d %6d %6d %6d %6d %6d %24.16g\n", twoj1, twoj2, twoj3, twoj4, twoj5, twoj6, twoj7, twoj8, twoj9, fx2 ); } return; } /******************************************************************************/ void gsl_sf_dawson_test ( ) /******************************************************************************/ /* Purpose: gsl_sf_dawson_test() tests gsl_sf_dawson(). Licensing: This code is distributed under the MIT license. Modified: 14 February 2022 Author: John Burkardt */ { double fx1; double fx2; int n_data; double x; printf ( "\n" ); printf ( "gsl_sf_dawson_test():\n" ); printf ( " gsl_sf_dawson() evaluates Dawson's integral.\n" ); printf ( "\n" ); printf ( " X Exact Computed\n" ); printf ( "\n" ); n_data = 0; for ( ; ; ) { dawson_values ( &n_data, &x, &fx1 ); if ( n_data <= 0 ) { break; } fx2 = gsl_sf_dawson ( x ); printf ( " %12.6g %12.6g %12.6g\n", x, fx1, fx2 ); } return; } /******************************************************************************/ void gsl_sf_hyperg_2F1_test ( ) /******************************************************************************/ /* Purpose: gsl_sf_hyperg_2F1_test() tests gsl_sf_hyperg_2F1(). Licensing: This code is distributed under the MIT license. Modified: 26 December 2023 Author: John Burkardt */ { double a; double b; double c; double fx1; double fx2; int n_data; double x; printf ( "\n" ); printf ( "gsl_sf_hyperg_2F1_test():\n" ); printf ( " Test gsl_sf_hyperg_2F1(), which evaluates the\n" ); printf ( " hypergeometric function 2F1 for real parameters and arguments.\n" ); printf ( "\n" ); printf ( " A B C X Hyper_2F1(A,B,C;X)\n" ); printf ( "\n" ); n_data = 0; for ( ; ; ) { hyper_2f1_values ( &n_data, &a, &b, &c, &x, &fx1 ); if ( n_data == 0 ) { break; } fx2 = gsl_sf_hyperg_2F1 ( a, b, c, x ); printf ( " %8f %8f %8f %8f %24.16g %24.16g\n", a, b, c, x, fx1, fx2 ); } return; } /******************************************************************************/ void gsl_spmatrix_test ( ) /******************************************************************************/ /* Purpose: gsl_spmatrix_test() tests gsl_spmatrix(). Licensing: This code is distributed under the MIT license. Modified: 29 September 2024 Author: Original C version by GSL reference manual. This version by John Burkardt. */ { /* Allocate space for the triplet sparse matrix format. */ gsl_spmatrix *A = gsl_spmatrix_alloc(5, 4); gsl_spmatrix *B; gsl_spmatrix *C; size_t i; size_t j; printf ( "\n" ); printf ( "gsl_spmatrix_test():\n" ); printf ( " Demonstrate gsl sparse matrix storage.\n" ); /* Build the sparse matrix. */ gsl_spmatrix_set ( A, 0, 2, 3.1 ); gsl_spmatrix_set ( A, 0, 3, 4.6 ); gsl_spmatrix_set ( A, 1, 0, 1.0 ); gsl_spmatrix_set ( A, 1, 2, 7.2 ); gsl_spmatrix_set ( A, 3, 0, 2.1 ); gsl_spmatrix_set ( A, 3, 1, 2.9 ); gsl_spmatrix_set ( A, 3, 3, 8.5 ); gsl_spmatrix_set ( A, 4, 0, 4.1 ); printf ( "\n" ); printf ( " Print all matrix elements:\n" ); for ( i = 0; i < 5; i++ ) { for ( j = 0; j < 4; j++ ) { printf ( "A(%zu,%zu) = %g\n", i, j, gsl_spmatrix_get(A, i, j) ); } } /* Print elements in triplet format */ printf ( "\n" ); printf ( " matrix in triplet (COO) format (i,j,Aij):\n" ); gsl_spmatrix_fprintf ( stdout, A, "%.1f" ); /* Convert to compressed column format. */ B = gsl_spmatrix_ccs ( A ); printf ( "\n" ); printf ( "matrix in compressed column sparse (CCS) format:\n" ); printf ( "i = [ " ); for ( i = 0; i < B->nz; ++i ) { printf ( "%d, ", B->i[i] ); } printf ( "]\n" ); printf ( "p = [ " ); for ( i = 0; i < B->size2 + 1; ++i ) { printf ( "%d, ", B->p[i] ); } printf ( "]\n" ); printf ( "d = [ " ); for ( i = 0; i < B->nz; ++i ) { printf ( "%g, ", B->data[i] ); } printf ( "]\n" ); /* Convert to compressed row format. */ C = gsl_spmatrix_crs ( A ); printf ( "\n" ); printf ( " matrix in compressed row sparse (CRS) format:\n" ); printf ( "i = [ " ); for ( i = 0; i < C->nz; ++i ) { printf ( "%d, ", C->i[i] ); } printf ( "]\n" ); printf ( "p = [ " ); for ( i = 0; i < C->size1 + 1; ++i ) { printf ( "%d, ", C->p[i] ); } printf ( "]\n" ); printf ( "d = [ " ); for ( i = 0; i < C->nz; ++i ) { printf ( "%g, ", C->data[i] ); } printf ( "]\n" ); /* Free memory. */ gsl_spmatrix_free ( A ); gsl_spmatrix_free ( B ); gsl_spmatrix_free ( C ); return; } /******************************************************************************/ void hyper_2f1_values ( int *n_data, double *a, double *b, double *c, double *x, double *fx ) /******************************************************************************/ /* Purpose: hyper_2f1_values() returns some values of the hypergeometric function 2F1. Discussion: In Mathematica, the function can be evaluated by: fx = Hypergeometric2F1 [ a, b, c, x ] Licensing: This code is distributed under the MIT license. Modified: 09 September 2007 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Shanjie Zhang, Jianming Jin, Computation of Special Functions, Wiley, 1996, ISBN: 0-471-11963-6, LC: QA351.C45 Daniel Zwillinger, CRC Standard Mathematical Tables and Formulae, 30th Edition, CRC Press, 1996, pages 651-652. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. The routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *A, *B, *C: the parameters. double *X: the argument. double *FX: the value of the function. */ { # define N_MAX 24 static double a_vec[N_MAX] = { -2.5, -0.5, 0.5, 2.5, -2.5, -0.5, 0.5, 2.5, -2.5, -0.5, 0.5, 2.5, 3.3, 1.1, 1.1, 3.3, 3.3, 1.1, 1.1, 3.3, 3.3, 1.1, 1.1, 3.3 }; static double b_vec[N_MAX] = { 3.3, 1.1, 1.1, 3.3, 3.3, 1.1, 1.1, 3.3, 3.3, 1.1, 1.1, 3.3, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7 }; static double c_vec[N_MAX] = { 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, -5.5, -0.5, 0.5, 4.5, -5.5, -0.5, 0.5, 4.5, -5.5, -0.5, 0.5, 4.5 }; static double fx_vec[N_MAX] = { 0.72356129348997784913, 0.97911109345277961340, 1.0216578140088564160, 1.4051563200112126405, 0.46961431639821611095, 0.95296194977446325454, 1.0512814213947987916, 2.3999062904777858999, 0.29106095928414718320, 0.92536967910373175753, 1.0865504094806997287, 5.7381565526189046578, 15090.669748704606754, -104.31170067364349677, 21.175050707768812938, 4.1946915819031922850, 1.0170777974048815592E+10, -24708.635322489155868, 1372.2304548384989560, 58.092728706394652211, 5.8682087615124176162E+18, -4.4635010147295996680E+08, 5.3835057561295731310E+06, 20396.913776019659426 }; static double x_vec[N_MAX] = { 0.25, 0.25, 0.25, 0.25, 0.55, 0.55, 0.55, 0.55, 0.85, 0.85, 0.85, 0.85, 0.25, 0.25, 0.25, 0.25, 0.55, 0.55, 0.55, 0.55, 0.85, 0.85, 0.85, 0.85 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *a = 0.0; *b = 0.0; *c = 0.0; *x = 0.0; *fx = 0.0; } else { *a = a_vec[*n_data-1]; *b = b_vec[*n_data-1]; *c = c_vec[*n_data-1]; *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ void nine_j_values ( int *n_data, double *j1, double *j2, double *j3, double *j4, double *j5, double *j6, double *j7, double *j8, double *j9, double *fx ) /******************************************************************************/ /* Purpose: nine_j_values() returns some values of the Wigner 9J function. Licensing: This code is distributed under the MIT license. Modified: 09 February 2007 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. The routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *J1, *J2, *J3, *J4, *J5, *J6, *J7, *J8, *J9, the arguments of the function. double *FX, the value of the function. */ { # define N_MAX 9 static double fx_vec[N_MAX] = { 0.0004270039294528318, -0.001228915451058514, -0.0001944260688400887, 0.003338419923885592, -0.0007958936865080434, -0.004338208690251972, 0.05379143536399187, 0.006211299937499411, 0.03042903097250921 }; static double j1_vec[N_MAX] = { 1.0, 1.5, 2.0, 1.0, 1.5, 2.0, 0.5, 1.0, 1.5 }; static double j2_vec[N_MAX] = { 8.0, 8.0, 8.0, 3.0, 3.0, 3.0, 0.5, 0.5, 0.5 }; static double j3_vec[N_MAX] = { 7.0, 7.0, 7.0, 2.0, 2.0, 2.0, 1.0, 1.0, 1.0 }; static double j4_vec[N_MAX] = { 6.5, 6.5, 6.5, 4.0, 4.0, 4.0, 2.0, 2.0, 2.0 }; static double j5_vec[N_MAX] = { 7.5, 7.5, 7.5, 1.5, 1.5, 1.5, 1.0, 1.0, 1.0 }; static double j6_vec[N_MAX] = { 7.5, 7.5, 7.5, 3.0, 3.0, 3.0, 1.5, 1.5, 1.5 }; static double j7_vec[N_MAX] = { 6.0, 6.0, 6.0, 3.5, 3.5, 3.5, 1.5, 1.5, 1.5 }; static double j8_vec[N_MAX] = { 10.0, 10.0, 10.0, 2.0, 2.0, 2.0, 0.5, 0.5, 0.5 }; static double j9_vec[N_MAX] = { 6.0, 6.0, 6.0, 2.0, 2.0, 2.0, 1.5, 1.5, 1.5 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *j1 = 0.0; *j2 = 0.0; *j3 = 0.0; *j4 = 0.0; *j5 = 0.0; *j6 = 0.0; *j7 = 0.0; *j8 = 0.0; *j9 = 0.0; *fx = 0.0; } else { *j1 = j1_vec[*n_data-1]; *j2 = j2_vec[*n_data-1]; *j3 = j3_vec[*n_data-1]; *j4 = j4_vec[*n_data-1]; *j5 = j5_vec[*n_data-1]; *j6 = j6_vec[*n_data-1]; *j7 = j7_vec[*n_data-1]; *j8 = j8_vec[*n_data-1]; *j9 = j9_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ int ode_func ( double t, const double y[], double f[], void *params ) /******************************************************************************/ /* Purpose: ode_func() evaluates the right hand side of the ODE. Discussion: This function defines the right hand side of the van der Pol equation: u'' + mu u' (u^2-1) + u = 0 rewritten in first order form as: u' = v v' = - u + mu v ( 1 - u^2) Licensing: This code is distributed under the MIT license. Modified: 01 May 2021 Author: Original C version by GSL reference manual. This version by John Burkardt. Input: double T: the current time. const double Y[]: the current solution value. void *PARAMS: user-specified parameters. Output: int FUNC: returns GSL_SUCCESS indicating a successful computation. double F[]: the right hand side of the ODE at T, Y. */ { /* Declare t to avoid unused parameter warning. */ (void) (t); double mu = *(double *) params; f[0] = y[1]; f[1] = - y[0] - mu * y[1] * ( y[0] * y[0] - 1.0 ); return GSL_SUCCESS; } /******************************************************************************/ int ode_jac ( double t, const double y[], double *dfdy, double dfdt[], void *params ) /******************************************************************************/ /* Purpose: ode_jac() evaluates the jacobian of the ODE. Discussion: This function defines the jacobian of the van der Pol equation: f(u,v) = ( v ) ( -u+mu*v*(1-u^2 ) df = ( 0 ) dt ( 0 ) df = ( 0 1 ) dy ( -1-2*mu*v*u mu*(1-u^2) ) Licensing: This code is distributed under the MIT license. Modified: 01 May 2021 Author: Original C version by GSL reference manual. This version by John Burkardt. Input: double T: the current time. const double Y[]: the current solution value. void *PARAMS: user-specified parameters. Output: int JAC: returns GSL_SUCCESS indicating a successful computation. double *DFDY[]: the Jacobian matrix. double DFDT[]: the time derivative of F. */ { /* Declare t to avoid unused parameter warning. */ (void) (t); double mu = *(double *) params; /* Allocate a 2x2 matrix for dfdy. */ gsl_matrix_view dfdy_mat = gsl_matrix_view_array ( dfdy, 2, 2 ); /* Let "m" be a pointer to this matrix. */ gsl_matrix * m = &dfdy_mat.matrix; /* Set dFdY */ gsl_matrix_set ( m, 0, 0, 0.0 ); gsl_matrix_set ( m, 0, 1, 1.0 ); gsl_matrix_set ( m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0 ); gsl_matrix_set ( m, 1, 1, -mu*(y[0]*y[0] - 1.0) ); /* Set dFdT */ dfdt[0] = 0.0; dfdt[1] = 0.0; return GSL_SUCCESS; } /******************************************************************************/ int rosenbrock_f ( const gsl_vector *x, void *params, gsl_vector *f ) /******************************************************************************/ /* Purpose: rosenbrock_f() evaluates the Rosenbrock function. Licensing: This code is distributed under the MIT license. Modified: 05 August 2005. Author: Original C version by GSL reference manual. This version by John Burkardt. */ { double a = ( ( struct rparams *) params )-> a; double b = ( ( struct rparams *) params )-> b; const double x0 = gsl_vector_get ( x, 0 ); const double x1 = gsl_vector_get ( x, 1 ); const double y0 = a * ( 1.0 - x0 ); const double y1 = b * ( x1 - x0 * x0 ); gsl_vector_set ( f, 0, y0 ); gsl_vector_set ( f, 1, y1 ); return GSL_SUCCESS; } /******************************************************************************/ void six_j_values ( int *n_data, double *j1, double *j2, double *j3, double *j4, double *j5, double *j6, double *fx ) /******************************************************************************/ /* Purpose: six_j_values() returns some values of the Wigner 6J function. Discussion: In Mathematica, the function can be evaluated by: SixJSymbol[{j1,j2,j3},{j4,j5,j6}] Licensing: This code is distributed under the MIT license. Modified: 07 February 2007 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. The routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *J1, *J2, *J3, *J4, *J5, *J6, the arguments of the function. double *FX, the value of the function. */ { # define N_MAX 15 static double fx_vec[N_MAX] = { 0.03490905138373300, -0.03743025039659792, 0.01890866390959560, 0.007342448254928643, -0.02358935185081794, 0.01913476955215437, 0.001288017397724172, -0.01930018366290527, 0.01677305949382889, 0.005501147274850949, -0.02135439790896831, 0.003460364451435387, 0.02520950054795585, 0.01483990561221713, 0.002708577680633186 }; static double j1_vec[N_MAX] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0 }; static double j2_vec[N_MAX] = { 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0 }; static double j3_vec[N_MAX] = { 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0 }; static double j4_vec[N_MAX] = { 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5 }; static double j5_vec[N_MAX] = { 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5 }; static double j6_vec[N_MAX] = { 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *j1 = 0.0; *j2 = 0.0; *j3 = 0.0; *j4 = 0.0; *j5 = 0.0; *j6 = 0.0; *fx = 0.0; } else { *j1 = j1_vec[*n_data-1]; *j2 = j2_vec[*n_data-1]; *j3 = j3_vec[*n_data-1]; *j4 = j4_vec[*n_data-1]; *j5 = j5_vec[*n_data-1]; *j6 = j6_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ void three_j_values ( int *n_data, double *j1, double *j2, double *j3, double *m1, double *m2, double *m3, double *fx ) /******************************************************************************/ /* Purpose: three_j_values() returns some values of the Wigner 3J function. Discussion: In Mathematica, the function can be evaluated by: ThreeJSymbol[{j1,m1},{j2,m2},{j3,m3}] Licensing: This code is distributed under the MIT license. Modified: 07 February 2007 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA. The user sets N_DATA to 0 before the first call. Output: int *N_DATA. The routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *J1, *J2, *J3, *M1, *M2, *M3, the arguments of the function. double *FX, the value of the function. */ { # define N_MAX 8 static double fx_vec[N_MAX] = { 0.2788866755113585, -0.09534625892455923, -0.06741998624632421, 0.1533110351679666, -0.1564465546936860, 0.1099450412156551, -0.05536235693131719, 0.01799835451137786 }; static double j1_vec[N_MAX] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0 }; static double j2_vec[N_MAX] = { 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5 }; static double j3_vec[N_MAX] = { 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5 }; static double m1_vec[N_MAX] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 }; static double m2_vec[N_MAX] = { -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5 }; static double m3_vec[N_MAX] = { 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *j1 = 0.0; *j2 = 0.0; *j3 = 0.0; *m1 = 0.0; *m2 = 0.0; *m3 = 0.0; *fx = 0.0; } else { *j1 = j1_vec[*n_data-1]; *j2 = j2_vec[*n_data-1]; *j3 = j3_vec[*n_data-1]; *m1 = m1_vec[*n_data-1]; *m2 = m2_vec[*n_data-1]; *m3 = m3_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: timestamp() prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE }