API — ICoreBlocks/ICoreMath/Calculus

The public contract of 3 header(s) under src/ICoreBlocks/ICoreMath/Calculus — 3 class/struct definition(s), 5 declaration(s). Each section shows the header's banner and its public (and protected-virtual) surface exactly as the file writes it.

ICoreNumericalJacobian.h#

src/ICoreBlocks/ICoreMath/Calculus/ICoreNumericalJacobian.h

ICoreNumericalJacobian#

ICoreNumericalJacobian.h:13 · class · 1 declaration(s)

Central-difference gradient of a scalar expression typed as a string, via Eigen::NumericalDiff's functor pattern -- the functor's operator() calls back into ICoreExpressionEvaluator::evaluate() wit...

class ICoreNumericalJacobian {
public:
    // expr: a scalar-valued expression using variables named x1, x2, ..., xn
    //       (n = number of rows/entries in x0).
    // x0: the point to differentiate around (a vector).
    // Returns a 1 x n row vector: d(expr)/d(xi) at x0, for each i.
    static ICoreMatrix gradient(const std::string& expr, const ICoreMatrix& x0);

};
};

ICoreOdeSolver.h#

src/ICoreBlocks/ICoreMath/Calculus/ICoreOdeSolver.h

ICoreOdeSolver#

ICoreOdeSolver.h:50 · class · nested Options · 1 declaration(s)

Core MATLAB's ode45 (the core-MATLAB board, C8.25): the initial-value problem y' = f(t, y), y(t0) = y0, integrated with adaptive step size.

class ICoreOdeSolver {
public:

    // f(t, y) -> dy. False means the call itself failed, which is a console
    // function handle whose body is in error (C2.8); the integration stops and
    // reports rather than reading the failure as a derivative. `dy` must come
    // back with exactly as many entries as `y`.
    using Derivative = std::function<bool(double t, const std::vector<double>& y,
                                          std::vector<double>& dy)>;

    // `odeset`'s options, as ode45 actually reads them -- the five that change
    // the ANSWER and nothing else. Each `has*` flag distinguishes "the user
    // said so" from "left at the default", which for MaxStep and InitialStep
    // is not the same question: MaxStep's default depends on the interval, and
    // an InitialStep that was never given is COMPUTED from y'(t0).
    struct Options {
        double relTol = 1e-3;              // RelTol, MATLAB's default
        std::vector<double> absTol{1e-6};  // AbsTol: one value, or one per equation
        double maxStep     = 0.0;
        bool   hasMaxStep  = false;
        double initialStep = 0.0;
        bool   hasInitialStep = false;
        int    refine      = 4;            // Refine, ode45's default (ode23's is 1)
    };

    // ode45(f, tspan, y0) with `opts`. `tspan` is two endpoints -- integrate
    // between them and report at the solver's own refined points -- or three
    // or more times, which asks for the answer at exactly those times and
    // turns Refine off, as MATLAB's does.
    //
    // `tout` comes back with one entry per output row and `yout` with one ROW
    // per output row, each of `y0.size()` entries: MATLAB's [t, y] shape, so
    // `y(:, k)` is the k-th component on both sides.
    //
    // MATLAB, when a step is refused at the smallest step it will take, warns
    // and answers the part it managed. This answers false with that `t` in the
    // reason instead: the console has no warning channel that survives into a
    // value, and half an integration reported as a whole one is the answer
    // that cannot be right. It is the one place the two deliberately differ.
    static bool solve(const Derivative& f,
                      const std::vector<double>& tspan,
                      const std::vector<double>& y0,
                      const Options& opts,
                      std::vector<double>& tout,
                      std::vector<std::vector<double>>& yout,
                      std::string* whyNot = nullptr);

    // MATLAB's `eps(x)` -- the gap to the next double above |x|, and the
    // smallest denormal at zero. ode45 rebuilds its floor step from it at
    // every step (`tinystep = 16*eps(t)`), so it is part of the algorithm
    // rather than a utility, and it is exposed because odeset's MinStep
    // validation needs the same number.
    static double epsAt(double x);
};
};

ICoreQuadrature.h#

src/ICoreBlocks/ICoreMath/Calculus/ICoreQuadrature.h

ICoreQuadrature#

ICoreQuadrature.h:31 · class · 3 declaration(s)

Core MATLAB's integral (the core-MATLAB board, C8.19): adaptive quadrature over a function you pass in.

class ICoreQuadrature {
public:

    // f(x) -> fx. False means the call itself failed, which is a console
    // function handle whose body is in error (C2.8); the integration stops
    // and reports rather than reading the failure as a number.
    using Integrand = std::function<bool(double x, double& fx)>;

    // f(x, y) -> fxy and f(x, y, z) -> f, for the multiple integrals; and the
    // two shapes a LIMIT of one takes -- MATLAB's inner limits are a number or
    // a function of the outer variables, and both spellings are one callable
    // here so the caller decides once and the algorithm never asks again.
    using Integrand2 = std::function<bool(double x, double y, double& fxy)>;
    using Integrand3 = std::function<bool(double x, double y, double z, double& f)>;
    using Curve      = std::function<bool(double x, double& y)>;
    using Surface    = std::function<bool(double x, double y, double& z)>;

    // integral(f, a, b) with MATLAB's defaults, or a tolerance of your own:
    // absTol 1e-10 and relTol 1e-6 are what `integralParseArgs` sets. Either
    // limit may be infinite, and a > b answers the negated integral, as
    // MATLAB's does.
    //
    // `initialIntervalCount` is `integralParseArgs`'s own field and is 10 for
    // every `integral` call. It is a parameter because MATLAB sets it to
    // **3** for the outer integral of `integral3` and for both integrals of
    // `integral2`'s iterated method -- a different first mesh is a different
    // answer, so the caller that needs 3 has to be able to say so.
    static bool integrate(const Integrand& f, double a, double b,
                          double absTol, double relTol, double& value,
                          std::string* whyNot = nullptr,
                          int initialIntervalCount = 10);

    // quadgk(f, a, b) -- the SAME Gauss-Kronrod (7,15) rule as `integral` and
    // a different loop around it: `quadgk.m`'s `vadapt` accepts subintervals
    // and stops, where `integralCalc.m` can decide its early tolerance was
    // never real and rebuild the mesh from the start. The two therefore
    // answer differently in their last digits, which is why this is a second
    // entry point and not `integrate` under another name.
    //
    // `waypoints` are the interior points MATLAB's Waypoints option puts in
    // the first mesh -- an integrand's kinks, so the mesh never straddles one.
    // `maxIntervalCount` is MATLAB's 650. `errbnd` is quadgk's second output.
    static bool gaussKronrod(const Integrand& f, double a, double b,
                             const std::vector<double>& waypoints,
                             double absTol, double relTol, int maxIntervalCount,
                             double& value, double& errbnd,
                             std::string* whyNot = nullptr);

    // quad(f, a, b, tol) -- adaptive SIMPSON with one Romberg step, and
    // quadl(f, a, b, tol) -- adaptive Lobatto with a Kronrod refinement.
    // MATLAB keeps both and recommends neither; they are here because a `.m`
    // that calls one has to cross, and they are not `integral` in disguise:
    // each has its own first mesh, its own local error test and its own
    // function-count ceiling, so each answers its own number.
    static bool adaptiveSimpson(const Integrand& f, double a, double b, double tol,
                                double& value, std::string* whyNot = nullptr);
    static bool adaptiveLobatto(const Integrand& f, double a, double b, double tol,
                                double& value, std::string* whyNot = nullptr);

    // Which two-dimensional algorithm runs. MATLAB's 'auto' is not a third
    // algorithm: it is Tiled for finite limits and Iterated for an infinite
    // one, and asking for Tiled with an infinite limit is an ERROR there
    // rather than a quiet fall back to Iterated.
    enum class Method { Auto, Tiled, Iterated };

    // integral2(f, xmin, xmax, ymin, ymax): the y limits are constants or
    // functions of x. Tiled is Shampine's TwoD -- a (3,7) Gauss-Kronrod
    // tensor product over a rectangle list ordered by adjusted error, under a
    // COSINE substitution in both directions that flattens all four edges.
    //
    // `improperLimits` is MATLAB's own `isImproper`, and the caller computes it
    // because only the caller can: it is true when a limit GIVEN AS A CONSTANT
    // is infinite, and a limit given as a function of x is never improper
    // however large it grows (`integral2.m` lines 25-45). It decides what
    // Method::Auto resolves to, and makes an explicit Method::Tiled an error.
    static bool integrate2(const Integrand2& f, double xmin, double xmax,
                           const Curve& ymin, const Curve& ymax,
                           bool improperLimits,
                           double absTol, double relTol, Method method,
                           double& value, std::string* whyNot = nullptr);

    // integral3(f, xmin, xmax, ymin, ymax, zmin, zmax): always the outer
    // `integral` over x of an `integral2` in (y, z), whatever `method` says --
    // `method` steers that inner two-dimensional call and nothing else.
    static bool integrate3(const Integrand3& f, double xmin, double xmax,
                           const Curve& ymin, const Curve& ymax,
                           const Surface& zmin, const Surface& zmax,
                           bool improperLimits,
                           double absTol, double relTol, Method method,
                           double& value, std::string* whyNot = nullptr);

    static double defaultAbsTol();
    static double defaultRelTol();

    // quad/quadl's own default, which is neither of the two above: one
    // tolerance of 1e-6, tested against the LOCAL error of a subinterval.
    static double defaultLegacyTol();
};
};