使用带有参数的 Boost C++ 计算 2D 数值双积分

计算科学 C++ 一体化 正交 图书馆
2021-12-14 06:59:47

我正在尝试计算双 Richardson 和 Wolf 积分,用于在 C++ 中使用 Boost 进行镜头聚焦(使用 Gauss Kronrod 方法)。

作为起点,我使用了这个问题中提供的示例。我的问题是我必须使用可变参数执行几次积分(在图像空间中循环),这些参数主要取决于图像平面中点的位置,最终取决于瞳孔的光学特性。

为了在每次迭代时轻松更新参数(我不想通过全局变量传递)并简化被积函数的编写,我使用了一个Integrand类,其中一些公共成员函数对应于要集成的函数和要更新的参数的私有属性。

#ifndef INTEGRAND_H
#define INTEGRAND_H
#include <cmath>
#include <complex>


using namespace boost::math::quadrature;
using complex = std::complex<double>;

    class integrand
    {
    public:
        integrand();
        void set_k(double k);
        void set_rp(double rp);
        void set_thetap(double thetap);
        void set_phip(double phip);
        complex ex(double theta,double phi);
        complex ey(double theta,double phi);
        complex ez(double theta,double phi);
    private:
        double _k;
        double _rp;
        double _thetap;
        double _phip;
    
    };
    
    #endif // INTEGRAND_H

作为一个例子,这里是我如何实现允许计算沿 x 轴的电场分量的成员函数:

complex integrand::ex(double theta,double phi){
    complex cos_e = (cos(theta)*cos(_thetap)+sin(theta)*sin(_thetap)*cos(phi-_phip),0);
    complex part1 = (sqrt(cos(theta))*sin(theta)*(cos(theta)+1-cos(theta)*sin(phi)*sin(phi)),0);
    complex part2 = (cos(_k*_rp*cos_e),sin(_k*_rp.*cos_e));
    double pupil = 1;
    return part1*part2*pupil;
}

main在更新参数并遵循上述问题中提出的方法后,我在计算循环中(在我的)中调用此成员函数

///compute I0
        auto f0 = [](double theta){
            auto g0 = [&](double phi){
                return I.ex(theta, phi);
            };
            return gauss_kronrod<complex,61>::integrate(g0,0,2*M_PI, 5);
        };
        complex I0 = gauss_kronrod<double,61>::integrate(f0,0,a, 5, 1e-9, &error);

在那里,在编译期间,我得到一个I is not captured错误,然后是the lambda has no capture defaultand Integrand I declared here(指向我实例化的行I)。我知道我定义或调用 integrand::ex 函数的方式有问题,但我不知道如何解决这个问题。

这种方法有什么问题?我应该完全不同吗?

1个回答

编译器告诉您他无法捕获I. 这只是一个 C++ 问题。你可以看看什么是 lambda 表达式。基本上,它只是函子的简写。

例如,考虑以下情况:

double phi = 2.0;
auto capture_by_value = [] (double x) {return x + phi;};
std::cout << capture_by_value(2.0)<<"\n";

错误消息将与您之前收到的消息基本相同。此外,它是不言自明的:您需要捕获变量phi

//Better yet: capture phi by reference:
auto capture_by_reference = [&phi] (double x)->double {return x + phi;}; 
std::cout << capture_by_reference(2.0)<<"\n"; 

请注意,通过引用捕获,您可以修改捕获的值。如果您通过复制捕获,即简单地写入,则情况并非如此[phi]