from typing import Union
import casadi as _cas
import numpy as _onp
from aerosandbox.numpy.array import length, roll, concatenate
from aerosandbox.numpy.calculus import diff
[docs]def integrate_discrete_intervals(
f: Union[_onp.ndarray, _cas.MX],
x: Union[_onp.ndarray, _cas.MX] = None,
multiply_by_dx: bool = True,
method: str = "trapezoidal",
method_endpoints: str = "lower_order",
):
"""
Given a set of sampled points (x_i, f_i) from a function, computes the integral of that function over each set of
adjacent points ("intervals"). Does this via a reconstruction approach, with several methods available.
In general, N points will yield N-1 integrals (one for each "interval" between points).
Args:
f: A 1D array of function values.
x: A 1D array of x-values where the function was evaluated. If not specified, defaults to the indices of f.
Should be the same length as f and should be monotonically increasing (i.e. x[i] < x[i+1], with no duplicated points).
multiply_by_dx: Whether to multiply the integral by the width of the segment. Defaults to True.
- If True, summing the integrals will yield the integral of the function over the entire domain (x[0] to x[-1])
- If False, you can think of the output as the "average function value" over each interval.
method: The integration method to use. Options are:
- "forward_euler"
- "backward_euler"
- "trapezoidal" (default)
- "forward_simpson"
- "backward_simpson"
- "cubic"
Note that some methods, like "cubic", approximate each segment interval by looking beyond just the integral itself (i.e., f(a) and f(b)),
and so are not possible near the endpoints of the array.
method_endpoints: The integration method to use at the endpoints, for those higher-order methods that require handling. Options are:
- "lower_order" (default)
- "ignore" (i.e. return the integral of the interior points only - note that this may result in a different number of integrals than segments!)
- "periodic"
"""
# Determine if an x-array was specified, and calculate dx.
x_is_specified = (x is not None)
if not x_is_specified:
x = _onp.arange(length(f))
dx = x[1:] - x[:-1]
method = str(method).lower().replace(" ", "_")
# Implement integration methods
if method in ["forward_euler", "forward", "euler_forward", "left", "left_riemann"]:
avg_f = f[:-1]
degree = 0 # Refers to the highest degree of the polynomial that the method is exact for.
remaining_endpoint_intervals = (0, 0)
elif method in ["backward_euler", "backward", "euler_backward", "right", "right_riemann"]:
avg_f = f[1:]
degree = 0
remaining_endpoint_intervals = (0, 0)
elif method in ["trapezoidal", "trapezoid", "trapz", "midpoint"]:
if method == "midpoint":
raise PendingDeprecationWarning(
"The 'midpoint' method will be deprecated at a future point, since 'trapezoidal' is the more accurate term here.")
avg_f = (f[1:] + f[:-1]) / 2
degree = 1
remaining_endpoint_intervals = (0, 0)
elif method in ["forward_simpson", "simpson_forward", "simpson"]:
x1 = x[:-2]
x2 = x[1:-1]
x3 = x[2:]
f1 = f[:-2]
f2 = f[1:-1]
f3 = f[2:]
h = x2 - x1
hp = x3 - x2
# q1 = 0 # Integration lower bound
# q2 = 1 # Integration upper bound
q3 = 1 + hp / h
avg_f = (
f1 - f3
+ 3 * q3 ** 2 * (f1 + f2)
- 2 * q3 * (2 * f1 + f2)
) / (
6 * q3 * (q3 - 1)
)
degree = 2
remaining_endpoint_intervals = (0, 1)
elif method in ["backward_simpson", "simpson_backward"]:
x1 = x[:-2]
x2 = x[1:-1]
x3 = x[2:]
f1 = f[:-2]
f2 = f[1:-1]
f3 = f[2:]
h = x3 - x2
hm = x2 - x1
q1 = -hm / h
# q2 = 0 # Integration lower bound
# q3 = 1 # Integration upper bound
avg_f = (
f2 - f1
+ 3 * q1 ** 2 * (f2 + f3)
- 2 * q1 * (2 * f2 + f3)
) / (
6 * q1 * (q1 - 1)
)
degree = 2
remaining_endpoint_intervals = (1, 0)
elif method in ["cubic", "cubic_spline"]:
x1 = x[:-3]
x2 = x[1:-2]
x3 = x[2:-1]
x4 = x[3:]
f1 = f[:-3]
f2 = f[1:-2]
f3 = f[2:-1]
f4 = f[3:]
h = x3 - x2
hm = x2 - x1
hp = x4 - x3
q1 = -hm / h
# q2 = 0 # Integration lower bound
# q3 = 1 # Integration upper bound
q4 = 1 + hp / h
avg_f = (
6 * q1 ** 3 * q4 ** 2 * (f2 + f3)
- 4 * q1 ** 3 * q4 * (2 * f2 + f3)
+ 2 * q1 ** 3 * (f2 - f4)
- 6 * q1 ** 2 * q4 ** 3 * (f2 + f3)
+ 3 * q1 ** 2 * q4 * (3 * f2 + f3)
+ 3 * q1 ** 2 * (f4 - f2)
+ 4 * q1 * q4 ** 3 * (2 * f2 + f3)
- 3 * q1 * q4 ** 2 * (3 * f2 + f3)
+ q1 * (f2 - f4)
+ 2 * q4 ** 3 * (f1 - f2)
+ 3 * q4 ** 2 * (f2 - f1)
+ q4 * (f1 - f2)
) / (
12 * q1 * q4 * (q1 - 1) * (q1 - q4) * (q4 - 1)
)
degree = 3
remaining_endpoint_intervals = (1, 1)
else:
raise ValueError(f"Invalid method '{method}'.")
if method_endpoints == "lower_order":
if degree >= 3:
method_endpoints = "simpson"
else:
method_endpoints = "trapezoidal"
if method_endpoints == "simpson":
if remaining_endpoint_intervals[0] != 0:
avg_f_left_intervals = integrate_discrete_intervals(
f=f[:2 + remaining_endpoint_intervals[0]],
x=x[:2 + remaining_endpoint_intervals[0]],
multiply_by_dx=False,
method="forward_simpson",
method_endpoints="ignore",
)
avg_f = concatenate((
avg_f_left_intervals,
avg_f
))
if remaining_endpoint_intervals[1] != 0:
avg_f_right_intervals = integrate_discrete_intervals(
f=f[-(2 + remaining_endpoint_intervals[1]):],
x=x[-(2 + remaining_endpoint_intervals[1]):],
multiply_by_dx=False,
method="backward_simpson",
method_endpoints="ignore",
)
avg_f = concatenate((
avg_f,
avg_f_right_intervals,
))
elif method_endpoints == "trapezoidal":
if remaining_endpoint_intervals[0] != 0:
avg_f_left_intervals = integrate_discrete_intervals(
f=f[:1 + remaining_endpoint_intervals[0]],
x=x[:1 + remaining_endpoint_intervals[0]],
multiply_by_dx=False,
method="trapezoidal",
method_endpoints="ignore",
)
avg_f = concatenate((
avg_f_left_intervals,
avg_f
))
if remaining_endpoint_intervals[1] != 0:
avg_f_right_intervals = integrate_discrete_intervals(
f=f[-(1 + remaining_endpoint_intervals[1]):],
x=x[-(1 + remaining_endpoint_intervals[1]):],
multiply_by_dx=False,
method="trapezoidal",
method_endpoints="ignore",
)
avg_f = concatenate((
avg_f,
avg_f_right_intervals,
))
else:
raise ValueError(f"Invalid method_endpoints '{method_endpoints}'.")
elif method_endpoints == "ignore":
pass
elif method_endpoints == "periodic":
raise NotImplementedError("Periodic integration is not yet implemented.")
else:
raise ValueError(f"Invalid method_endpoints '{method_endpoints}'.")
if multiply_by_dx:
if x_is_specified:
return avg_f * dx
else:
return avg_f
else:
return avg_f
[docs]def integrate_discrete_squared_curvature(
f: Union[_onp.ndarray, _cas.MX],
x: Union[_onp.ndarray, _cas.MX] = None,
method: str = "hybrid_simpson_cubic",
):
"""
Given a set of sampled points (x_i, f_i) from a function f(x), computes the following quantity:
int_{x[0]}^{x[-1]} (f''(x))^2 dx
This is useful for regularization of smooth curves (i.e., encouraging smooth functions as optimization results).
Performs this through one of several reconstruction-based methods, specified by `method`:
* "cubic": On each interval, reconstructs a piecewise cubic polynomial. This cubic is the unique polynomial
that passes through the two points at the endpoints of the interval, plus the next point beyond each endpoint
of the interval (i.e., 4 points in total). Numerically, this cubic is obtained using Bernstein polynomial
reconstruction, so it is numerically stable. This cubic is then analytically differentiated twice, squared,
and integrated over the interval. At the ends of the overall array, where this "look beyond" strategy is not
possible, a one-sided cubic is used instead (i.e., looks beyond the interval at one end only and uses two
extra points from this side).
* "simpson": On each interval, makes two unique quadratic reconstructions:
* One reconstruction that uses the two points at the endpoints of the interval, plus the next point beyond
the right endpoint of the interval (i.e., 3 points in total).
* One reconstruction that uses the two points at the endpoints of the interval, plus the next point beyond
the left endpoint of the interval (i.e., 3 points in total).
These two quadratics are then analytically differentiated twice, squared, and integrated over the
interval. This requires much less calculation, since the quadratics have uniform curvature over the
interval, causing a lot of things to simplify. The result is then computed by combining the results of this
process for the two quadratic reconstructions.
This is similar to a Simpson's rule integration, balanced between the two sides of the interval. In
frequency-domain testing, this method appears to be more accurate than the "cubic" strategy at every
frequency, with less computational effort. Thus, it should be preferred to the "cubic" strategy.
* "hybrid_simpson_cubic": First, starts out by estimating the first derivative of the function at each point
in the array (including endpoints) using a quadratic reconstruction. (See `numpy.gradient()` for more
information or source code on this; this code uses `numpy.gradient()` directly for this step.) Then,
reconstructs a cubic polynomial on each interval, with the following boundary conditions:
* The cubic passes through the two points at the endpoints of the interval.
* The cubic has the same first derivative as the precomputed derivatives at the endpoints of the interval.
This cubic is then analytically differentiated twice, squared, and integrated over the interval.
In frequency-domain testing, this method is also more accurate than the "cubic" strategy at every
frequency. Compared to the "simpson" strategy, it is more accurate at high frequencies and less accurate
at low frequencies. Because the goal of this function is to be used as a regularization term,
which should be more sensitive to high-frequency oscillations, this method is preferred to the "simpson"
strategy. This method is also preferred as its estimate tends to err high rather than low, which serves
well as a regularization strategy. (It is still convergent to the true value in the high-sample-rate limit.)
"""
# Determine if an x-array was specified, and calculate dx.
x_is_specified = (x is not None)
if not x_is_specified:
x = _onp.arange(length(f))
if method in ["cubic", "cubic_spline"]:
x1 = x[:-3]
x2 = x[1:-2]
x3 = x[2:-1]
x4 = x[3:]
f1 = f[:-3]
f2 = f[1:-2]
f3 = f[2:-1]
f4 = f[3:]
h = x3 - x2
hm = x2 - x1
hp = x4 - x3
dfm = f2 - f1
df = f3 - f2
dfp = f4 - f3
### The following section computes the integral of the squared second derivative of the cubic spline interpolant
### for the "middle" intervals (i.e. not the first or last intervals).
### Code is generated by sympy; here s_i variables represent common subexpressions.
s0 = hm ** 2
s1 = hp ** 2
s2 = h + hm
s3 = h ** 2
s4 = hp ** 6
s5 = h ** 6
s6 = hp ** 5
s7 = h ** 3
s8 = 3 * s7
s9 = hp ** 4
s10 = h ** 4
s11 = 4 * s10
s12 = hp ** 3
s13 = 3 * h ** 5
s14 = hm ** 6
s15 = hm ** 5
s16 = hm ** 4
s17 = hm ** 3
s18 = hm * s10 * s12
s19 = hp * s10 * s17
s20 = s12 * s17
s21 = s0 * s12 * s8
s22 = s1 * s17 * s8
s23 = 3 * s0 * s1 * s10 + s21 + s22
s24 = 2 * h
s25 = 6 * s3
s26 = 7 * s7
s27 = 3 * s3
s28 = -s20 * s27
s29 = 3 * h
middle_intervals = 4 * (df ** 2 * (
s0 * s27 * s9 + s0 * s29 * s6 + s0 * s4 + s1 * s14 + s1 * s15 * s29 + s1 * s16 * s27 - s12 * s16 * s29 - 2 * s16 * s9 - s17 * s29 * s9 + s23 + s28) + df * dfm * (
2 * h * s17 * s9 - hm * s1 * s13 - hm * s24 * s4 - hm * s25 * s6 - hm * s26 * s9 + 3 * s0 * s3 * s9 + s1 * s17 * s7 - 6 * s18 + s21 - s28) + df * dfp * (
2 * h * s12 * s16 - hp * s0 * s13 - hp * s14 * s24 - hp * s15 * s25 - hp * s16 * s26 + s0 * s12 * s7 + 3 * s1 * s16 * s3 - 6 * s19 + s22 - s28) + dfm ** 2 * (
s1 * s5 + s11 * s9 + s12 * s13 + s3 * s4 + s6 * s8) + dfm * dfp * (
hm * hp * s5 - s18 - s19 - 2 * s20 * s3 - s23) + dfp ** 2 * (
s0 * s5 + s11 * s16 + s13 * s17 + s14 * s3 + s15 * s8)) / (
h * s0 * s1 * s2 ** 2 * (h + hp) ** 2 * (hp + s2) ** 2)
### Now we compute the integral for the first interval.
h_f = h[slice(0, 1)]
hm_f = hm[slice(0, 1)]
hp_f = hp[slice(0, 1)]
df_f = df[slice(0, 1)]
dfm_f = dfm[slice(0, 1)]
dfp_f = dfp[slice(0, 1)]
s0 = h_f ** 2
s1 = hp_f ** 2
s2 = h_f + hm_f
s3 = hp_f ** 6
s4 = df_f * dfm_f
s5 = hm_f ** 6
s6 = 2 * dfp_f
s7 = df_f * s6
s8 = h_f ** 6
s9 = dfm_f * dfp_f
s10 = 4 * s9
s11 = df_f ** 2
s12 = hm_f ** 2
s13 = dfm_f ** 2
s14 = dfp_f ** 2
s15 = hm_f ** 3
s16 = hp_f ** 5
s17 = 3 * s11
s18 = hp_f ** 4
s19 = hm_f ** 4
s20 = 4 * s19
s21 = hm_f ** 5
s22 = hp_f ** 3
s23 = h_f ** 3
s24 = 6 * s13
s25 = h_f ** 4
s26 = h_f ** 5
s27 = 3 * s14
s28 = df_f * dfp_f
s29 = -dfm_f
s30 = df_f * h_f
s31 = 15 * df_f
s32 = 2 * df_f
s33 = dfm_f * hm_f
s34 = 3 * dfm_f
s35 = -s34
s36 = dfp_f * hp_f
s37 = 3 * s12
s38 = 5 * s11
s39 = 3 * s0
s40 = -3 * s28 + s9
s41 = 18 * s11
s42 = -s7
first_interval = 4 * (
-2 * h_f * hm_f * s3 * s4 - h_f * hp_f * s5 * s7 + hm_f * hp_f * s10 * s8 - hm_f * s0 * s16 * s34 * (
4 * df_f + s29) - 9 * hp_f * s0 * s21 * s28 + s0 * s13 * s3 + s0 * s14 * s5 + s0 * s15 * s22 * (
27 * s11 - 21 * s4 + s40) + s1 * s11 * s5 + 4 * s1 * s13 * s8 + s1 * s15 * s23 * (
-12 * s28 - 14 * s4 + s41 + 9 * s9) + s1 * s19 * s39 * (
s38 - s4 + s40) + 3 * s1 * s21 * s30 * (-dfp_f + s32) + s1 * s25 * s37 * (
s10 + s13 + s17 - 7 * s4 + s42) + 6 * s1 * s26 * s33 * (
dfm_f + dfp_f - s32) + s11 * s12 * s3 + s11 * s18 * s20 + s12 * s14 * s8 + 6 * s12 * s16 * s30 * (
df_f + s29) + s12 * s18 * s39 * (s13 + s38 - 9 * s4) + s12 * s22 * s23 * (
s24 - 42 * s4 + s41 + s42 + 3 * s9) + 13 * s13 * s18 * s25 + 12 * s13 * s22 * s26 + s14 * s20 * s25 + s15 * s16 * s17 + s15 * s18 * s30 * (
-7 * dfm_f + s31) - s15 * s25 * s36 * (
-8 * dfm_f + s31) + s15 * s26 * s27 + s16 * s23 * s24 + s17 * s21 * s22 - 4 * s18 * s23 * s33 * (
7 * df_f + s35) - s19 * s22 * s30 * (dfp_f - s31 + s34) - s19 * s23 * s36 * (
16 * df_f + s35) + s21 * s23 * s27 + s22 * s25 * s33 * (
-30 * df_f + 15 * dfm_f + s6) - s26 * s36 * s37 * (s32 + s35)) / (
hm_f * s0 * s1 * s2 ** 2 * (h_f + hp_f) ** 2 * (hp_f + s2) ** 2)
### Now we compute the integral for the last interval.
h_l = h[slice(-1, None)]
hm_l = hm[slice(-1, None)]
hp_l = hp[slice(-1, None)]
df_l = df[slice(-1, None)]
dfm_l = dfm[slice(-1, None)]
dfp_l = dfp[slice(-1, None)]
s0 = h_l ** 2
s1 = hm_l ** 2
s2 = h_l + hm_l
s3 = hp_l ** 6
s4 = 2 * dfm_l
s5 = df_l * s4
s6 = hm_l ** 6
s7 = df_l * dfp_l
s8 = h_l ** 6
s9 = dfm_l * dfp_l
s10 = 4 * s9
s11 = df_l ** 2
s12 = hp_l ** 2
s13 = dfm_l ** 2
s14 = dfp_l ** 2
s15 = hm_l ** 3
s16 = hp_l ** 5
s17 = 3 * s11
s18 = hm_l ** 4
s19 = hp_l ** 4
s20 = 4 * s19
s21 = hm_l ** 5
s22 = hp_l ** 3
s23 = h_l ** 3
s24 = 3 * s13
s25 = h_l ** 4
s26 = h_l ** 5
s27 = 6 * s14
s28 = df_l * dfm_l
s29 = -dfp_l
s30 = df_l * h_l
s31 = 15 * df_l
s32 = 2 * df_l
s33 = dfp_l * hp_l
s34 = 3 * dfp_l
s35 = -s34
s36 = dfm_l * hm_l
s37 = 3 * s12
s38 = 5 * s11
s39 = 3 * s0
s40 = -3 * s28 + s9
s41 = 18 * s11
s42 = -s5
last_interval = 4 * (
-h_l * hm_l * s3 * s5 - 2 * h_l * hp_l * s6 * s7 + hm_l * hp_l * s10 * s8 - 9 * hm_l * s0 * s16 * s28 - hp_l * s0 * s21 * s34 * (
4 * df_l + s29) + s0 * s13 * s3 + s0 * s14 * s6 + s0 * s15 * s22 * (
27 * s11 + s40 - 21 * s7) + s1 * s11 * s3 + 4 * s1 * s14 * s8 + 3 * s1 * s16 * s30 * (
-dfm_l + s32) + s1 * s19 * s39 * (s38 + s40 - s7) + s1 * s22 * s23 * (
-12 * s28 + s41 - 14 * s7 + 9 * s9) + s1 * s25 * s37 * (
s10 + s14 + s17 + s42 - 7 * s7) + 6 * s1 * s26 * s33 * (
dfm_l + dfp_l - s32) + s11 * s12 * s6 + s11 * s18 * s20 + s12 * s13 * s8 + s12 * s15 * s23 * (
s27 + s41 + s42 - 42 * s7 + 3 * s9) + s12 * s18 * s39 * (
s14 + s38 - 9 * s7) + 6 * s12 * s21 * s30 * (
df_l + s29) + s13 * s20 * s25 + 12 * s14 * s15 * s26 + 13 * s14 * s18 * s25 + s15 * s16 * s17 - s15 * s19 * s30 * (
dfm_l - s31 + s34) + s15 * s25 * s33 * (
-30 * df_l + 15 * dfp_l + s4) + s16 * s23 * s24 + s17 * s21 * s22 + s18 * s22 * s30 * (
-7 * dfp_l + s31) - 4 * s18 * s23 * s33 * (7 * df_l + s35) - s19 * s23 * s36 * (
16 * df_l + s35) + s21 * s23 * s27 + s22 * s24 * s26 - s22 * s25 * s36 * (
-8 * dfp_l + s31) - s26 * s36 * s37 * (s32 + s35)) / (
hp_l * s0 * s1 * s2 ** 2 * (h_l + hp_l) ** 2 * (hp_l + s2) ** 2)
### Now, we stitch together the intervals.
res = concatenate((
first_interval,
middle_intervals,
last_interval,
))
return res
elif method in ["simpson"]:
### Forward Simpson for intervals 0 to N-2
x2 = x[:-2]
x3 = x[1:-1]
x4 = x[2:]
f2 = f[:-2]
f3 = f[1:-1]
f4 = f[2:]
h = x3 - x2
hp = x4 - x3
df = f3 - f2
dfp = f4 - f3
res_forward_simpson = 4 * (df * hp - dfp * h) ** 2 / (h * hp ** 2 * (h + hp) ** 2)
### Backward Simpson for intervals 1 to N-1
x1 = x[:-2]
x2 = x[1:-1]
x3 = x[2:]
f1 = f[:-2]
f2 = f[1:-1]
f3 = f[2:]
h = x3 - x2
hm = x2 - x1
dfm = f2 - f1
df = f3 - f2
res_backward_simpson = 4 * (df * hm - dfm * h) ** 2 / (h * hm ** 2 * (h + hm) ** 2)
### Fuse them together
first_interval = res_forward_simpson[slice(0, 1)]
a = res_backward_simpson[slice(None, -1)]
b = res_forward_simpson[slice(1, None)]
# middle_intervals = (a + b) / 2
middle_intervals = ((a ** 2 + b ** 2) / 2 + 1e-100) ** 0.5 # This is more accurate across all frequencies
last_interval = res_backward_simpson[slice(-1, None)]
res = concatenate((
first_interval,
middle_intervals,
last_interval,
))
return res
elif method in ["hybrid_simpson_cubic"]:
from aerosandbox.numpy.calculus import gradient
dfdx = gradient(
f,
x,
edge_order=2
)
h = x[1:] - x[:-1]
df = f[1:] - f[:-1]
dfdx1 = dfdx[:-1]
dfdx2 = dfdx[1:]
res = (
4 * (dfdx1 ** 2 + dfdx1 * dfdx2 + dfdx2 ** 2) / h
+ 12 * df / h ** 2 * (df / h - dfdx1 - dfdx2)
)
return res
else:
raise ValueError(f"Invalid method '{method}'.")
if __name__ == '__main__':
import aerosandbox as asb
import aerosandbox.numpy as np
from scipy import integrate, interpolate
import sympy as s
np.random.seed(0)
# degree = 4
# coeffs = np.random.randn(degree + 1)
# a = 1
# b = 3
#
#
# def f(x):
# out = 0
# for i in range(degree + 1):
# out += coeffs[i] * x ** i
#
# return out
b = 2
def f(x):
sin = np.sin if isinstance(x, np.ndarray) else s.sin
return sin(2 * np.pi * x * 1) + 1
print("\n\nTest 1: Integration")
exact = integrate.quad(
f,
a,
b,
epsrel=1e-15,
)[0]
x_vals = np.cosspace(a, b, 100)
f_vals = f(x_vals)
approx_intervals = integrate_discrete_intervals(
f=f_vals,
x=x_vals,
multiply_by_dx=True,
# method="trapz",
method="cubic",
# method_endpoints="ignore",
)
integral = np.sum(approx_intervals)
print(f"exact: {exact}")
print(f"approx: {integral}")
print(f"error: {integral - exact}")
print("\n\nTest 2: Squared curvature")
x_sym = s.symbols("x")
f_sym = f(x_sym)
df2dx_func = s.lambdify(x_sym, s.diff(f_sym, x_sym, 2))
exact = integrate.quad(
lambda x: df2dx_func(x) ** 2,
a,
b,
epsrel=1e-15,
)[0]
print(f"exact: {exact}")
approx = integrate_discrete_squared_curvature(
f=f_vals,
x=x_vals,
# method="simpson"
method="hybrid_simpson_cubic",
)
integral = np.sum(approx)
print(f"\nintegrate_discrete_squared_curvature: {integral}")
print(f"error: {integral - exact}")
approx = integrate_discrete_intervals(
f=np.gradient(
f_vals,
x_vals,
n=2
) ** 2,
x=x_vals,
)
integral = np.sum(approx)
print(f"\nintegrate_discrete_intervals + np.gradient: {integral}")
print(f"error: {integral - exact}")
print("\n\nTest 3: Squared curvature with global-spline reconstruction")
x = np.arange(0, 100 + 1)
f = np.cos(np.pi * x / 2)
f_interp = interpolate.InterpolatedUnivariateSpline(
x=x,
y=f,
k=3
)
exact = integrate.quad(
lambda x: f_interp.derivative(2)(x) ** 2,
x[0],
x[-1],
epsrel=1e-8,
)[0]
approx = integrate_discrete_squared_curvature(
f=f,
x=x,
method="hybrid_simpson_cubic",
# method="simpson"
)
integral = np.sum(approx)
print(f"exact: {exact}")
print(f"\nintegrate_discrete_squared_curvature: {integral}")
print(f"error: {integral - exact}")
x_plot = np.linspace(x[0], x[-1], 10000)
import matplotlib.pyplot as plt
import aerosandbox.tools.pretty_plots as p
# p.qp(x_plot, f_interp.derivative(2)(x_plot) ** 2)