Оболочка C++ для алгоритма поиска корня GSL с производной

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

Некоторое время назад я нашел очень удобный шаблон (оболочка GSL C++) который отлично работает для одной функции, например интегрировать, и я активно использую его.

Теперь мне интересно, можно ли расширить этот подход, чтобы предоставить две функции для GSL, а именно саму функцию и ее производную.

Изменить: Решение

template <typename F, typename G>
class gsl_root_deriv : public gsl_function_fdf
{
private:
    const F&    _f;
    const G&    _df;

    static double invoke_f(double x, void* params){
        return static_cast<gsl_root_deriv*>(params) -> _f(x);
    }

    static double invoke_df(double x, void* params){
        return static_cast<gsl_root_deriv*>(params) -> _df(x);
    }

    static void     invoke_fdf(double x, void* params, double* f, double* df){
        (*f)    = static_cast<gsl_root_deriv*>(params)  ->  _f(x);
        (*df)   = static_cast<gsl_root_deriv*>(params)  ->  _df(x);
    }

public:
    gsl_root_deriv(const F& f_init, const G& df_init)
    : _f(f_init), _df(df_init)
    {
        f               = &gsl_root_deriv::invoke_f;
        df          = &gsl_root_deriv::invoke_df;
        fdf         = &gsl_root_deriv::invoke_fdf;
        params  = this;
    }
};

И минимальный пример, напоминающий пример из GSL:

#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <memory>

#include "gsl_root_deriv.h"

int main()
{
auto f = [](double x) -> double { return 0.25 * x*x - 1.0; };
auto df = [](double x) -> double { return 0.5 * x; };

gsl_root_deriv<decltype(f),decltype(df)> Fp(f,df);

int status(0), iter(0), max_iter(100);

const gsl_root_fdfsolver_type* T( gsl_root_fdfsolver_newton);

std::unique_ptr<gsl_root_fdfsolver,void(*)(gsl_root_fdfsolver*)>
    s(gsl_root_fdfsolver_alloc(T),gsl_root_fdfsolver_free);

double x_0(0.0), x(5.0);

gsl_root_fdfsolver_set( s.get(), &Fp, x );

do
{
    iter++;
    std::cout << "Iteration " << iter << std::endl;
    status = gsl_root_fdfsolver_iterate( s.get() );
    x_0 = x;
    x = gsl_root_fdfsolver_root( s.get() );
    status = gsl_root_test_delta( x, x_0, 0.0, 1.0e-3 );
} while( status == GSL_CONTINUE && iter < max_iter );

std::cout << "Converged to root " << x << std::endl;

return 0;
}

С уважением,
-- Клаус


c++ wrapper gsl
person Community    schedule 06.07.2014    source источник
comment
Если вам так понравилось, вы можете проголосовать за этот ответ: p   -  person Vivian Miranda    schedule 07.07.2014
comment
Нужен второй параметр шаблона. Я забыл об этом!   -  person Vivian Miranda    schedule 07.07.2014
comment
Еще кое-что. Если у вас есть более эффективный способ вычисления fdf (f и df одновременно), вы можете добавить в конструктор третий аргумент. Здесь я просто показал базовую реализацию.   -  person Vivian Miranda    schedule 07.07.2014
comment
Честно говоря, я не понял вашего комментария. Кроме того, я должен признаться, что не понимаю, почему нужен второй аргумент шаблона. Я имею в виду, что и f, и df (из структуры gsl_function_fdf) являются указателями на функции, возвращающие двойное значение, поэтому они одного и того же «типа» или нет? Так почему же одного имени типа F недостаточно? Не стесняйтесь отвечать: «Мне пора почитать книгу» ;) С уважением --   -  person    schedule 07.07.2014
comment
GSL запрашивает третью функцию fdf, если у пользователя есть оптимизированный способ одновременного вычисления f и df. Если у вас есть лучшая реализация в этом случае, возможно, стоит добавить третью функцию в конструктор.   -  person Vivian Miranda    schedule 07.07.2014
comment
Это хороший вопрос C++, чтобы спросить @stackoverflow, зачем нужен второй параметр шаблона. Я тоже запутался.   -  person Vivian Miranda    schedule 07.07.2014

Ответы (1)


arrow_upward
1
arrow_downward

Вам нужно будет сделать некоторые модификации

Структура Gsl fdf следующая

struct gsl_function_fdf_struct 
{
  double (* f) (double x, void * params);
  double (* df) (double x, void * params);
  void (* fdf) (double x, void * params, double * f, double * df);
  void * params;
};

typedef struct gsl_function_fdf_struct gsl_function_fdf ;

Если вы поняли, что на самом деле делает обертка, то вы увидите, что обобщение довольно простое.

class gsl_function_fdf_pp : public gsl_function_fdf
{
   public:
   gsl_function_fdf_pp
   (
     std::function<double(double)> const& kf, 
     std::function<double(double)> const& kdf
   ) : _f(kf), _df(kdf)
   {
      f   = &gsl_function_fdf_pp::invoke;
      df  = &gsl_function_fdf_pp::invoke2;
      fdf = &gsl_function_fdf_pp::invoke3;
      params=this;
   }     
   private:
   std::function<double(double)> _f;
   std::function<double(double)> _df;

   static double invoke(double x, void *params) 
   {
     return static_cast<gsl_function_fdf_pp*>(params)->_f(x);
   }

   static double invoke2(double x, void *params) 
   {
     return static_cast<gsl_function_fdf_pp*>(params)->_df(x);
   }

   static void invoke3(double x, void * params, double* f, double* df)
   {
     (*f)  = static_cast<gsl_function_fdf_pp*>(params)->_f(x);
     (*df) = static_cast<gsl_function_fdf_pp*>(params)->_df(x);
   }
};

Версия шаблона.

template< typename F, typename U >  class gsl_function_fdf_pp : public gsl_function_fdf 
{
  public:
  gsl_function_fdf_pp(const F& kf, const U& kdf) : _f(kf), _df(kdf)
  {
    f   = &gsl_function_fdf_pp::invoke;
    df  = &gsl_function_fdf_pp::invoke2;
    fdf = &gsl_function_fdf_pp::invoke3;
    params=this;
  }
  private:
  const F& _f;
  const U& _df;

  static double invoke(double x, void *params) 
  {
    return static_cast<gsl_function_fdf_pp*>(params)->_f(x);
  }

  static double invoke2(double x, void *params) 
  {
    return static_cast<gsl_function_fdf_pp*>(params)->_df(x);
  }

  static void invoke3(double x, void * params, double* f, double* df)
  {
    (*f)  = static_cast<gsl_function_fdf_pp*>(params)->_f(x);
    (*df) = static_cast<gsl_function_fdf_pp*>(params)->_df(x);
  }
};

EDIT2: после исправления 2 небольших проблем этот пример работает

int main()
{
  auto f  = [](double x)->double{ return x*x; };
  auto df = [](double x)->double{ return 2.0 * x; };

  gsl_function_fdf_pp<decltype(f),decltype(df)> Fp(f,df);

  return 0;
}
person Vivian Miranda    schedule 07.07.2014
comment
Уважаемый Vinicius, Большое спасибо за ваш ответ. Я копался в каталоге gsl и успешно нашел структуру gsl_function_fdf в math.h. Я полагаю, что до сих пор я понял идею класса (шаблона): вы получаете наследование от структуры gsl_function_fdf и модифицируете ее в соответствии с интерфейсом C++. Из-за наследования объект gsl_function_fdf_pp is является объектом gsl_function_fdf, и gsl может его обрабатывать. Единственное, что я не совсем понимаю: почему вы отбрасываете указатель params? Не могли бы вы уточнить это? С уважением, Клаус - person ; 07.07.2014
comment
Потому что статические функции-члены не могут указывать на нестатические указатели-члены (например, на указатель this!). Поэтому мы используем указатель params void для ссылки на текущий мгновенный класс! - person Vivian Miranda; 07.07.2014
comment
Спасибо. Но шаблон не работает, смотрите минимальный пример в редактировании моего первого поста. С уважением - person ; 07.07.2014
comment
В имени конструктора была опечатка, и я забыл добавить второй параметр шаблона. - person Vivian Miranda; 07.07.2014
comment
Хорошо, теперь я думаю, что вы добавили слишком много запятой. ;) С уважением - person ; 07.07.2014