Touuki的博客 > 程序开发 > C++数值积分 – 基于GSL的简单封装类

C++数值积分 – 基于GSL的简单封装类

作者Touuki发布于 分类目录 程序开发 566

5月28日追记　其实现在最新版本的Boost库数值积分已经比较完善了，处理常见的绝大多数积分都没有问题，因此C++的数值积分也可以考虑使用Boost，它的接口更加友好、使用也更方便，不需要像这篇文章一样再进行封装了。

　　C/C++中的数值积分在不考虑引用库Licence类型情况下，可以使用GNU Scientific Library (GSL)这个库来进行处理。如果只是需要计算一两种固定种类的积分，那么直接按照文档中的例子来写没什么问题。但是如果需要计算的积分种类很多的话，比如一个积分中还嵌套着其他好几个积分时，按照文档的写法写出来代码会变得非常丑陋、不可读以及难以维护。因为每算一个积分都需要构造一个函数double f (double x, void * params)，这会造成两个问题：一是函数定义和使用不在同一个位置，由于这个函数是专门为GSL积分所准备的，因此通常只用一次就不会再用到了，这样读代码和维护的时候非常麻烦；二是传参的时候需要自己将一个或多个参数装配成void * params的形式传递进函数，之后再在函数中将自定义参数从void * params中释放出来。

　　为了解决这两个问题，我决定利用C++的新特性对GSL的数值积分进行一下简单的封装。第一个问题非常适合用lambda表达式来解决，这样定义和使用都在同一个位置，更利于理解和维护。第二个问题其实是可变参数的问题，最适合我们这种情况的是可变模板参数，这样我们既不用手工装配也不用手工释放参数，lambda表达式传进来参数就可以直接使用。

　　所以我们的模板类可能是这个样子的：

template <typename... T> class Quadgk {
  private:
    double (*func)(double x, T... params);
    double (*gslFunc)(double x, void *params);
  public:
    Quadgk(double (*_func)(double x, T... params));
    double operator()(double a, double b, T... params);
};

　　不过这样我们遇到了一个难题，由于我们最终是调用GSL进行数值积分运算，因此传给GSL的是double (*gslFunc)(double x, void *params)这种形式的函数，而我们lambda表达式的形式是double (*func)(double x, T... params)那这就意味着要在gslFunc中调用func。但是怎么把T...通过void *传递进去，又该怎么调用这个含有可变模板参数的函数呢？

　　最终我们在这篇博客『Modern C++: Variadic template parameters and tuples』中找到了答案，至于原理是啥我是看不懂了，反正直接抄过来能用就对了……所以最终封装好的类长这样：

/*
 * Copyright (C) 2020 Touuki
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
 * USA.
 */

#include <gsl/gsl_integration.h>
#include <tuple>
#include <utility>

template <typename... T> class Quadgk {
  private:
    gsl_integration_workspace *w;
    double epsabs;
    double epsrel;
    double (*func)(double x, T... params);
    double (*gslFunc)(double x, void *params);

  public:
    Quadgk(double (*_func)(double x, T... params), size_t _limit = 1000,
           double _epsabs = 0., double _epsrel = 1e-8)
        : w(gsl_integration_workspace_alloc(_limit)), epsabs(_epsabs),
          epsrel(_epsrel), func(_func),
          gslFunc([](double x, void *params) -> double {
              void **p = reinterpret_cast<void **>(params);
              Quadgk *self = reinterpret_cast<Quadgk *>(*p);
              std::tuple<T...> tuple =
                  *reinterpret_cast<std::tuple<T...> *>(*(p + 1));
              return self->callFunc(x, tuple, std::index_sequence_for<T...>());
          }){};
    ~Quadgk() { gsl_integration_workspace_free(w); };
    void setEpsabs(double _epsabs) { epsabs = _epsabs; };
    void setEpsrel(double _epsrel) { epsrel = _epsrel; };
    template <std::size_t... Is>
    double callFunc(double x, const std::tuple<T...> &amp;tuple,
                    std::index_sequence<Is...>) const {
        return func(x, std::get<Is>(tuple)...);
    };
    double operator()(double a, double b, T... params) {
        std::tuple<T...> tuple(params...);
        void *p[] = {this, &amp;tuple};
        gsl_function f = {gslFunc, &amp;p};
        double result, abserr;
        if (isinf(a)) {
            if (isinf(b)) {
                gsl_integration_qagi(&amp;f, epsabs, epsrel, w->limit, w, &amp;result,
                                     &amp;abserr);
            } else {
                gsl_integration_qagil(&amp;f, b, epsabs, epsrel, w->limit, w,
                                      &amp;result, &amp;abserr);
            }
        } else {
            if (isinf(b)) {
                gsl_integration_qagiu(&amp;f, a, epsabs, epsrel, w->limit, w,
                                      &amp;result, &amp;abserr);
            } else {
                gsl_integration_qags(&amp;f, a, b, epsabs, epsrel, w->limit, w,
                                     &amp;result, &amp;abserr);
            }
        }
        return result;
    };
};

template <typename... T> class QuadSine {
  private:
    gsl_integration_qawo_table *wf;
    gsl_integration_workspace *workspace;
    gsl_integration_workspace *cycle_workspace;
    double epsabs;
    double epsrel;
    double (*func)(double x, T... params);
    double (*gslFunc)(double x, void *params);

  public:
    QuadSine(double (*_func)(double x, T... params),
             gsl_integration_qawo_enum _sine, size_t _n = 3,
             size_t _limit = 1000, double _epsabs = 1e-8, double _epsrel = 1e-8)
        : wf(gsl_integration_qawo_table_alloc(1., 1., _sine, _n)),
          workspace(gsl_integration_workspace_alloc(_limit)),
          cycle_workspace(gsl_integration_workspace_alloc(_limit)),
          epsabs(_epsabs), epsrel(_epsrel), func(_func),
          gslFunc([](double x, void *params) -> double {
              void **p = reinterpret_cast<void **>(params);
              QuadSine *self = reinterpret_cast<QuadSine *>(*p);
              std::tuple<T...> tuple =
                  *reinterpret_cast<std::tuple<T...> *>(*(p + 1));
              return self->callFunc(x, tuple, std::index_sequence_for<T...>());
          }){};
    ~QuadSine() {
        gsl_integration_qawo_table_free(wf);
        gsl_integration_workspace_free(workspace);
        gsl_integration_workspace_free(cycle_workspace);
    };
    void setEpsabs(double _epsabs) { epsabs = _epsabs; };
    void setEpsrel(double _epsrel) { epsrel = _epsrel; };
    template <std::size_t... Is>
    double callFunc(double x, const std::tuple<T...> &amp;tuple,
                    std::index_sequence<Is...>) const {
        return func(x, std::get<Is>(tuple)...);
    };
    double operator()(double a, double b, double omega, T... params) {
        std::tuple<T...> tuple(params...);
        void *p[] = {this, &amp;tuple};
        gsl_function f = {gslFunc, &amp;p};
        double result, abserr;
        if (isinf(a)) {
            throw std::invalid_argument("a can't be INFINITY.");
        } else {
            if (isinf(b)) {
                if (wf->omega != omega) {
                    gsl_integration_qawo_table_set(wf, omega, 1., wf->sine);
                }
                gsl_integration_qawf(&amp;f, a, epsabs, workspace->limit, workspace,
                                     cycle_workspace, wf, &amp;result, &amp;abserr);
            } else {
                double L = b - a;
                if (wf->omega != omega) {
                    gsl_integration_qawo_table_set(wf, omega, L, wf->sine);
                } else if (wf->L != L) {
                    gsl_integration_qawo_table_set_length(wf, L);
                }
                gsl_integration_qawo(&amp;f, a, epsabs, epsrel, workspace->limit,
                                     workspace, wf, &amp;result, &amp;abserr);
            }
        }
        return result;
    };
};

　　这样使用起来就非常简单方便了，这里给一个嵌套积分的例子：

#include <gsl/gsl_specfunc.h>
#include <iostream>

int main(int argc, char const *argv[]) {
    Quadgk<double> quad([](double chi, double eta) {
        static Quadgk<double> intK(
            [](double x, double nu) { return gsl_sf_bessel_Knu(nu, x); });
        double u = 2. * chi / (eta - 2. * chi);
        double u_p = 2. * u / (3. * eta);
        return (1. / (1. + u) + 1. + u) * gsl_sf_bessel_Knu(2. / 3., u_p) -
               intK(u_p, INFINITY, 1. / 3.);
    });
    double eta = 10;
    std::cout << quad(0., eta / 2., eta) << std::endl;
    return 0;
};

发表评论

电子邮件地址不会被公开。 必填项已用*标注