Source code for simpson_backward_interpolant_derivation
import sympy as s
from sympy import init_printing
init_printing()
# "Backward" -> gets the integral from x2 to x3, by analogy to backward Euler
# Define the symbols
x1, x2, x3 = s.symbols('x1 x2 x3', real=True)
f1, f2, f3 = s.symbols('f1 f2 f3', real=True)
# Mapping from x-space to q-space has x=x2 -> q=0, x=x3 -> q=1.
# Define the Bernstein basis polynomials
c1, c2, c3 = s.symbols('c1 c2 c3', real=True)
# Can solve for c2 and c3 exactly
f = c1 * b1 + c2 * b2 + c3 * b3
# factors = [f1, f2, f3, f4]
# Solve for c2 and c3
integral = (c1 + c2 + c3) / 3 # God I love Bernstein polynomials
# integral = s.simplify(integral)
print(s.pretty(integral, num_columns=100))
print(f"Parsimony: {parsimony}")