c++ - 如何将多个集成订单整合到我的Integrator类中？

Question

我有一个有效的Integrator类，它将计算单个变量的基本函数的确定积分。我已经测试了一些基本功能的集成，并且看来工作正常。

我现在想扩展此类，以便能够执行相同功能的多个积分……这就是我遇到的障碍……

这是我的集成器类和一些基本用法示例:

Integrator.h

#pragma once

#include <algorithm>
#include <utility>
#include <functional>

struct Limits {
    double lower;
    double upper;

    Limits() : lower{ 0 }, upper{ 0 } {}
    Limits(double a, double b) : lower{ a }, upper{ b } {
        if (a > b) std::swap(lower, upper);
    }

    void applyLimits(double a, double b) {
        lower = a;
        upper = b;
        if (a > b) std::swap(lower, upper);
    }
};

class Integrator {
private:
    Limits limits_;
    std::function<double(double)> integrand_;

    double dx_;
    double dy_;  
    double integral_; 
    int step_size_;

public:
    Integrator(Limits limits, int stepSize, std::function<double(double)> integrand, double dy = 0) 
        : limits_{ limits }, 
        step_size_{ stepSize }, 
        integrand_{ integrand }, 
        dx_{ 0 }, dy_{ 0 } 
    {}
    ~Integrator() = default;

    constexpr double dx() const { return this->dx_; }
    constexpr double dy() const { return this->dy_; }
    constexpr double integral() const { return this->integral_; }

    Limits limits() const { return limits_; }    
    std::function<double(double)>* integrand() { return &this->integrand_; }

    // This is always a 1st order of integration!
    constexpr double evaluate() {
        double distance = limits_.upper - limits_.lower;      // Distance is defined as X0 to XN. (upperLimit - lowerLimit) 
        dx_ = distance / step_size_;                          // Calculate the amount of iterations by dividing 
                                                              // the x-distance by the dx stepsize
        integral_ = 0;                                        // Initialize area to zero
        for (auto i = 0; i < step_size_; i++) {               // For each dx step or iteration calculate the area at Xi
            dy_ = integrand_(limits_.lower + i * dx_);
            double area = dy_ * dx_;                          // Where the width along x is defines as dxStepSize*i 
            integral_ += area;                                // and height(dy) is f(x) at Xi. Sum all of the results
        }

        return integral_;
    }
};

main.cpp

#include <iostream>
#include <exception>
#include <cmath>

#include "Integrator.h"

constexpr double PI = 3.14159265358979;

constexpr double funcA(double x) {
    return x;
}

constexpr double funcB(double x) {
    return (x*x);
}

constexpr double funcC(double x) {
    return ((0.5*(x*x)) + (3*x) - (1/x));
}

double funcD(double x) {
    return sin(x);
}

int main() {
    try {    
        std::cout << "Integration of f(x) = x from a=3.0 to b=5.0\nwith an expected output of 8\n";
        Integrator integratorA(Limits(3.0, 5.0), 10000, &funcA);
        std::cout << integratorA.evaluate() << '\n';        

        std::cout << "\n\nIntegration of f(x) = x^2 from a=2.0 to b=20.0\nwith an expected output of 2664\n";
        Integrator integratorB(Limits(2.0, 20.0), 10000, &funcB);
        std::cout << integratorB.evaluate() << '\n';

        std::cout << "\n\nIntegration of f(x) = (1\\2)x^2 + 3x - (1\\x) from a=1.0 to b=10.0\nwith an expected output of 312.6974\n";
        Integrator integratorC(Limits(1.0, 10.0), 10000, &funcC);
        std::cout << integratorC.evaluate() << '\n';

        std::cout << "\n\nIntegration of f(x) = sin(x) from a=0.0 to b=" <<PI<< "\nwith an expected output of 2\n";
        Integrator integratorD(Limits(0.0, PI), 10000, &funcD);
        std::cout << integratorD.evaluate() << '\n';

    } catch (const std::exception& e) {
        std::cerr << e.what() << std::endl;
        return EXIT_FAILURE;
    }

    return EXIT_SUCCESS;
}

输出

Integration of f(x) = x from a=3.0 to b=5.0
with an expected output of 8
7.9998


Integration of f(x) = x^2 from a=2.0 to b=20.0
with an expected output of 2664
2663.64


Integration of f(x) = (1\2)x^2 + 3x - (1\x) from a=1.0 to b=10.0
with an expected output of 312.6974
312.663


Integration of f(x) = sin(x) from a=0.0 to b=3.14159
with an expected output of 2
2

我当时正在考虑向此类添加另一个函数，类似于evaluate()函数...目前看起来像这样:

double integrate(Limits limits, double dy) {
    double total = 0;
    dy_ = dy;

    for (int i = 0; i < step_size_; i++) {
        double yi = limits_.lower*i*dy_;
        double dx = static_cast<double>(yi - limits.lower) / stepSize;
        double innerArea = 0;

        for (int j = 0; j < step_size_; j++) {
            Integrator inner(limits, step_size_, integrand_, dy_);
            innerArea += inner.evaluate();
        }
        double outerArea = innerArea * dy_;
        total += outerArea;
    }

    integral_ = total;
    return integral_;
}

这就是我感到困惑或困惑的地方...当涉及到内部和外部积分的积分限制时，我不确定如何正确实现此功能。

例如，下面的积分:



每次计算的迭代，内部积分的上限都基于y……这必须动态完成。外部积分从[3,5]而不是[1,y]出发是直截了当的。

我认为我走在正确的轨道上，但是上述方法中的某些方法完全不可用...我从预期或预期值中得到了完全错误的值...

任何和所有的建议和技巧都将受到欢迎!

编辑-注意-我在上面提供了错误的图片，该图片已更新...

预期的输出应为:65.582以及正确提供的函数f(x) = 1/2x^2 + 3x - (1/x)。当我尝试计算双积分时，我最终得到了...

这是驱动程序或main.cpp中添加的代码...

std::cout << "\n\nTesting Double Integration of f(x) = (1\\2)x^2 + 3x - (1\\x) from [3,5] and [1,y]\nwith an expected output of 65.582\n";
Integrator integratorE(Limits(3, 5), 1000, &funcC);
double dy = integratorE.limits().upper - integratorE.limits().lower;
integratorE.integrate(Limits(1, integratorE.integral()), dy);
std::cout << integratorE.integral() << '\n';

但是，它没有在控制台上打印任何内容。

编辑

我没有得到足够的输出，因为我没有等待足够长的时间。迭代由step_size定义为1000。这最终将生成1000^1000总迭代...我在Integrator对象的构造中忽略了这一点。我已在代码中对此进行了更改，使其step_size为100。现在我的应用程序正在输出2.68306e+189的值，这显然是错误的!当我将step_size增加到500时，它给了我6.62804e+190顺序的东西，但这仍然是错误的。

Answer 1

返回并再次观看视频后，...我开始分解类的integrate()函数中的双循环结构。

我从构造函数和此函数的签名中都删除了一些不需要的参数。我删除了传递dy的依赖性，因为我可以在内部计算和存储此值。

我已经对integrate成员函数进行了大修。我现在在适当的时间使用针对dy的适当积分限制来计算dx和step_size。

而不是在此函数中创建Integrator实例并使用该实例的evaluate()函数。我完全删除了此行为，因为我不需要这样做，因为此类存储了一个称为integrand的集成函数实例，其中这是一个std::function<T>对象。有了这个，我可以通过将y传递到该积分中来计算当前的xi。然后，我可以用它来计算求和的内部区域。

我更新的函数如下所示:

double integrate(double lower = 0.0, double upper = 0.0) {
    // Since we are not using the inner upper limit directly
    // make sure that it is still greater than the lower limit
    if (upper <= lower) {
        upper = lower + 1;
    }
    Limits limits(lower, upper);

    double outerSum = 0;
    dy_ = static_cast<double>(limits_.upper - limits_.lower) / step_size_;

    for (int i = 0; i < step_size_; i++) {
        double yi = limits_.lower+i*dy_;
        double dx_ = static_cast<double>(yi - limits.lower) / step_size_;
        double innerSum = 0;

        for (int j = 0; j < step_size_; j++) {
            double xi = limits.lower + dx_ * j;
            double fx = integrand_(xi);                
            double innerArea = fx*dx_;
            innerSum += innerArea;
        }
        double outerArea = innerSum * dy_;
        outerSum += outerArea;
    }

    integral_ = outerSum;
    return integral_;
}

这是我的主类中此函数的用法:

std::cout << "\n\nTesting Double Integration of f(x) = (1\\2)x^2 + 3x - (1\\x) from [3,5] and [1,y]\nwith an expected output of 65.582\n";
Integrator integratorE(Limits(3, 5), 100, &funcC);
//double dy = integratorE.limits().upper - integratorE.limits().lower;
integratorE.integrate(1);
std::cout << integratorE.integral() << '\n';

这给了我以下输出:

Testing Double Integration of f(x) = (1\2)x^2 + 3x - (1\x) from [3,5] and [1,y]
with an expected output of 65.582
64.6426

具有step_size迭代的100和以下内容的输出:

Testing Double Integration of f(x) = (1\2)x^2 + 3x - (1\x) from [3,5] and [1,y]
with an expected output of 65.582
65.3933

具有step_size迭代的500。

因此，就目前该类而言，我可以使用evaluate()对单个变量执行单次确定积分，并且我可以使用integrate(lower,upper)对单个变量进行至少两次定额积分。