численные методы — интегрировать функцию трех переменных переполнение стека
Я потратил некоторое время, просматривая Интернет, чтобы найти решение этой проблемы, может быть, оно есть, но ничто из того, что я увидел, не помогло мне.
У меня есть функция!
double integrand(double r, double phi, double theta)
Это я хочу объединить с некоторыми определенными границами по трем измерениям. Я нашел несколько строк кода в интернете, которые реализуют числовые схемы с определенными переменными и интегралами. Я подумал про себя: «Хорошо, я просто интегрируюсь по одному измерению за другим».
С точки зрения алгоритма, я хотел сделать следующее:
double firstIntegral(double r, double phi) {
double result = integrationFunction(integrand,lower_bound,upper_bound);
return result;
}
И просто сделайте это снова еще два раза. Это легко работает в таких языках, как Matlab, где я могу создать обработчик функций где угодно, но я не знаю, как это сделать в C ++. Сначала мне нужно определить функцию, которую некоторые r и phi будут вычислять подынтегральные функции (r, phi, theta) для любой тэты, и сделать ее в C ++ функцией только одной переменной, но я не знаю, как это сделать.
Как я могу вычислить тройной интеграл моей функции с тремя переменными в C ++, используя одномерную процедуру интеграции (или что-нибудь еще на самом деле …)?
Решение
Это очень медленная и неточная версия для интегралов по декартовым координатам, которая должна работать с C ++ 11.
Он использует std :: function и lambdas для реализации численного интегрирования. Не было предпринято никаких шагов для оптимизации этого.
Решение на основе шаблонов может быть намного быстрее (на несколько порядков), чем это, потому что оно может позволить компилятору встроить и упростить часть кода.
#include<functional>
#include<iostream>
static double integrand(double /*x*/, double y, double /*z*/)
{
return y;
}
double integrate_1d(std::function<double(double)> const &func, double lower, double upper)
{
static const double increment = 0.001;
double integral = 0.0;
for(double x = lower; x < upper; x+=increment) {
integral += func(x) * increment;
}
return integral;
}
double integrate_2d(std::function<double(double, double)> const &func, double lower1, double upper1, double lower2, double upper2)
{
static const double increment = 0.001;
double integral = 0.0;
for(double x = lower2; x < upper2; x+=increment) {
auto func_x = [=](double y){ return func(x, y);};
integral += integrate_1d(func_x, lower1, upper1) * increment;
}
return integral;
}
double integrate_3d(std::function<double(double, double, double)> const &func,
double lower1, double upper1,
double lower2, double upper2,
double lower3, double upper3)
{
static const double increment = 0.001;
double integral = 0.0;
for(double x = lower3; x < upper3; x+=increment) {
auto func_x = [=](double y, double z){ return func(x, y, z);};
integral += integrate_2d(func_x, lower1, upper1, lower2, upper2) * increment;
}
return integral;
}int main()
{
double integral = integrate_3d(integrand, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
std::cout << "Triple integral: " << integral << std::endl;
return 0;
}
Другие решения
Вы можете использовать функторы
#include <iostream>
struct MyFunctorMultiply
{
double m_coeff;
MyFunctorMultiply(double coeff)
{
m_coeff = coeff;
}
double operator()(double value)
{
return m_coeff * value;
}
};
struct MyFunctorAdd
{
double m_a;
MyFunctorAdd(double a)
{
m_a = a;
}
double operator()(double value)
{
return m_a + value;
}
};
template<class t_functor>
double calculate(t_functor functor, double value, double other_param)
{
return functor(value) - other_param;
}
int main()
{
MyFunctorMultiply multiply2(2.);
MyFunctorAdd add3(3.);
double result_a = calculate(multiply2, 4, 1); // should obtain 4 * 2 - 1 = 7
double result_b = calculate(add3, 5, 6); // should obtain 5 + 3 - 6 = 2
std::cout << result_a << std::endl;
std::cout << result_b << std::endl;
}
Если ваша задача заключается только в том, чтобы получить правильный прототип для передачи функции интеграции, вы можете очень хорошо использовать альтернативные механизмы передачи данных, более простой из которых — использование глобальных переменных.
Предполагая, что порядок интеграции включен theta, затем phi, затем r, напишите три функции одного аргумента:
It(theta) вычисляет подынтегральное выражение из аргумента theta прошло явно и глобальный phi а также r,
Ip(phi) вычисляет границы theta из аргумента phi прошло явно и глобальный r; это также копирует phi аргумент глобальной переменной и вызывает integrationFunction(It, lower_t, upper_t),
Ir(r) вычисляет границы phi из аргумента r передано явно; это также копирует r аргумент глобальной переменной и вызывает integrationFunction(Ip, lower_p, upper_p),
Теперь вы готовы позвонить integrationFunction(Ir, lower_r, upper_r),
Это также может быть integrationFunction поддерживает «context«аргумент, где вы можете хранить то, что вы хотите.