integrate function
double
integrate(
- UnaryFunction<
double> function, - double a,
- double b, {
- int depth = 6,
- double epsilon = 1e-6,
- Iterable<
double> poles = const [], - void onWarning(
- IntegrateWarning type,
- double x
Returns the numerical integration of the provided function
from a
to
b
, that is the result of int(f(x), dx=a..b).
epsilon
is the maximum error tolerance to be accepted. depth
is the
maximum recursion depth, a warning is raised if it is too shallow.
poles
is a list of points at which the function
should not be
evaluated. The integration is automatically split over the generated
intervals: [a..p_1,
p_1..p_2, ...
p_n, b].
In case of an integration problem, onWarning
is evaluated with the
IntegrateWarning and the x position. The default implementation throws
an IntegrateError exception, but a custom handler can continue the
evaluation.
Implementation
double integrate(
UnaryFunction<double> function,
double a,
double b, {
int depth = 6,
double epsilon = 1e-6,
Iterable<double> poles = const [],
void Function(IntegrateWarning type, double x)? onWarning,
}) {
onWarning ??= (type, x) => throw IntegrateError._(type, x);
// Validate boundary condition.
if (a.isNaN) {
throw ArgumentError.value(a, 'a', 'Invalid lower bound');
} else if (b.isNaN) {
throw ArgumentError.value(b, 'b', 'Invalid upper bound');
} else if (a > b) {
return -integrate(function, b, a,
depth: depth, epsilon: epsilon, poles: poles, onWarning: onWarning);
} else if (a == b) {
return 0.0;
}
// Break up at the poles and merge with bounds, if necessary.
if (poles.isNotEmpty) {
final normalized =
poles.where((pole) => a <= pole && pole <= b).toSet().toList();
if (normalized.isNotEmpty) {
normalized.sort();
final expanded = normalized
.expand((pole) => [pole - epsilon, pole + epsilon])
.toList(growable: true);
a < expanded.first ? expanded.insert(0, a) : expanded.removeAt(0);
expanded.last < b ? expanded.add(b) : expanded.removeLast();
var result = 0.0;
for (var i = 0; i < expanded.length; i += 2) {
result += integrate(function, expanded[i], expanded[i + 1],
depth: depth, epsilon: epsilon, onWarning: onWarning);
}
return result;
}
}
// Deal with infinite bounds:
// https://en.wikipedia.org/wiki/Numerical_integration#Integrals_over_infinite_intervals
if (a == double.negativeInfinity && b == double.infinity) {
return integrate(function, a, 0,
depth: depth, epsilon: epsilon, onWarning: onWarning) +
integrate(function, 0, b,
depth: depth, epsilon: epsilon, onWarning: onWarning);
} else if (a == double.negativeInfinity) {
return integrate((t) {
final omt = 1.0 - t, t2 = t * t;
return function(b - omt / t) / t2;
}, epsilon, 1, depth: depth, epsilon: epsilon, onWarning: onWarning);
} else if (b == double.infinity) {
return integrate((t) {
final omt = 1.0 - t, t2 = t * t;
return function(a + omt / t) / t2;
}, epsilon, 1, depth: depth, epsilon: epsilon, onWarning: onWarning);
}
// Solve the actual integral:
// https://en.wikipedia.org/wiki/Adaptive_quadrature
var result = 0.0;
final queue = QueueList.from([
_Quadrature.simpson(
function, depth, epsilon, a, function(a), b, function(b))
]);
while (queue.isNotEmpty) {
final full = queue.removeLast();
final left = _Quadrature.simpson(function, full.depth - 1,
full.epsilon / 2.0, full.a, full.fa, full.m, full.fm);
final right = _Quadrature.simpson(
function, left.depth, left.epsilon, full.m, full.fm, full.b, full.fb);
if (left.epsilon == full.epsilon || left.a == left.m) {
onWarning(IntegrateWarning.doesNotConverge, full.m);
result += full.w;
continue;
}
if (full.depth <= 0) {
onWarning(IntegrateWarning.depthTooShallow, full.m);
}
final delta = left.w + right.w - full.w;
if (full.depth <= 0 || delta.abs() <= 15 * full.epsilon) {
result += left.w + right.w + delta / 15;
} else {
queue.addLast(right);
queue.addLast(left);
}
}
return result;
}