Model the bending of a beam (equation 10.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.

We consider the domain $[0, L]$, with $h=\frac{L}{N}$, $x_k=kh$ with $k=0, \dots, N$.

Let's compute $A_{ij}$ with a similar form as in the previous exercise:

$$A = \alpha \begin{bmatrix} D & E & & & \hphantom{\ddots} \\ E^{\top} & \ddots & \ddots & & \\ & \ddots & \ddots & \ddots & \\ & & \ddots & \ddots & E \\ \phantom{\ddots} & & & E^{\top} & D \\ \end{bmatrix} $$

With $\alpha$ being related to the material's properties, considered uniform:

$$\alpha = E_{\rm Y} I$$

And the submatrices are giben by:

$$D = \begin{bmatrix} <\varphi''_0 | \varphi''_0> & <\varphi''_0 | \varphi''_1>\\ <\varphi''_1 | \varphi''_0> & <\varphi''_1 | \varphi''_1> \end{bmatrix} $$$$E = \begin{bmatrix} <\varphi''_0 | \varphi''_2> & <\varphi''_0 | \varphi''_3>\\ <\varphi''_1 | \varphi''_2> & <\varphi''_1 | \varphi''_3> \end{bmatrix} $$

We can use sympy again to solve those:

import sympy

h = sympy.symbols("h", real=True)
x = sympy.symbols("x", real=True)

p1000 = [1, 0, -3/(h*h), 2/(h*h*h)]
p0100 = [0, 1, -2/h, 1/(h*h)]
p0010 = [0, 0, 3/(h*h), -2/(h*h*h)]
p0001 = [0, 0, -1/h, 1/(h*h)]

x_pow = [1, x, x*x, x*x*x]
poly1000 = sum([a*b for a, b in zip(p1000, x_pow)])
poly0100 = sum([a*b for a, b in zip(p0100, x_pow)])
poly0010 = sum([a*b for a, b in zip(p0010, x_pow)])
poly0001 = sum([a*b for a, b in zip(p0001, x_pow)])

poly1000_deriv = sympy.diff(poly1000, x)
poly0100_deriv = sympy.diff(poly0100, x)
poly0010_deriv = sympy.diff(poly0010, x)
poly0001_deriv = sympy.diff(poly0001, x)

poly1000_deriv2 = sympy.diff(poly1000_deriv, x)
poly0100_deriv2 = sympy.diff(poly0100_deriv, x)
poly0010_deriv2 = sympy.diff(poly0010_deriv, x)
poly0001_deriv2 = sympy.diff(poly0001_deriv, x)

D11 = 2*sympy.integrate(poly1000_deriv2*poly1000_deriv2, (x, 0, h))
D22 = 2*sympy.integrate(poly0100_deriv2*poly0100_deriv2, (x, 0, h))

E11 = sympy.integrate(poly1000_deriv2*poly0010_deriv2, (x, 0, h))
E12 = sympy.integrate(poly1000_deriv2*poly0001_deriv2, (x, 0, h))
E21 = sympy.integrate(poly0100_deriv2*poly0010_deriv2, (x, 0, h))
E22 = sympy.integrate(poly0100_deriv2*poly0001_deriv2, (x, 0, h))

display(sympy.Eq(sympy.Symbol(r"D_{11}"), D11))
display(sympy.Eq(sympy.Symbol(r"D_{22}"), D22))

display(sympy.Eq(sympy.Symbol(r"E_{11}"), E11))
display(sympy.Eq(sympy.Symbol(r"E_{12}"), E12))
display(sympy.Eq(sympy.Symbol(r"E_{21}"), E21))
display(sympy.Eq(sympy.Symbol(r"E_{22}"), E22))

$\displaystyle D_{11} = \frac{24}{h^{3}}$

$\displaystyle D_{22} = \frac{8}{h}$

$\displaystyle E_{11} = - \frac{12}{h^{3}}$

$\displaystyle E_{12} = \frac{6}{h^{2}}$

$\displaystyle E_{21} = - \frac{6}{h^{2}}$

$\displaystyle E_{22} = \frac{2}{h}$

Therefore:

$$D = \begin{bmatrix} \frac{24}{h^3} & 0\\ 0 & \frac{8}{h} \end{bmatrix} $$$$E = \begin{bmatrix} -\frac{12}{h^3} & \frac{6}{h^2}\\ -\frac{6}{h^2} & \frac{2}{h} \end{bmatrix} $$

The vector $b$ is trivial, as the force is only applied on the last element. Considering a perfect impulse centered on the last spline, only one item in $b$ is non-zero:

$$b=\begin{bmatrix} 0\\ \vdots \\ 0\\ F\\ 0 \end{bmatrix} $$

We can simulate the system by inverting $A$. In the following Python code, there is no optimization for inverting $A$ by its diagonals. This is not necessary when the number of elements remains small.

Here are the settings for the simulation:

$$ \begin{cases} L = 1\\ N = 8\\ EI = 1\\ F = 0.1\\ \end{cases} $$

The step size $h$ is inferred as $h=L/N$. When plotting the result, we can have a quick approximation of the solution by taking $a_{2k}$ as the values. For knowing the solution's values in between those points, we use the Hermite spline interpolations.

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (1000/72., 4.8)


def hermite1000(x, h):
    mask = (x >= -h) & (x <= h)
    return mask * np.polyval([2/(h*h*h), -3/(h*h), 0, 1], np.abs(x))


def hermite0100(x, h):
    mask = (x >= -h) & (x <= h)
    return mask * np.polyval([1/(h*h), -2/h, 1, 0], np.abs(x))*np.sign(x)


def solve_beam(h, n, EI, F):
    n_tot = 2*(n+1) # from 0 to n, 2 basis func per point
    A = np.zeros((n_tot, n_tot))
    b = np.zeros((n_tot,))

    for k in range(n+1):
        A[2*k, 2*k] = 24/(h*h*h)
        A[2*k+1, 2*k+1] = 8/h

    for k in range(n):
        A[2*k, 2*k+2] = -12/(h*h*h)
        A[2*k+2, 2*k] = -12/(h*h*h)

        A[2*k, 2*k+3] = 6/(h*h)
        A[2*k+3, 2*k] = 6/(h*h)

        A[2*k+1, 2*k+2] = -6/(h*h)
        A[2*k+2, 2*k+1] = -6/(h*h)

        A[2*k+1, 2*k+3] = 2/h
        A[2*k+3, 2*k+1] = 2/h

    A *= EI

    # force on the last element, centered
    b[-2] = F

    # fixed position and slope
    alpha1 = 0
    alpha2 = 0

    b[0] = alpha1
    b[1] = alpha2

    b[2:] = b[2:] - A[2:, 0] * alpha1 - A[2:, 1] * alpha2
    A[:, 0] = 0
    A[0, :] = 0
    A[:, 1] = 0
    A[1, :] = 0
    A[0, 0] = 1
    A[1, 1] = 1

    return np.dot(np.linalg.inv(A), b)


L = 1
n = 8
n_smooth = 256
EI = 1
F = 0.1

h = L/n

a = solve_beam(h, n, EI, F)

x_interp = np.linspace(0, L, n_smooth)
y_interp = np.zeros_like(x_interp)
for k in range(n+1):
    y_interp += a[2*k] * hermite1000(x_interp-h*k, h)
    y_interp += a[2*k+1] * hermite0100(x_interp-h*k, h)

# rough solution
plt.figure()
plt.subplot(121)
plt.plot(np.linspace(0, L, n+1), a[0::2])
plt.plot([0, L], [0, 0], "k--")
plt.title(r'rough $y(x)$, provided by $a_{2k}$')

# interpolated solution
plt.subplot(122)
plt.plot(x_interp, y_interp)
plt.plot([0, L], [0, 0], "k--")
plt.title('exact $y(x)$, interpolated')

plt.show()

We can modify $F$ and animate the solution:

import matplotlib.animation
import matplotlib
from IPython.display import HTML
import math

n_curve = 256
F_curve = np.sin(np.linspace(0, 2*math.pi, n_curve))
F_i = 0

def update(frame, ln):
    global L, n, n_smooth, EI, F_curve, F_i
    F_i += 1
    F_i %= n_curve
    a = solve_beam(h, n, EI, F_curve[F_i])

    x_interp = np.linspace(0, L, n_smooth)
    y_interp = np.zeros_like(x_interp)
    for k in range(n+1):
        y_interp += a[2*k] * hermite1000(x_interp-h*k, h)
        y_interp += a[2*k+1] * hermite0100(x_interp-h*k, h)
    ln.set_data(x_interp, y_interp)
    
fig, ax = plt.subplots()
y_init = np.zeros_like(x_interp)
y_init[0] = 0.002
y_init[-1] = -y_init[0]
ln = ax.plot(x_interp, y_init)

def init():
    return ln,

ani = matplotlib.animation.FuncAnimation(fig, update, init_func=init, fargs=ln, frames=n_curve, interval=10)

# plt.show()
plt.close()
HTML(ani.to_html5_video())