Add derivations

aerosandbox/numpy/derivative_discrete_derivations/quadratic.py

@@ -0,0 +1,97 @@
+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
+
+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,
+    })
+    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))
+
+
+# integral = (c1 + c2 + c3) / 3 # God I love Bernstein polynomials
+
+
+
+# integral = s.simplify(integral)
+
+parsimony = len(str(dfdx_q1))
+print(s.pretty(dfdx_q1, num_columns=100))
+print(f"Parsimony: {parsimony}")

aerosandbox/numpy/integrate_discrete_derivations/cubic_spline_interpolant_derivation.py

@@ -2,6 +2,8 @@
 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
 x1, x2, x3, x4 = s.symbols('x1 x2 x3 x4', real=True)
 f1, f2, f3, f4 = s.symbols('f1 f2 f3 f4', real=True)

aerosandbox/numpy/integrate_discrete_derivations/simpson_backward_interpolant_derivation.py

@@ -2,7 +2,7 @@
 from sympy import init_printing
 init_printing()
 
-# "Forward" -> gets the integral from x1 to x2, by analogy to forward Euler
+# "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)