c++包装器的GSL根查找算法与导数
C++ wrapper for GSL root finding algorithm with derivative
所以,当我很高兴在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;
}
亲切的问候,
——克劳斯
你需要做一些修改
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:在修复了两个小问题之后,这个例子可以运行了
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;
}
相关文章:
- 算法问题:查找从堆栈中弹出的所有序列
- 这种用于查找连续子数组中最大和的递归算法有什么优势吗?
- 用于查找网格中最短路径的算法
- 查找最短路径算法
- C++ 查找算法:如何找到元素的最后一次出现?
- 用于在链表中查找节点的算法
- 用于查找ABS的最小值的算法(A [I] A [J] -K)
- 大 O 表示用于查找素数的算法
- 算法在容器中查找具有给定值的元素之一的成员
- 联合查找算法(查找周期),不使用 malloc
- 优化查找特定斐波那契数的算法
- 我如何在C 中制作算法,以在不重复的情况下查找集合的变化(即n元素,选择k)
- 使用STL算法查找集合中的所有匹配项
- 使用算法查找自定义数据向量中的最大值和最小值
- 使用 STL 算法查找集合中的前两个不相邻元素
- 使用Ford-Fulkerson算法查找边
- 使用 Kruskal 算法查找最小生成树的错误
- C++用指针算法查找字符串长度
- 此算法查找所有路径总和的时间复杂度是多少?
- 算法:查找给定范围内的数字计数