Разделение функций C++ с несколькими переменными для квадратурного интегрирования

Я хочу численно интегрировать с boost::math::quadrature::trapezoidal(g, a, b, 1e-6); Здесь я интегрирую функцию g(x). Проблема в том, что я должен выполнить двойной интеграл. Кроме того, у меня есть 4 переменные в функции, которую я хочу интегрировать. 2 из них я передаю при интегрировании (m, n), а другие 2 - это переменные интегрирования (r, z). Это интеграл, который я хочу вычислить:

$$ \int_0^b\int_0^af(r,z)\sin{(\frac{n\pi}{a}z)}J_0(\frac{\alpha_{ 0,м}}{b}r)dzdr $$

Я видел этот пример Выполнение двумерного численного интегрирования с Boost Cpp и замечает, что он использует лямбда-функции для разделения основного интеграла на 2. до сих пор мне это удавалось

double integrate(int m, int n)
{ 
  auto f1 = [](double r, double z, int m, int n) { return integrand(r,z,m,n); };
  auto f = [&](double r, m) {
         auto g = [&](double z, n) {
            return f1(r, z);
              };
         //return gauss_kronrod<double, 61>::integrate(g, 0, a, 5);
         return boost::math::quadrature::trapezoidal(g, 0, a, 1e-6);
       };
  double error;
  //double Q = gauss_kronrod<double, 15>::integrate(f, 0, b, 5, 1e-9, &error);
  double Q =  boost::math::quadrature::trapezoidal(f, 0, b, 1e-6);
  //std::cout << Q << ", error estimated at " << error <<std::endl;
  return Q;
}

Реализация функции $f(r,z)$ и остальной части интеграла следующая

double initial(double r, double z, int m, int n)
{
return std::sin(M_PI*n*z/a)*std::cyl_bessel_j(0, boost::math::cyl_bessel_j_zero(0,m)*r/b);
}
double integrand(double r,double z,int n,int m)
{
  return initial(r,z,m,n)*std::sin(M_PI*n*z/a)*std::cyl_bessel_j(0, boost::math::cyl_bessel_j_zero(0,m)*r/b);
}

Обычно Initial не нуждается в них и n переменных, но в этом случае мне нужно провести некоторые тесты.

Проблема в том, что я действительно не понимаю, как разделить мою функцию, как в примере для моей проблемы, и выполнить интеграцию, потому что boost принимает только 1 переменную функцию.

Пожалуйста помоги


c++ integration
person Carlos Andrés del Valle    schedule 29.04.2021    source источник
comment
Чтобы увидеть, что происходит, создайте отдельные функции для всех лямбд. Кроме того, вам не хватает m и n, когда вы возвращаете f1 из лямбда-функции g.   -  person Biswajit Banerjee    schedule 29.04.2021
comment
Привет, дело в том, что boost принимает только 1 переменную функцию, поэтому g и f должны быть только 1 переменной, не знаю, что делать с m и n.   -  person    schedule 29.04.2021
comment
Это похоже на вопрос программирования о том, как передавать параметры в лямбда-выражения, поэтому здесь он кажется не по теме.   -  person Federico Poloni    schedule 29.04.2021
comment
@CarlosAndrésdelValle: вам не нужно вызывать функцию повышения с двумя аргументами. Просто предварительно умножьте аргумент на m или n, прежде чем отправлять его на повышение.   -  person Biswajit Banerjee    schedule 29.04.2021

Ответы (1)


arrow_upward
1
arrow_downward

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

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

Вот псевдокод:

double integrand(double r,double z, int m,int n, double a, double b)
{
   //this is the function you want to integrate
}

double integrate(int m, int n)
{
    double a=1.0;
    double b=1.0;

    auto inner_integral = [m,n,a,b](double z)
    {
        auto f = [z,m,n,a,b](double r) { return integrand(r,z,m,n,a,b);} 
        return trapezoidal(f,0,a);
    }
    return trapezoidal(inner_integral,0,b);
};

Вам, вероятно, не нужно записывать лямбда-захват, т. е. эталонный захват с &, скорее всего, тоже будет работать (auto inner_integral = [&](double z){...}).

person davidhigh    schedule 29.04.2021
comment
Спасибо за Ваш ответ. У меня вопрос, компилятор кричит о пропаже ; или , перед последним возвратом, и говорит, что функция не имеет оператора возврата. Как я мог это решить? Это что ; напоследок }? - person ; 29.04.2021
comment
Да, и это правильно. После последней лямбда-скобки перед возвратом должен быть символ ;. - person davidhigh; 30.04.2021