Beam Deflection Analysis¶

This notebook analyses beam deflection using the Euler-Bernoulli beam theory. We derive the governing ODE symbolically, obtain the exact closed-form solution for a simply supported beam under uniform load, and then solve the same problem numerically using the finite difference method.

A convergence study validates the numerical approach, and we compute derived quantities (bending moment, shear force) both symbolically and numerically.

In [49]:
load("numerics")$
load("ax-plots")$

Euler-Bernoulli Beam Equation [S+N]¶

The Euler-Bernoulli beam equation governs the deflection $y(x)$ of a slender beam under transverse loading:

$$EI \frac{d^4 y}{dx^4} = w(x)$$

where $EI$ is the flexural rigidity (product of Young's modulus and second moment of area) and $w(x)$ is the distributed load per unit length. For a simply supported beam of length $L$ with uniform load $w_0$, the boundary conditions are:

$$y(0) = 0, \quad y''(0) = 0, \quad y(L) = 0, \quad y''(L) = 0$$

In [50]:
/* Beam equation: EI * y''''(x) = w(x) */
/* Simply supported beam of length L with uniform load w0 */
depends(y, x)$
beam_eqn : EI * diff(y, x, 4) = w0;
print("Governing equation:", beam_eqn)$

EI*'diff(y,x,4) = w0
Governing equation: EI*'diff(y,x,4) = w0

Symbolic Solution for Simple Cases¶

For a uniform load $w(x) = w_0$, the 4th-order ODE can be solved by direct integration four times, introducing four constants of integration $C_1, C_2, C_3, C_4$. Applying the four boundary conditions for a simply supported beam determines all constants uniquely.

In [51]:
/* General solution by integrating 4 times */
y_general : integrate(integrate(integrate(integrate(w0/EI, x), x), x), x)$
/* Add integration constants */
y_general : y_general + C1*x^3/6 + C2*x^2/2 + C3*x + C4$
print("General solution:", y_general)$

/* Apply BCs: y(0) = 0, y''(0) = 0, y(L) = 0, y''(L) = 0 */
eqs : [subst(x=0, y_general) = 0,
       subst(x=0, diff(y_general, x, 2)) = 0,
       subst(x=L, y_general) = 0,
       subst(x=L, diff(y_general, x, 2)) = 0]$
print("Boundary conditions:")$
for eq in eqs do print("  ", eq)$

sol : solve(eqs, [C1, C2, C3, C4])$
print("Integration constants:", sol[1])$

y_exact : subst(sol[1], y_general)$
y_exact : ratsimp(y_exact);
print("Exact deflection y(x) =", y_exact)$

General solution: (w0*x^4)/(24*EI)+(C1*x^3)/6+(C2*x^2)/2+C3*x+C4
Boundary conditions:
   C4 = 0
   C2 = 0
   (L^4*w0)/(24*EI)+(C1*L^3)/6+(C2*L^2)/2+C3*L+C4 = 0
   (L^2*w0)/(2*EI)+C1*L+C2 = 0
Integration constants:
                      [C1 = -((L*w0)/(2*EI)),C2 = 0,C3 = (L^3*w0)/(24*EI),
                       C4 = 0]
(w0*x^4-2*L*w0*x^3+L^3*w0*x)/(24*EI)
Exact deflection y(x) = (w0*x^4-2*L*w0*x^3+L^3*w0*x)/(24*EI)
In [52]:
/* Substitute numeric values and plot the exact deflection */
params : [EI=1, w0=1, L=10]$
y_plot : subst(params, y_exact)$

ax_draw2d(
  color="blue", line_width=2, name="Exact deflection",
  explicit(y_plot, x, 0, 10),
  title="Simply Supported Beam \u2014 Exact Solution",
  xlabel="Position x", ylabel="Deflection y",
  grid=true
)$
No description has been provided for this image

Finite Difference Method¶

We discretise the 4th derivative using the central difference stencil:

$$y''''(x_i) \approx \frac{y_{i-2} - 4y_{i-1} + 6y_i - 4y_{i+1} + y_{i+2}}{h^4}$$

For $n$ interior points on the interval $[0, L]$ with spacing $h = L/(n+1)$, this produces an $n \times n$ pentadiagonal linear system $K\mathbf{y} = \mathbf{f}$. The displacement BCs $y_0 = y_{n+1} = 0$ are absorbed by dropping those columns. The moment BCs $y''(0) = 0$ and $y''(L) = 0$ introduce ghost points $y_{-1}$ and $y_{n+2}$; since $y''(0) \approx (y_{-1} - 2y_0 + y_1)/h^2 = 0$ with $y_0 = 0$, we get $y_{-1} = -y_1$, and similarly $y_{n+2} = -y_n$. Substituting these into the stencil modifies the first and last diagonal entries from $6/h^4$ to $5/h^4$.

In [53]:
/* 4th derivative central difference stencil:              */
/* y''''(x_i) ~ (y_{i-2} - 4*y_{i-1} + 6*y_i              */
/*               - 4*y_{i+1} + y_{i+2}) / h^4              */
/* For n interior points, this gives an (n x n) system     */

/* Display the stencil symbolically */
stencil : (y[i-2] - 4*y[i-1] + 6*y[i] - 4*y[i+1] + y[i+2]) / h^4;
print("FD stencil for y'''':", stencil)$

(y[i+2]-4*y[i+1]+6*y[i]-4*y[i-1]+y[i-2])/h^4
FD stencil for y'''': (y[i+2]-4*y[i+1]+6*y[i]-4*y[i-1]+y[i-2])/h^4
In [54]:
/* Build the pentadiagonal stiffness matrix K */
n : 20$  /* interior points */
h_val : float(subst(params, L) / (n + 1))$
print("Grid spacing h =", h_val)$

K : np_zeros([n, n])$
coeffs : [1, -4, 6, -4, 1]$
for i : 0 thru n - 1 do (
  for k : -2 thru 2 do block([j : i + k],
    if j >= 0 and j < n then
      np_set(K, i, j, np_ref(K, i, j) + coeffs[k + 3] / h_val^4)
  )
)$

/* Apply moment BCs: y''(0)=0 gives y_{-1} = -y_1 (ghost point) */
/* y''(L)=0 gives y_{n+2} = -y_n (ghost point)                  */
/* Each adds -1/h^4 to the corresponding diagonal entry          */
np_set(K, 0, 0, np_ref(K, 0, 0) - 1/h_val^4)$
np_set(K, n-1, n-1, np_ref(K, n-1, n-1) - 1/h_val^4)$

print("Stiffness matrix K: shape", np_shape(K))$

Grid spacing h = 0.47619047619047616
Stiffness matrix K: shape [20,20]
In [55]:
/* Build the RHS force vector and solve */
/* RHS: w0/EI at each interior point */
f : np_full([n, 1], float(subst(params, w0/EI)))$

/* Solve K * y = f */
y_fd : np_solve(K, f)$
print("FD solution computed, shape:", np_shape(y_fd))$
print("Max deflection (FD):", np_max(y_fd))$

/* Compare to exact midpoint deflection */
y_mid_exact : float(subst(append(params, [x=5]), y_exact))$
print("Exact midpoint deflection:", y_mid_exact)$

FD solution computed, shape: [20,1]
Max deflection (FD): 130.08982882645722
Exact midpoint deflection: 130.20833333333334
In [56]:
/* Plot FD solution vs exact solution */
x_fd : makelist(float(i * h_val), i, 1, n)$
y_fd_list : np_to_list(np_col(y_fd, 0))$

ax_draw2d(
  color="blue", line_width=2, name="Exact",
  explicit(subst(params, y_exact), x, 0, 10),
  color="red", marker_size=5, name="FD (n=20)",
  points(x_fd, y_fd_list),
  title="Beam Deflection: Exact vs Finite Difference",
  xlabel="Position x", ylabel="Deflection y",
  grid=true, showlegend=true
)$
No description has been provided for this image

Convergence Study¶

As we refine the mesh (increase $n$), the finite difference error should decrease. For the central difference stencil of the 4th derivative, we expect the truncation error to scale as $O(h^2)$, which we verify by computing the midpoint error for several values of $n$.

In [57]:
/* Convergence study: solve for several mesh sizes */
nn_list : [5, 10, 20, 40, 80]$
err_list : []$
hh_list : []$

for nn in nn_list do block([hh, KK, ff, yy, mid_idx, y_mid_fd, err],
  hh : float(subst(params, L) / (nn + 1)),
  /* Build stiffness matrix */
  KK : np_zeros([nn, nn]),
  for i : 0 thru nn - 1 do (
    for k : -2 thru 2 do block([j : i + k],
      if j >= 0 and j < nn then
        np_set(KK, i, j, np_ref(KK, i, j) + coeffs[k + 3] / hh^4)
    )
  ),
  /* Apply moment BCs: ghost points from y''(0)=0 and y''(L)=0 */
  np_set(KK, 0, 0, np_ref(KK, 0, 0) - 1/hh^4),
  np_set(KK, nn-1, nn-1, np_ref(KK, nn-1, nn-1) - 1/hh^4),
  ff : np_full([nn, 1], float(subst(params, w0/EI))),
  yy : np_solve(KK, ff),
  /* Midpoint index (closest to L/2) */
  mid_idx : round(nn / 2) - 1,
  y_mid_fd : np_ref(yy, mid_idx, 0),
  err : abs(y_mid_fd - y_mid_exact),
  hh_list : endcons(hh, hh_list),
  err_list : endcons(err, err_list),
  print("n =", nn, ", h =", hh, ", midpoint error =", err)
)$

n = 5 , h = 1.6666666666666667 , midpoint error = 14.467592592593562
n = 10 , h = 0.9090909090909091 , midpoint error = 0.43577681396661205
n = 20 , h = 0.47619047619047616 , midpoint error = 0.1185045068764623
n = 40 , h = 0.24390243902439024 , midpoint error = 0.031011187625580305
n = 80 , h = 0.12345679012345678 , midpoint error = 0.007940162648679916
In [58]:
/* Plot convergence (log-log) */
log_h : makelist(log(hh_list[i]) / log(10.0), i, 1, length(hh_list))$
log_err : makelist(log(max(err_list[i], 1e-16)) / log(10.0), i, 1, length(err_list))$

/* Estimate convergence rate from last two points */
rate : (last(log_err) - log_err[length(log_err)-1])
      / (last(log_h) - log_h[length(log_h)-1])$
print("Estimated convergence rate (slope):", rate)$

ax_draw2d(
  color="blue", line_width=2, marker_size=6,
  name="FD error",
  lines(log_h, log_err),
  points(log_h, log_err),
  color="gray", line_width=1, name="O(h^2) reference",
  explicit(2*x + last(log_err) - 2*last(log_h), x,
           log_h[1] - 0.2, last(log_h) + 0.2),
  title="Convergence: log(error) vs log(h)",
  xlabel="log10(h)", ylabel="log10(error)",
  grid=true, showlegend=true
)$

Estimated convergence rate (slope): 2.000969477212739
No description has been provided for this image

The blue line tracks the grey $O(h^2)$ reference closely for $n \ge 10$, confirming second-order convergence. The coarsest mesh ($n = 5$) sits well above the reference because with only 5 interior points the mesh is too coarse for the asymptotic error estimate to apply.

Bending Moment and Shear Force [S+N]¶

From beam theory, the internal bending moment and shear force are related to the deflection by:

$$M(x) = EI \, y''(x), \quad V(x) = EI \, y'''(x)$$

We compute these symbolically from the exact solution and verify against numerical differences of the FD solution.

In [59]:
/* Symbolic bending moment and shear force */
M_exact : ratsimp(subst(params, EI * diff(y_exact, x, 2)));
V_exact : ratsimp(subst(params, EI * diff(y_exact, x, 3)));

print("M(x) =", M_exact)$
print("V(x) =", V_exact)$

/* Plot bending moment and shear force */
ax_draw2d(
  color="blue", line_width=2, name="M(x) — Bending moment",
  explicit(M_exact, x, 0, 10),
  color="red", line_width=2, name="V(x) — Shear force",
  explicit(V_exact, x, 0, 10),
  title="Bending Moment and Shear Force",
  xlabel="Position x", ylabel="Value",
  grid=true, showlegend=true
)$

(x^2-10*x)/2
x-5
M(x) = (x^2-10*x)/2
V(x) = x-5
No description has been provided for this image
In [60]:
/* Compute M numerically from FD solution using second differences */
/* M_i = EI * (y_{i-1} - 2*y_i + y_{i+1}) / h^2 */
/* Use the n=20 solution, including boundary values y_0=0, y_{n+1}=0 */
y_full : append([0.0], y_fd_list, [0.0])$  /* add boundary values */
x_full : makelist(float(i * h_val), i, 0, n + 1)$

M_fd : []$
M_exact_pts : []$
for i : 1 thru n do block([M_i],
  M_i : float(subst(params, EI)) * (y_full[i] - 2*y_full[i+1] + y_full[i+2]) / h_val^2,
  M_fd : endcons(M_i, M_fd),
  M_exact_pts : endcons(float(subst(x = x_fd[i], M_exact)), M_exact_pts)
)$

ax_draw2d(
  color="blue", line_width=2, name="M(x) exact",
  explicit(M_exact, x, 0, 10),
  color="red", marker_size=4, name="M(x) FD",
  points(x_fd, M_fd),
  title="Bending Moment: Exact vs Finite Difference",
  xlabel="Position x", ylabel="M(x)",
  grid=true, showlegend=true
)$
No description has been provided for this image

Summary¶

This notebook demonstrated the complementary use of symbolic and numerical methods for beam analysis:

  1. Symbolic ODE solution — Maxima integrates the 4th-order beam equation directly and applies boundary conditions to obtain the exact deflection in closed form.
  2. Finite difference discretisation — The 4th-derivative central difference stencil produces a pentadiagonal linear system solved numerically with np_solve.
  3. Convergence study — Refining the mesh confirms $O(h^2)$ convergence, validating the numerical method against the exact solution.
  4. Derived quantities — Bending moment $M(x)$ and shear force $V(x)$ are computed both symbolically (by differentiation) and numerically (by finite differences), showing good agreement.

The symbolic-numeric bridge allows engineers to derive governing equations in general form and then instantiate them for specific load cases, material properties, and boundary conditions.