Use better basis functions

_code/pset_06/scripts/beam2.py

+import math
+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)
+# used to specify the applied force
+f = Symbol("f", 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.i_left = integrate(self.left, x)
+        self.i_right = integrate(self.right, x)
+        self.cumulative_i_right = self.i_left.subs(x, h) + self.i_right
+        self.ii_left = integrate(self.i_left, x)
+        self.ii_right = integrate(self.i_right, x)
+        self.cumulative_ii_right = self.ii_left.subs(x, h) + integrate(self.cumulative_i_right, x)
+
+    def left_value(self, xval):
+        return self.left.subs(x, xval)
+
+    def right_value(self, xval):
+        return self.right.subs(x, xval)
+
+    def integral(self, index, n_nodes):
+        if index == 0:
+            return Piecewise(
+                (0, x < 0),
+                (self.i_right, x < h),
+                (self.i_right.subs(x, h), True)
+            )
+        elif index == n_nodes - 1:
+            return Piecewise(
+                (0, x < (index - 1) * h),
+                (self.i_left.subs(x, x - (index - 1) * h), x < index * h),
+                (self.i_left.subs(x, h), True)
+            )
+        return Piecewise(
+            (0, x < (index - 1) * h),
+            (self.i_left.subs(x, x - (index - 1) * h), x < index * h),
+            (self.cumulative_i_right.subs(x, x - index * h), x < (index + 1) * h),
+            (self.cumulative_i_right.subs(x, h), True)
+        )
+
+    def double_integral(self, index, n_nodes):
+        if index == 0:
+            return Piecewise(
+                (0, x < 0),
+                (self.ii_right, x < h),
+                (self.ii_right.subs(x, h) + self.i_right.subs(x, h) * (x - h), True)
+            )
+        elif index == n_nodes - 1:
+            return Piecewise(
+                (0, x < (index - 1) * h),
+                (self.ii_left.subs(x, x - (index - 1) * h), x < index * h),
+                (self.ii_left.subs(x, h) + self.i_left.subs(x, h) * (x - index * h), True)
+            )
+        return Piecewise(
+            (0, x < (index - 1) * h),
+            (self.ii_left.subs(x, x - (index - 1) * h), x < index * h),
+            (self.cumulative_ii_right.subs(x, x - index * h), x < (index + 1) * h),
+            (self.cumulative_ii_right.subs(x, h) + self.cumulative_i_right.subs(x, h) * (x - (index + 1) * h), True)
+        )
+
+
+def hat_basis_function():
+    p_1 = x / h
+    p_2 = 1 - x / h
+    print(latex(p_1))
+    print(latex(p_2))
+    print("")
+
+    return LocalBasisFunction(p_1, p_2)
+
+#def integrate_curvature(basis_function, coefficients, s0):
+
+
+# Computes the matrix A_{i, j} = integrate(f1_i * f2_j, (x, 0, h)). Both arguments must be
+# LocalBasisFunctions. Note that A is tridiagonal.
+def compute_A(f1, f2):
+    result = {}
+
+    # Compute the diagonal elements. The first and last are special cases since the basis function
+    # is chopped in half.
+    first_diag = integrate(f1.right * f2.right, (x, 0, h))
+    last_diag = integrate(f1.left * f2.left, (x, 0, h))
+    diag = first_diag + last_diag
+    result[0, 0] = first_diag
+    result[k, k] = diag
+    result[n, n] = last_diag
+
+    # Compute off diagonal elements. These entries will always be the same if f1 and f2 are the
+    # same. But they need not be in general.
+    above_diag = integrate(f1.right * f2.left, (x, 0, h))
+    below_diag = integrate(f1.left * f2.right, (x, 0, h))
+    result[k, k + 1] = above_diag
+    result[k, k - 1] = below_diag
+
+    return result
+
+def compute_b(basis_function, force, s0, n_nodes, f_start, length):
+    do_integral = lambda i: integrate(basis_function.double_integral(i, n_nodes), (x, f_start, length))
+    common_term = s0 * force * integrate(x, (x, f_start, length))
+    result = [ force * do_integral(i) + common_term for i in range(0, n_nodes) ]
+    return result
+
+def make_block(f1, f2, n_nodes):
+    block = compute_A(f1, f2)
+    block = render_matrix(block, n_nodes)
+    if n_nodes < 8:
+        print(latex(block))
+        print("")
+    return block
+
+def render_matrix(entries, size):
+    result = zeros(size, size)
+    for i in range(0, size):
+        result[i, i] = entries[(k, k)]
+    for i in range(0, size - 1):
+        result[i, i + 1] = entries[(k, k + 1)]
+    for i in range(1, size):
+        result[i, i - 1] = entries[(k, k - 1)]
+    if (0, 0) in entries:
+        result[0, 0] = entries[(0, 0)]
+    if (n, n) in entries:
+        result[size - 1, size - 1] = entries[(n, n)]
+    if (0, 1) in entries:
+        result[0, 1] = entries[(0, 1)]
+    if (1, 0) in entries:
+        result[1, 0] = entries[(1, 0)]
+    return result
+
+# basis_function is a LocalBasisFunction objects
+# coefficients is a list of coefficients
+def plot_result(name, basis_function, coefficients, s0):
+    n_nodes = len(coefficients)
+    n_subdivs = 10
+
+    integral_sum = s0
+    double_integral_sum = s0 * x
+    for i in range(0, n_nodes):
+        integral_sum += coefficients[i] * basis_function.integral(i, n_nodes)
+        double_integral_sum += coefficients[i] * basis_function.double_integral(i, n_nodes)
+
+    xvals = []
+    yvals = []
+    iyvals = []
+    iiyvals = []
+    for xval in np.linspace(0, n_nodes - 1, n_subdivs * (n_nodes - 1) + 1, endpoint=True):
+        i = math.floor(xval)
+        if i > n_nodes - 2:
+            i = n_nodes - 2
+        y = coefficients[i] * basis_function.right_value(xval - i).subs(h, 1)
+        y += coefficients[i + 1] * basis_function.left_value(xval - i).subs(h, 1)
+        iy  = integral_sum.subs([(x, xval), (h, 1)])
+        iiy = double_integral_sum.subs([(x, xval), (h, 1)])
+        xvals.append(xval)
+        yvals.append(y)
+        iyvals.append(iy)
+        iiyvals.append(iiy)
+
+    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, iiyvals, label="displacement")
+    ax1.plot(xvals, iyvals, label="slope")
+    ax1.plot(xvals, yvals, label="curvature")
+    ax1.set_xlabel('h')
+    ax1.set_title(name)
+    ax1.legend()
+    #plt.show(fig1)
+    fig1.savefig(f"../../../assets/img/06_beam2_{name}.png", transparent=True)
+
+if __name__ == "__main__":
+    # define simulation parameters
+    # these are the ones I modified to produce different charts
+    n_nodes = 5
+    slope_0 = 0.5
+    force = -0.5
+    f_start = 0.666
+
+    # these I left constant
+    length = 1.0
+    # n nodes means n - 1 intervals
+    # there are 2 * n_intervals degrees of freedom (we fix the coefficients for the first node)
+    n_intervals = n_nodes - 1
+    hval = length / n_intervals
+    ei = 1.0
+
+    print(f"{n_nodes} nodes")
+    print(f"{n_intervals} elements")
+    print(f"h = {hval}")
+    print(f"E*I = {ei}")
+    print(f"force = {force}")
+    print(f"force range = ({f_start}, {length}")
+    print(f"fixed slope = {slope_0}")
+    print("")
+
+    # construct basis functions
+    basis_function = hat_basis_function()
+    plot_result("curvature_hats", basis_function, [0, 1, 0, 0], 0)
+
+    # construct A
+    A = make_block(basis_function, basis_function, n_nodes)
+    A = A.subs(h, hval)
+    A = np.array(A).astype(np.float64)
+    if n_nodes < 8:
+        print(A)
+        print("")
+
+    # compute b
+    b = compute_b(basis_function, f, slope_0, n_nodes, f_start, length)
+    #if n_nodes < 8:
+    #    print(latex(b))
+    b = [ expr.subs([(h, hval), (f, force)]) for expr in b ]
+    b = np.array(b).astype(np.float64)
+    if n_nodes < 8:
+        print(latex(b))
+        print("")
+
+    a = np.linalg.solve(A, b)
+    if n_nodes < 8:
+        print(a)
+        print("")
+
+    plot_result(f"{n_nodes}_nodes", basis_function, a, slope_0)

_psets/06.md

@@ -493,7 +493,7 @@ solution.
 
 ![psi](../assets/img/06_5_nodes.png)
 
-It's unexpectedly wobbly. One way to get rid of this is to use only two nodes. But this kind of
+It's unexpectedly wobbly. One way to get rid of this is to use only two nodes. But this kinda
 defeats the purpose of FEM.
 
 ![psi](../assets/img/06_2_nodes.png)
@@ -501,3 +501,57 @@ defeats the purpose of FEM.
 Going to more nodes makes things worse.
 
 ![psi](../assets/img/06_100_nodes.png)
+
+Ultimately I think the problem is our choice of basis functions. $$\zeta$$ has a massive
+discontinuity in its second derivative, which is exactly what matters for the potential energy of
+the beam. So let's try something different.
+
+In particular, since the second derivative is what's most important, let's model it directly.
+
+$$
+\frac{d^2 u(x)}{dx^2} = \sum_{i=1}^n a_i \varphi_i(x)
+$$
+
+We can find $$u(x)$$ by integrating. This produces an expression that's no longer a sum of local
+basis functions, but that doesn't matter for this problem.
+
+$$
+\begin{aligned}
+\frac{d u(x)}{dx} &= \sum_{i=1}^n a_i \int_0^x \varphi_i(x') dx' + s_0 \\
+u(x) &= \sum_{i=1}^n a_i \int_0^x \int_0^{x'} \varphi_i(x'') dx'' dx' + s_0 x + u_0 \\
+\end{aligned}
+$$
+
+Here $$s_0$$ is a constant of integration that determines the slope at $$x = 0$$, and $$u_0$$ is a
+constant of integration that determines the displacement at $$x = 0$$. The problem is invariant
+under translation up and down, so without loss of generality I'll assume $$u_0 = 0$$. (If we want to
+change this later, we just need to add the new $$u_0$$ to all the displacements we calculate for
+$$u_0 = 0$$.) Here are what these basis functions look like.
+
+![psi](../assets/img/06_beam2_curvature_hats.png)
+
+Plugging this in,
+
+$$
+\begin{aligned}
+V = &\int_0^L \frac{1}{2} E I \left( \sum_{j=1}^n a_j \varphi_j(x) \right)^2 dx \\
+&-\int_0^L \left( \sum_{j=1}^n a_j \int_0^x \int_0^{x'} \varphi_j(x'') dx'' dx' + s_0 x \right) f(x) dx
+\end{aligned}
+$$
+
+So
+
+$$
+\begin{aligned}
+\frac{\partial V}{\partial a_i}
+= &\int_0^L E I \left( \sum_{j=1}^n a_j \varphi_j(x) \right) \varphi_i(x) dx \\
+&-\int_0^L \left( \int_0^x \int_0^{x'} \varphi_i(x'') dx'' dx' + s_0 x \right) f(x) dx \\
+= &\sum_{j=1}^n a_j E I \int_0^L \varphi_j(x) \varphi_i(x) dx \\
+&-\int_0^L \left( \int_0^x \int_0^{x'} \varphi_i(x'') dx'' dx' \right) f(x) dx \\
+&-s_0 \int_0^L x f(x) dx \\
+\end{aligned}
+$$
+
+These results are much better.
+
+![psi](../assets/img/06_beam2_5_nodes.png)