From da584d7ca7f128d65d70ff6ebe2bd3c54e881440 Mon Sep 17 00:00:00 2001
From: Peter Sharpe <peterdsharpe@gmail.com>
Date: Mon, 1 Jan 2024 19:48:32 +0100
Subject: [PATCH] update derivations

---
 .../quadratic_2nd_derivative.py               | 104 ++++++++++++
 ..._squared_curvature_derivation_bernstein.py | 159 ++++++++++++++++++
 ...e_squared_curvature_derivation_lagrange.py | 123 ++++++++++++++
 3 files changed, 386 insertions(+)
 create mode 100644 aerosandbox/numpy/derivative_discrete_derivations/quadratic_2nd_derivative.py
 create mode 100644 aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_derivation_bernstein.py
 create mode 100644 aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_derivation_lagrange.py

diff --git a/aerosandbox/numpy/derivative_discrete_derivations/quadratic_2nd_derivative.py b/aerosandbox/numpy/derivative_discrete_derivations/quadratic_2nd_derivative.py
new file mode 100644
index 00000000..8c820151
--- /dev/null
+++ b/aerosandbox/numpy/derivative_discrete_derivations/quadratic_2nd_derivative.py
@@ -0,0 +1,104 @@
+import sympy as s
+from sympy import init_printing
+init_printing()
+
+# Reconstructs a quadratic interpolant from x1...x3, then gets the derivative at x2
+
+# Define the symbols
+x1, x2, x3 = s.symbols('x1 x2 x3', real=True)
+f1, f2, f3 = s.symbols('f1 f2 f3', real=True)
+
+# hm = x2 - x1
+# hp = x3 - x2
+hm, hp = s.symbols("hm hp")
+
+q = s.symbols('q')  # Normalized space for a Bernstein basis.
+# Mapping from x-space to q-space has x=x2 -> q=0, x=x3 -> q=1.
+
+q1 = 0
+# q2 = s.symbols('q2', real=True) # (x2 - x1) / (x3 - x1)
+q2 = hm / (hm + hp)
+q3 = 1
+
+# Define the Bernstein basis polynomials
+b1 = (1 - q) ** 2
+b2 = 2 * q * (1 - q)
+b3 = q ** 2
+
+c1, c2, c3 = s.symbols('c1 c2 c3', real=True)
+
+# Can solve for c2 and c3 exactly
+c1 = f1
+c3 = f3
+
+f = c1 * b1 + c2 * b2 + c3 * b3
+
+f2_quadratic = f.subs(q, q2)#.simplify()
+
+factors = [q2]
+# factors = [f1, f2, f3, f4]
+
+# Solve for c2 and c3
+sol = s.solve(
+    [
+        f2_quadratic - f2,
+    ],
+    [
+        c2,
+    ],
+)
+c2 = sol[c2].factor(factors).simplify()
+
+
+f = c1 * b1 + c2 * b2 + c3 * b3
+dfdq = f.diff(q).simplify()
+# dqdx = 1 / (x3 - x1)
+dqdx = 1 / (x3 - x1)
+dfdx = dfdq * dqdx
+
+df2dx = (
+    dfdx.diff(q) / (x3 - x1)
+)
+
+dfm, dfp = s.symbols("dfm dfp")
+
+def simplify(expr):
+    import copy
+    original_expr = copy.copy(expr)
+    expr = expr.subs({
+        f3 - f2: dfp,
+        f2 - f1: dfm,
+        f3 - f1: dfp + dfm,
+        x3 - x1: hm + hp,
+    })
+    expr = expr.subs({
+        f3 - f2: dfp,
+        f2 - f1: dfm,
+        f3 - f1: dfp + dfm,
+        x3 - x1: hm + hp,
+        f1 - 2 * f2 + f3: dfp - dfm,
+    })
+    expr = expr.factor([
+        hp,
+        hm
+    ]).simplify()
+    if expr != original_expr:
+        expr = simplify(expr)
+    return expr
+
+dfdx_q1 = simplify(dfdx.subs(q, q1))
+dfdx_q2 = simplify(dfdx.subs(q, q2))
+dfdx_q3 = simplify(dfdx.subs(q, q3))
+
+df2dx = simplify(df2dx)
+
+
+# integral = (c1 + c2 + c3) / 3 # God I love Bernstein polynomials
+
+
+
+# integral = s.simplify(integral)
+
+parsimony = len(str(df2dx))
+print(s.pretty(df2dx, num_columns=100))
+print(f"Parsimony: {parsimony}")
\ No newline at end of file
diff --git a/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_derivation_bernstein.py b/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_derivation_bernstein.py
new file mode 100644
index 00000000..c9aab800
--- /dev/null
+++ b/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_derivation_bernstein.py
@@ -0,0 +1,159 @@
+import sympy as s
+from sympy import init_printing
+
+init_printing()
+
+# Gets the integral from x2 to x3, looking at the cubic spline interpolant from x1...x4
+
+# Define the symbols
+hm, h, hp = s.symbols("hm h hp", real=True)
+dfm, df, dfp = s.symbols("dfm df dfp", real=True)
+
+x1, x2, x3, x4 = s.symbols('x1 x2 x3 x4', real=True)
+# f1, f2, f3, f4 = s.symbols('f1 f2 f3 f4', real=True)
+
+f1 = 0
+f2 = dfm
+f3 = f2 + df
+f4 = f3 + dfp
+
+
+def simplify(expr, maxdepth=10, _depth=0):
+    import copy
+    original_expr = copy.copy(expr)
+    print(f"Depth: {_depth} | Parsimony: {len(str(expr))}")
+    maps = {
+        f4 - f3: dfp,
+        f3 - f2: df,
+        f2 - f1: dfm,
+        f4 - f2: dfp + df,
+        f3 - f1: df + dfm,
+        f4 - f1: dfp + df + dfm,
+        x4 - x3: hp,
+        x3 - x2: h,
+        x2 - x1: hm,
+        x4 - x2: hp + h,
+        x3 - x1: h + hm,
+        x4 - x1: hp + h + hm,
+    }
+    expr = expr.factor([h, hm, hp])
+    expr = expr.subs(maps)
+    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
+
+
+q = s.symbols('q')  # Normalized space for a Bernstein basis.
+# Mapping from x-space to q-space has x=x2 -> q=0, x=x3 -> q=1.
+
+q2 = 0
+q3 = 1
+q1 = q2 - hm / h
+q4 = q3 + hp / h
+# q1, q4 = s.symbols('q1 q4', real=True)
+
+# Define the Bernstein basis polynomials
+b1 = (1 - q) ** 3
+b2 = 3 * q * (1 - q) ** 2
+b3 = 3 * q ** 2 * (1 - q)
+b4 = q ** 3
+
+c1, c2, c3, c4 = s.symbols('c1 c2 c3 c4', real=True)
+
+# Can solve for c2 and c3 exactly
+c1 = f2
+c4 = f3
+
+f = c1 * b1 + c2 * b2 + c3 * b3 + c4 * b4
+
+f1_cubic = f.subs(q, q1)  # .simplify()
+f4_cubic = f.subs(q, q4)  # .simplify()
+
+# factors = [q1, q4]
+# factors = [f1, f2, f3, f4]
+
+# Solve for c2 and c3
+sol = s.solve(
+    [
+        f1_cubic - f1,
+        f4_cubic - f4,
+    ],
+    [
+        c2,
+        c3,
+    ],
+)
+c2_sol = simplify(sol[c2].factor([dfm, df, dfp]))
+c3_sol = simplify(sol[c3].factor([dfm, df, dfp]))
+
+# f = c1 * b1 + c2_sol * b2 + c3_sol * b3 + c4 * b4
+
+dfdx = f.diff(q) / h
+df2dx = dfdx.diff(q) / h
+
+res = s.integrate(df2dx ** 2, (q, q2, q3)) * h
+res = simplify(res.factor([dfm, df, dfp]))
+res = simplify(res)
+
+res = simplify(res.subs({c2: c2_sol, c3: c3_sol}).factor([dfm, df, dfp]))
+
+cse = s.cse(
+    res,
+    symbols=s.numbered_symbols("s"),
+    list=False)
+
+parsimony = len(str(res))
+# print(s.pretty(res, num_columns=100))
+print(f"Parsimony: {parsimony}")
+
+for i, (var, expr) in enumerate(cse[0]):
+    print(f"{var} = {expr}")
+print(f"res = {cse[1]}")
+
+if __name__ == '__main__':
+    x = s.symbols('x', real=True)
+
+    a = 0
+    b = 4
+
+    def example_f(x):
+        return x ** 3
+
+    h_val = b - a
+    hm_val = 1
+    hp_val = 3
+    df_val = example_f(b) - example_f(a)
+    dfm_val = example_f(a) - example_f(a - hm_val)
+    dfp_val = example_f(b + hp_val) - example_f(b)
+
+    subs = {
+            h: h_val,
+            hm: hm_val,
+            hp: hp_val,
+            df: df_val,
+            dfm: dfm_val,
+            dfp: dfp_val,
+        }
+
+    exact = s.N(
+        s.integrate(
+            example_f(x).diff(x, 2) ** 2,
+            (x, a, b)
+        )
+    )
+
+    eqn = s.N(
+        res.subs(subs)
+    )
+
+    print(f"exact: {exact}")
+    print(f"eqn: {eqn}")
+    print(f"ratio: {exact/eqn}")
\ No newline at end of file
diff --git a/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_derivation_lagrange.py b/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_derivation_lagrange.py
new file mode 100644
index 00000000..d96532b9
--- /dev/null
+++ b/aerosandbox/numpy/integrate_discrete_derivations/discrete_squared_curvature_derivation_lagrange.py
@@ -0,0 +1,123 @@
+import sympy as s
+from sympy import init_printing
+
+init_printing()
+
+# Gets the integral from x2 to x3, looking at the cubic spline interpolant from x1...x4
+
+# Define the symbols
+hm, h, hp = s.symbols("hm h hp", real=True)
+dfm, df, dfp = s.symbols("dfm df dfp", real=True)
+
+x1, x2, x3, x4 = s.symbols('x1 x2 x3 x4', real=True)
+f1, f2, f3, f4 = s.symbols('f1 f2 f3 f4', real=True)
+
+x = s.symbols('x', real=True)
+
+
+def simplify(expr, maxdepth=10, _depth=0):
+    import copy
+    original_expr = copy.copy(expr)
+    print(f"Depth: {_depth} | Parsimony: {len(str(expr))}")
+    maps = {
+        f4 - f3: dfp,
+        f3 - f2: df,
+        f2 - f1: dfm,
+        f4 - f2: dfp + df,
+        f3 - f1: df + dfm,
+        f4 - f1: dfp + df + dfm,
+        x4 - x3: hp,
+        x3 - x2: h,
+        x2 - x1: hm,
+        x4 - x2: hp + h,
+        x3 - x1: h + hm,
+        x4 - x1: hp + h + hm,
+    }
+    print("hi1")
+    # expr = expr.factor(list(maps.keys()))
+    print("hi2")
+    expr = expr.subs(maps)
+    print("hi3")
+    # expr = expr.simplify()
+    print("hi4")
+    if expr != original_expr:
+        if _depth < maxdepth:
+            expr = simplify(expr, maxdepth=maxdepth, _depth=_depth + 1)
+    return expr
+
+
+f = (
+        f1 * (x - x2) * (x - x3) * (x - x4) / ((x1 - x2) * (x1 - x3) * (x1 - x4)) +
+        f2 * (x - x1) * (x - x3) * (x - x4) / ((x2 - x1) * (x2 - x3) * (x2 - x4)) +
+        f3 * (x - x1) * (x - x2) * (x - x4) / ((x3 - x1) * (x3 - x2) * (x3 - x4)) +
+        f4 * (x - x1) * (x - x2) * (x - x3) / ((x4 - x1) * (x4 - x2) * (x4 - x3))
+)
+f = simplify(f)
+
+dfdx = f.diff(x)
+dfdx = simplify(dfdx)
+
+df2dx = f.diff(x, 2)
+df2dx = simplify(df2dx)
+
+res = s.integrate(df2dx ** 2, (x, x2, x3))
+res = simplify(res.factor([f1, f2, f3, f4]))
+
+parsimony = len(str(res))
+print(s.pretty(res, num_columns=100))
+print(f"Parsimony: {parsimony}")
+
+#
+# q = s.symbols('q')  # Normalized space for a Bernstein basis.
+# # Mapping from x-space to q-space has x=x2 -> q=0, x=x3 -> q=1.
+#
+# q2 = 0
+# q3 = 1
+# # q1 = q2 - hm / h
+# # q4 = q3 + hp / h
+# q1, q4 = s.symbols('q1 q4', real=True)
+#
+# # Define the Bernstein basis polynomials
+# b1 = (1 - q) ** 3
+# b2 = 3 * q * (1 - q) ** 2
+# b3 = 3 * q ** 2 * (1 - q)
+# b4 = q ** 3
+#
+# c1, c2, c3, c4 = s.symbols('c1 c2 c3 c4', real=True)
+#
+# # Can solve for c2 and c3 exactly
+# c1 = f2
+# c4 = f3
+#
+# f = c1 * b1 + c2 * b2 + c3 * b3 + c4 * b4
+#
+# f1_cubic = f.subs(q, q1)#.simplify()
+# f4_cubic = f.subs(q, q4)#.simplify()
+#
+# factors = [q1, q4]
+# # factors = [f1, f2, f3, f4]
+#
+# # Solve for c2 and c3
+# sol = s.solve(
+#     [
+#         f1_cubic - f1,
+#         f4_cubic - f4,
+#     ],
+#     [
+#         c2,
+#         c3,
+#     ],
+# )
+# c2 = sol[c2].factor(factors).simplify()
+# c3 = sol[c3].factor(factors).simplify()
+#
+# f = c1 * b1 + c2 * b2 + c3 * b3 + c4 * b4
+#
+# integral = (c1 + c2 + c3 + c4) / 4 # God I love Bernstein polynomials
+# # integral = s.simplify(integral)
+#
+# integral = integral.factor(factors).simplify()
+#
+# parsimony = len(str(integral))
+# print(s.pretty(integral, num_columns=100))
+# print(f"Parsimony: {parsimony}")