Блог

Я хочу использовать GSL для интеграцииhttp://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html Однако я не нахожу удобного способа, как интегрированная функция

(функция f в примере http://www.gnu.org/software/gsl/manual/html_node/Numerical-integration-examples.html )

может сообщать об ошибке интегратору. Я хочу интегрировать функцию, которая сама является результатом интеграции, которая может завершиться неудачей. Это мой пример программы

 #include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

double f (double x, void * params) {
    GSL_ERROR("test error",GSL_FAILURE);
    return 0.0;
}



int main (void)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

double result, error;

gsl_function F;
F.function = amp;f;

gsl_set_error_handler_off();
int status = gsl_integration_qags (amp;F, 0, 1, 0, 1e-7, 1000,
                     w, amp;result, amp;error); 

printf ("status          = %dn", status);
status  = GSL_FAILURE;
printf ("status          = %dn", status);


gsl_integration_workspace_free (w);

return 0;
}
  

в результате на выходе
статус = 0
статус = -1

Я думаю, интегратору лучше остановиться и вернуть мой код ошибки. Как я могу этого добиться?

Большое вам спасибо за вашу помощь!!!

2011-04-27: Я также попробовал этот вариант, после того как Брайан Гоф сказал мне,

 #include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

double f (double x, void * params) {
    GSL_ERROR("test error",GSL_FAILURE);
    return GSL_NAN;
}



int main (void)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

double result, error;

gsl_function F;
F.function = amp;f;

gsl_set_error_handler_off();
int status = gsl_integration_qags (amp;F, 0, 1, 0, 1e-7, 1000,
                     w, amp;result, amp;error); 

printf ("status          = %dn", status);
status  = GSL_FAILURE;
printf ("status          = %dn", status);


gsl_integration_workspace_free (w);

return 0;
}
  

это тоже не помогло. Теперь я заполню отчет об ошибке.

Благодаря Xuebin Wu из списка рассылки GSL проблема решена:

Привет,

GSL_ERROR сам по себе является макросом, он выглядит как

       gsl_error (reason, __FILE__, __LINE__, gsl_errno);
      return gsl_errno;
  

Функция уже возвращается до того, как вы вернете NAN, потому что GSL_ERROR
был вызван. Выключая обработчик, просто позвольте первой строке выполнить
ничего. Обработчик ошибок по умолчанию прерывает программу после печати
сообщение об ошибке.

Я не думаю, что это ошибка. Возможно, вы можете написать свой собственный обработчик ошибок для решения вашей проблемы. Например, вы можете использовать «goto» для перехода из gsl_integration_qags или задать некоторую глобальную переменную, чтобы указать, что результат интеграции неверен.

PS: Я считаю, что этот макрос — то, что вам нужно,

Макрос: GSL_ERROR_VAL (причина, gsl_errno, значение) Этот макрос такой же, как GSL_ERROR, но возвращает определенное пользователем значение value вместо кода ошибки. Его можно использовать для математических функций, которые возвращают значение с плавающей запятой.

В следующем примере показано, как вернуть NaN при математической особенности, используя макрос GSL_ERROR_VAL,

  if (x == 0)
   {
     GSL_ERROR_VAL("argument lies on singularity",
                   GSL_ERANGE, GSL_NAN);
   }
  

Итак, я скорректировал код в соответствии с

 #include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>

double f (double x, void * params) {
//  return GSL_NAN;
    GSL_ERROR_VAL ("argument lies on singularity", GSL_ERANGE, GSL_NAN);
}



int main (void)
{
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

double result, error;

gsl_function F;
F.function = amp;f;

gsl_set_error_handler_off();
int status = gsl_integration_qags (amp;F, 0, 1, 0, 1e-7, 1000,
                     w, amp;result, amp;error); 

printf ("status          = %dn", status);
status  = GSL_FAILURE;
printf ("status          = %dn", status);


gsl_integration_workspace_free (w);

return 0;
}
  

и все работает так, как ожидалось…

Немного по-хакерски, но я бы, вероятно, попросил вашу функцию сохранить какой-нибудь флаг. Когда он обнаруживает ошибку, он устанавливает флаг и возвращает ноль для всех последующих вычислений. Затем, после того как вы интегрировали его, вы можете проверить этот флаг, чтобы увидеть, является ли результат действительным.

Как насчет того, чтобы написать оболочку для функции, которая возвращает указатель на структуру, содержащую результаты функции и состояние ошибки? Или, если вы используете c , эта инкапсуляция может быть выполнена с использованием объектов ….