Derive some awful basis functions

My approach was fundamentally flawed since the basis functions aren't
local. That is, I was imagining taking a piecewise linear approximation
of curvature and integrating to get slope and displacement. But I'd have
to integrate from the beginning of the interval to do so.

_code/pset_06/scripts/beam.py

+import numpy as np
+from sympy import *
+from sympy.assumptions.assume import global_assumptions
+import matplotlib.pyplot as plt
+
+# used as a generic variable, i.e. for polynomials, integration, etc.
+x = Symbol("x", real=True)
+# defines the width of a quantized segment
+h = Symbol("h", real=True, positive=True)
+# used to refer to the number of segments (so there are n + 1 points)
+n = Symbol("n", real=True, positive=True)
+# used as a generic index for an interior segment (i.e. not touching the boundary)
+k = Symbol("k", real=True, positive=True)
+
+# Assumes both functions passed in are nonzero only on [0, h].
+class LocalBasisFunction:
+    def __init__(self, left, right):
+        self.left = left
+        self.right = right
+        self.d_left = diff(self.left, x)
+        self.d_right = diff(self.right, x)
+        self.dd_left = diff(self.d_left, x)
+        self.dd_right = diff(self.d_right, x)
+
+    def left_value(self, xval):
+        return self.left.subs(x, xval)
+
+    def left_slope(self, xval):
+        return self.d_left.subs(x, xval)
+
+    def left_curvature(self, xval):
+        return self.dd_left.subs(x, xval)
+
+    def right_value(self, xval):
+        return self.right.subs(x, xval)
+
+    def right_slope(self, xval):
+        return self.d_right.subs(x, xval)
+
+    def right_curvature(self, xval):
+        return self.dd_right.subs(x, xval)
+
+def derive_basis_function():
+    poly_matrix = Matrix([
+        [1, 0, 0,    0   ],
+        [1, h, h**2, h**3],
+        [0, 0, 2,    0   ],
+        [0, 0, 2,    6*h ]
+    ])
+    poly_matrix_inverse = poly_matrix**-1
+    print(latex(poly_matrix_inverse))
+    print("")
+
+    c_1 = simplify(poly_matrix_inverse * Matrix([0, -1, 0, 1]))
+    c_2 = simplify(poly_matrix_inverse * Matrix([-1, 0, 1, 0]))
+    p_1 = c_1.dot(Matrix([1, x, x**2, x**3]))
+    p_2 = c_2.dot(Matrix([1, x, x**2, x**3]))
+    print("coeffs")
+    print(latex(c_1))
+    print(latex(c_2))
+    print("polys")
+    print(latex(p_1))
+    print(latex(p_2))
+    print("poly derivs")
+    print(latex(diff(p_1, x)))
+    print(latex(diff(p_2, x)))
+    print("poly 2nd derivs")
+    print(latex(diff(p_1, x, 2)))
+    print(latex(diff(p_2, x, 2)))
+    print("")
+    print("sum")
+    print(latex(simplify(p_1 + p_2)))
+
+    return LocalBasisFunction(p_1, p_2)
+
+def plot_result(basis_function, coefficients):
+    n_subdivs = 10
+    xvals = []
+    yvals = []
+    dyvals = []
+    ddyvals = []
+    for i in range(0, len(coefficients) - 1):
+        for xval in np.linspace(0, 1, n_subdivs, endpoint=False):
+            y = coefficients[i] * basis_function.right_value(xval).subs(h, 1)
+            y += coefficients[i + 1] * basis_function.left_value(xval).subs(h, 1)
+            dy = coefficients[i] * basis_function.right_slope(xval).subs(h, 1)
+            dy += coefficients[i + 1] * basis_function.left_slope(xval).subs(h, 1)
+            ddy = coefficients[i] * basis_function.right_curvature(xval).subs(h, 1)
+            ddy += coefficients[i + 1] * basis_function.left_curvature(xval).subs(h, 1)
+            xvals.append(i + xval)
+            yvals.append(y)
+            dyvals.append(dy)
+            ddyvals.append(ddy)
+
+    fig1 = plt.figure()
+    left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
+    ax1 = fig1.add_axes([left, bottom, width, height])
+    ax1.plot(xvals, yvals, label="value")
+    ax1.plot(xvals, dyvals, label="slope")
+    ax1.plot(xvals, ddyvals, label="curvature")
+    ax1.set_xlabel('h')
+    ax1.set_ylabel("displacement")
+    ax1.legend()
+    plt.show(fig1)
+    fig1.savefig("../../../assets/img/06_bad_basis_functions.png", transparent=True)
+
+if __name__ == "__main__":
+    basis_function = derive_basis_function()
+    plot_result(basis_function, [0.1, 0.1, 0.1, 0.1, 0.1])

_psets/06.md

@@ -340,3 +340,74 @@ $$\psi_0$$ and $$\psi_n$$ from the set of basis functions.
 Model the bending of a beam (equation 9.29) under an applied load. Use Hermite polynomial
 interpolation, and boundary conditions fixing the displacement and slope at one end, and applying a
 force at the other end.
+
+Now let's come up with a set of basis functions. We need their second derivatives, so they need to
+be quadratic or higher order. The whole problem is about curvature, so I'd rather have a piecewise
+linear approximation of the second derivative than a piecewise constant one. So I'll use cubic basis
+functions.
+
+But which ones? Essentially I want hat functions for the second derivative, and then I can integrate
+twice to figure out what the basis functions are. Of course I'll need some values to determine the
+constants of integration. In the end this boils down to redoing the Hermite interpolation derivation
+using the function values and second derivatives, instead of function values and first derivatives.
+
+So, for a polynomial
+
+$$
+u(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3
+$$
+
+the values $$u(0)$$, $$u(h)$$, $$\ddot{u}(0)$$, and $$\ddot{u}(h)$$ are given by
+
+$$
+\begin{bmatrix} u(0) \\ u(h) \\ \ddot{u}(0) \\ \ddot{u}(h) \end{bmatrix} =
+\begin{bmatrix}
+1 & 0 & 0 & 0 \\
+1 & h & h^{2} & h^{3} \\
+0 & 0 & 2 & 0 \\
+0 & 0 & 2 & 6 h
+\end{bmatrix}
+\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{bmatrix}
+$$
+
+Inverting this matrix, we find that
+
+$$
+\begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{bmatrix} =
+\begin{bmatrix}
+1 & 0 & 0 & 0 \\
+-\frac{1}{h} & \frac{1}{h} & - \frac{h}{3} & - \frac{h}{6} \\
+0 & 0 & \frac{1}{2} & 0 \\
+0 & 0 & - \frac{1}{6 h} & \frac{1}{6 h}
+\end{bmatrix}
+\begin{bmatrix} u(0) \\ u(h) \\ \ddot{u}(0) \\ \ddot{u}(h) \end{bmatrix}
+$$
+
+I want shape functions corresponding to $$(0, -1, 0, 1)$$ and $$(-1, 0, 1, 0)$$. The resulting
+polynomials are
+
+$$
+\begin{aligned}
+p_1(x) &= -\left(\frac{1}{h} + \frac{h}{6} \right) x + \frac{1}{6 h} x^3 \\
+p_2(x) &= -1 + \left( \frac{1}{h} - \frac{h}{3} \right) x + \frac{1}{2} x^2 - \frac{1}{6h} x^3
+\end{aligned}
+$$
+
+The second derivatives of these functions are
+
+$$
+\begin{aligned}
+\ddot{p_1}(x) &= \frac{x}{h} \\
+\ddot{p_2}(x) &= 1 - \frac{x}{h}
+\end{aligned}
+$$
+
+So if we use $$p_1$$ as the left half of a basis function, and $$p_2$$ as the right, each basis
+function determines the curvature and value at one node.
+
+$$
+\frac{6 \left(-1 + \frac{2 x}{h}\right)}{h^{2}} \\
+\frac{2 \left(-2 + \frac{3 x}{h}\right)}{h} \\
+\frac{6 \left(1 - \frac{2 x}{h}\right)}{h^{2}} \\
+\frac{2 \left(-1 + \frac{3 x}{h}\right)}{h}
+$$