1

我是一个非常糟糕的程序员，所以要小心我犯了一个非常愚蠢的错误的可能性很大。我正在编写 C 并且使用gsl_multiroot_fsolver_hybrid错误，任何帮助都会非常有用！

我的代码包含以下两个功能：

/* ---  --- */
int quasilognormal_f (const gsl_vector * x, void *params,
               gsl_vector * f)
 {

   const double x0 = gsl_vector_get (x, 0);
   const double x1 = gsl_vector_get (x, 1);
   const double x2 = gsl_vector_get (x, 2);


   qparams->N =x0;
   qparams->A =x1;
   qparams->S =x2;
                printf("INNER POINT 1 -- Attention!\n"); 
                printf("N,A,S = %e %e %e \n",
                    x0, x1, x2);
   const double y0 = norm_quasilognormal(qparams) - 1.0;
   //const double y0 = 0;
                    printf("INNER POINT 2\n"); 
   const double y1 = mean_quasilognormal(qparams) - 0.0;
   //const double y1 =0;
                     printf("INNER POINT 3\n"); 
   const double y2 = var_quasilognormal(qparams) - kappa_var(gcosmo-> theta);
                     printf("INNER POINT 4\n");

   gsl_vector_set (f, 0, y0);
   gsl_vector_set (f, 1, y1);
   gsl_vector_set (f, 2, y2);

   return GSL_SUCCESS;
 }

int COMPUTE_QUASILOGNORMAL_PARAMETERS()
 {

 gsl_vector * root;
 const gsl_multiroot_fsolver_type *T;
 gsl_multiroot_fsolver *s;

   int status;
   size_t i, iter = 0;

   const size_t n = 3;

   gsl_multiroot_function F;

F.f = &quasilognormal_f;
F.params = NULL;
F.n = 3;

   double x_init[3] = {1.0, 0.0,gcosmo->ausgelagertes_lognormal_sigma};
   gsl_vector *x = gsl_vector_alloc (n);


   gsl_vector_set (x, 0, x_init[0]);
   gsl_vector_set (x, 1, x_init[1]);
   gsl_vector_set (x, 2, x_init[2]);


   T = gsl_multiroot_fsolver_hybrid;
   s = gsl_multiroot_fsolver_alloc (T, 3);

   gsl_multiroot_fsolver_set (s, &F, x);
   printf("THE SOLVER WAS SET!\n");
   do
     {
       iter++;
            printf("ABOUT TO ITERATE SOLVER\n");
                        printf("iter = %i\n",
                    iter);


                printf("x_init = %e %e %e \n",
                    x_init[0], x_init[1], x_init[2]);

                    printf("N,A,S = %e %e %e \n",
                    root[0], root[1], root[2]);
       status = gsl_multiroot_fsolver_iterate (s);

                printf("THE SOLVER WAS ITERATED\n");
                        printf("iter = %i\n",
                    iter);

                printf("x_init = %e %e %e \n",
                    x_init[0], x_init[1], x_init[2]);

                    printf("N,A,S = %e %e %e \n",
                    root[0], root[1], root[2]);
       if (status)   /* check if solver is stuck */
         break;

       status =
         gsl_multiroot_test_residual (s->f, 1e-7);
     }
   while (status == GSL_CONTINUE && iter < 1000);

   printf ("status = %s\n", gsl_strerror (status));

          root = gsl_multiroot_fsolver_root(s);

          printf("x = % .3e \n",
                root[0], root[1], root[2]);


   gsl_multiroot_fsolver_free (s);
   gsl_vector_free (x);
   return 0;
 }     
/* ---  --- */

函数int COMPUTE_QUASILOGNORMAL_PARAMETERS应该求解在int quasilognormal_f中定义的三个方程组。这会出错，因为一旦第二次 (iter ==1) 迭代开始，参数 N、A 和 S在int quasilognormal_f中就是NAN。这是我从程序输出中得出的结论：

INNER POINT 1 -- Attention!
N,A,S = 1.000000e+00 0.000000e+00 2.071317e-01 
INNER POINT 2
INNER POINT 3
INNER POINT 4
INNER POINT 1 -- Attention!
N,A,S = 1.000000e+00 0.000000e+00 2.071317e-01 
INNER POINT 2
INNER POINT 3
INNER POINT 4
INNER POINT 1 -- Attention!
N,A,S = 1.000000e+00 1.490116e-08 2.071317e-01 
INNER POINT 2
INNER POINT 3
INNER POINT 4
INNER POINT 1 -- Attention!
N,A,S = 1.000000e+00 0.000000e+00 2.071318e-01 
INNER POINT 2
INNER POINT 3
INNER POINT 4
THE SOLVER WAS SET!
ABOUT TO ITERATE SOLVER
iter = 1
x_init = 1.000000e+00 0.000000e+00 2.071317e-01 
N,A,S = 1.000000e+00 0.000000e+00 2.071317e-01 
INNER POINT 1 -- Attention!
N,A,S = nan nan nan

我不知道该怎么做，因为求解器对我来说是一个黑匣子，我尝试按照 GSL 手册进行操作。

4

0 回答 0