qsys_quadrature.h File Reference

Adaptive quadrature for the qsys functions whose MATLAB originals call integral(), and the cumulative trapezoid rule for the one that calls cumtrapz(). More...

#include <cstddef>
#include <vector>
#include "line/num/number.h"
#include "line/util/error.h"

Include dependency graph for qsys_quadrature.h:

Go to the source code of this file.

Namespaces
namespace	line
namespace	line::qsys

Detailed Description

Adaptive quadrature for the qsys functions whose MATLAB originals call integral(), and the cumulative trapezoid rule for the one that calls cumtrapz().

Three MATLAB files in matlab/src/api/qsys reduce a per-class mean response time to a definite integral of a non-elementary integrand – qsys_mg1_fb, qsys_mg1_psjf and qsys_mg1k_loss – and one, qsys_mg1_srpt, integrates on a fixed uniform grid with cumtrapz/trapz. This header carries the two rules they need so that each ported function stays a 1:1 image of its MATLAB file. It is the numerical counterpart of qsys_types.h, which carries the shared return type; nothing here corresponds to a MATLAB file of its own.

The adaptive rule is the Gauss-Kronrod 7/15 pair on recursively bisected subintervals, with the local error estimated as |K - G| and the tolerance split evenly between halves. That is the same rule MATLAB's integral() uses (MATLAB applies it to a transformed interval and controls the error globally), so the two agree to the requested relative tolerance on the smooth integrands here; the ported functions therefore claim agreement at the tolerance MATLAB was asked for, not beyond it.

ARITHMETIC. Both rules are inherently inexact – the Kronrod nodes are irrational and the trapezoid rule has a discretization error – so both are gated on num_traits<T>::has_transcendental. The node and weight constants are the QUADPACK values, but they enter through num_traits<T>::from_double and are therefore carried at double precision: a Real<D> instantiation gains exact accumulation and no cancellation in the sums, but the quadrature error floor stays near 1e-16 relative because the nodes themselves do. Any qsys function that routes through this header inherits that floor, and its tests assert at the tolerance MATLAB's integral() was asked for, never below it.

Definition in file qsys_quadrature.h.