Source code for discrete_squared_curvature_derivation_hybrid

import sympy as s
from sympy import init_printing

init_printing()

[docs]def simplify(expr, maxdepth=10, _depth=0):
    import copy
    original_expr = copy.copy(expr)
    print(f"Depth: {_depth} | Parsimony: {len(str(expr))}")
    expr = expr.simplify()
    if expr != original_expr:
        if len(str(expr)) < len(str(original_expr)):
            if _depth < maxdepth:
                return simplify(expr, maxdepth=maxdepth, _depth=_depth + 1)
            else:
                return expr
        else:
            return original_expr
    else:
        return expr
##### Now, a new problem:
##### Do a cubic reconstruction of an interval using values and first derivatives at the endpoints.

# x1, x2 = s.symbols("x1 x2", real=True)
[docs]h = s.symbols("h", real=True)

[docs]x1 = 0
[docs]x2 = h

f1, f2 = s.symbols("f1 f2", real=True)
dfdx1, dfdx2 = s.symbols("dfdx1 dfdx2", real=True)

[docs]q = s.symbols("q", real=True)

[docs]q1 = 0
[docs]q2 = 1
[docs]dqdx = 1 / (x2 - x1)

# Define the Bernstein basis polynomials
[docs]b1 = (1 - q) ** 3
[docs]b2 = 3 * q * (1 - q) ** 2
[docs]b3 = 3 * q ** 2 * (1 - q)
[docs]b4 = q ** 3

c1, c2, c3, c4 = s.symbols('c1 c2 c3 c4', real=True)

# Can solve for c1 and c4 exactly
[docs]c1 = f1
[docs]c4 = f2

f = c1 * b1 + c2 * b2 + c3 * b3 + c4 * b4
dfdx = f.diff(q) * dqdx
# Solve for c2 and c3
[docs]sol = s.solve(
    [
        dfdx.subs(q, q1) - dfdx1,
        dfdx.subs(q, q2) - dfdx2,
    ],
    [
        c2,
        c3,
    ],
)
[docs]c2_sol = simplify(sol[c2])
[docs]c3_sol = simplify(sol[c3])

[docs]f = c1 * b1 + c2_sol * b2 + c3_sol * b3 + c4 * b4

[docs]dfdx = f.diff(q) * dqdx
[docs]df2dx = dfdx.diff(q) * dqdx

res = s.integrate(df2dx ** 2, (q, 0, 1)) * h
[docs]res = simplify(res)

[docs]cse = s.cse(
    res,
    symbols=s.numbered_symbols("s"),
    list=False
)

[docs]parsimony = len(str(res))
print("\nSquared Curvature:")
print(f"Parsimony: {parsimony}")

for i, (var, expr) in enumerate(cse[0]):
    print(f"{var} = {expr}")
print(f"res = {cse[1]}")