Visualizing Bending Moment Diagrams with Matplotlib for Structural Engineers

title: "Visualizing Bending Moment Diagrams with Matplotlib for Structural Engineers"
published: true
description: "A practical Python tutorial to generate accurate bending moment and shear force diagrams for simple and continuous beams — no FEA software needed."

tags: python, engineering, matplotlib, tutorial

If you've ever worked on a structural project, you know that bending moment diagrams (BMDs) are the backbone of beam design. Every reinforced concrete beam, every steel section check starts here.

The problem? Most structural engineers either rely on black-box software that gives them the diagram without the intuition, or they draw it by hand — which works but doesn't scale.

In this tutorial, we'll build a fully functional BMD and SFD (Shear Force Diagram) generator in pure Python using NumPy and Matplotlib. No FEA software, no license fee, no black box.

By the end, you'll be able to generate diagrams for:

Simply supported beams with point loads
Beams with uniformly distributed loads (UDL)
Continuous beams with multiple spans

💡 This is inspired by the calculation workflows we use at Structalis, a structural engineering consultancy based in Paris, where we deal with reinforced concrete, steel, and timber structures daily under Eurocode standards.

1. Prerequisites

pip install numpy matplotlib

You should be comfortable with basic Python. Structural mechanics knowledge helps but isn't strictly required — we'll cover the theory as we go.

2. Theory in 60 Seconds

For a beam in static equilibrium, at any cross-section at position x:

Shear Force V(x) = algebraic sum of all vertical forces to the left of x
Bending Moment M(x) = algebraic sum of moments of all forces to the left of x

The key relationships:

dV/dx = -q(x)        (distributed load intensity)
dM/dx = V(x)         (shear is the derivative of moment)

Sign convention (Eurocode / French standard):

Positive moment → sagging (tension at the bottom)
Positive shear → left side up, right side down

3. Simply Supported Beam — Point Load

Let's start with the classic case: a simply supported beam of span L, with a single point load P at position a from the left support.

    P
    ↓
A───────────────B
|←──a──→|←─b──→|
←──────L────────→
Ra                Rb

Reactions

Ra = P × b / L
Rb = P × a / L

Shear Force and Bending Moment expressions

For 0 ≤ x < a:
  V(x) = Ra
  M(x) = Ra × x

For a ≤ x ≤ L:
  V(x) = Ra - P
  M(x) = Ra × x - P × (x - a)

Python implementation

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches

def simply_supported_point_load(L, P, a):
    """
    Simply supported beam with a single point load.

    Parameters:
        L (float): Beam span in meters
        P (float): Point load in kN (positive = downward)
        a (float): Distance from left support to load in meters

    Returns:
        x, V, M arrays
    """
    b = L - a

    # Support reactions
    Ra = P * b / L
    Rb = P * a / L

    # Discretize the beam
    x = np.linspace(0, L, 1000)
    V = np.where(x < a, Ra, Ra - P)
    M = np.where(x < a, Ra * x, Ra * x - P * (x - a))

    print(f"Ra = {Ra:.2f} kN  |  Rb = {Rb:.2f} kN")
    print(f"Mmax = {max(M):.2f} kN·m at x = {a:.2f} m")

    return x, V, M


def plot_diagrams(x, V, M, title="Beam Diagrams"):
    """Plot SFD and BMD side by side."""
    fig, axes = plt.subplots(3, 1, figsize=(12, 9), sharex=True)
    fig.suptitle(title, fontsize=14, fontweight='bold')

    L = x[-1]

    # --- Beam schema ---
    ax0 = axes[0]
    ax0.plot([0, L], [0, 0], 'k-', linewidth=5, solid_capstyle='round')
    ax0.plot(0, 0, '^', markersize=15, color='steelblue')
    ax0.plot(L, 0, 'o', markersize=12, color='steelblue')
    ax0.set_xlim(-0.1, L + 0.1)
    ax0.set_ylim(-1, 1)
    ax0.axis('off')
    ax0.set_title('Beam Schema', loc='left', fontsize=10)

    # --- Shear Force Diagram ---
    ax1 = axes[1]
    ax1.fill_between(x, V, 0,
                     where=(V >= 0), color='steelblue', alpha=0.4, label='V > 0')
    ax1.fill_between(x, V, 0,
                     where=(V < 0), color='tomato', alpha=0.4, label='V < 0')
    ax1.plot(x, V, color='steelblue', linewidth=1.5)
    ax1.axhline(0, color='black', linewidth=0.8, linestyle='--')
    ax1.set_ylabel('V (kN)')
    ax1.set_title('Shear Force Diagram', loc='left', fontsize=10)
    ax1.legend(fontsize=8)
    ax1.grid(True, alpha=0.3)

    # --- Bending Moment Diagram ---
    ax2 = axes[2]
    ax2.fill_between(x, M, 0,
                     where=(M >= 0), color='seagreen', alpha=0.4, label='M > 0 (sagging)')
    ax2.fill_between(x, M, 0,
                     where=(M < 0), color='orange', alpha=0.4, label='M < 0 (hogging)')
    ax2.plot(x, M, color='seagreen', linewidth=1.5)
    ax2.axhline(0, color='black', linewidth=0.8, linestyle='--')
    ax2.set_ylabel('M (kN·m)')
    ax2.set_xlabel('x (m)')
    ax2.set_title('Bending Moment Diagram', loc='left', fontsize=10)
    ax2.legend(fontsize=8)
    ax2.grid(True, alpha=0.3)

    # Annotate max moment
    idx_max = np.argmax(np.abs(M))
    ax2.annotate(f'Mmax = {M[idx_max]:.1f} kN·m',
                 xy=(x[idx_max], M[idx_max]),
                 xytext=(x[idx_max] + 0.3, M[idx_max] * 0.8),
                 arrowprops=dict(arrowstyle='->', color='black'),
                 fontsize=9, color='darkgreen')

    plt.tight_layout()
    plt.savefig('bmd_point_load.png', dpi=150, bbox_inches='tight')
    plt.show()


# --- Run it ---
L = 6.0   # meters
P = 50.0  # kN
a = 2.0   # meters from left support

x, V, M = simply_supported_point_load(L, P, a)
plot_diagrams(x, V, M, title=f"Simply Supported Beam — Point Load P={P} kN at a={a} m")

Output:

Ra = 33.33 kN  |  Rb = 16.67 kN
Mmax = 66.67 kN·m at x = 2.00 m

4. Simply Supported Beam — Uniformly Distributed Load (UDL)

Now the more common case in practice: a uniformly distributed load q (kN/m) over the full span.

q (kN/m)
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
A──────────────────B
Ra                 Rb

Reactions and equations:

Ra = Rb = q × L / 2

V(x) = Ra - q × x
M(x) = Ra × x - q × x² / 2

Mmax = q × L² / 8  (at midspan x = L/2)

def simply_supported_udl(L, q):
    """
    Simply supported beam with uniformly distributed load.

    Parameters:
        L (float): Beam span in meters
        q (float): Distributed load in kN/m (positive = downward)

    Returns:
        x, V, M arrays
    """
    Ra = Rb = q * L / 2

    x = np.linspace(0, L, 1000)
    V = Ra - q * x
    M = Ra * x - q * x**2 / 2

    Mmax = q * L**2 / 8
    print(f"Ra = Rb = {Ra:.2f} kN")
    print(f"Mmax = {Mmax:.2f} kN·m at midspan (x = {L/2:.2f} m)")

    return x, V, M


# --- Run it ---
L = 8.0   # meters
q = 15.0  # kN/m

x, V, M = simply_supported_udl(L, q)
plot_diagrams(x, V, M, title=f"Simply Supported Beam — UDL q={q} kN/m, L={L} m")

Output:

Ra = Rb = 60.00 kN
Mmax = 120.00 kN·m at midspan (x = 4.00 m)

This parabolic moment diagram is what you'll see in virtually every floor beam design. At Structalis, this is the baseline we verify before moving to any reinforced concrete or steel section check under Eurocode 2 or 3.

5. Combined Loads — Point Load + UDL

Real structures rarely have a single load type. Let's combine both.

def combined_loads(L, q, P, a):
    """
    Simply supported beam with UDL + point load.

    By superposition (linear elasticity):
    M_total(x) = M_udl(x) + M_point(x)
    """
    b = L - a

    # Reactions
    Ra_udl = q * L / 2
    Ra_P   = P * b / L
    Ra     = Ra_udl + Ra_P

    Rb_udl = q * L / 2
    Rb_P   = P * a / L
    Rb     = Rb_udl + Rb_P

    x = np.linspace(0, L, 1000)

    # Shear (by superposition)
    V_udl = Ra_udl - q * x
    V_P   = np.where(x < a, Ra_P, Ra_P - P)
    V     = V_udl + V_P

    # Moment (by superposition)
    M_udl = Ra_udl * x - q * x**2 / 2
    M_P   = np.where(x < a, Ra_P * x, Ra_P * x - P * (x - a))
    M     = M_udl + M_P

    print(f"Ra = {Ra:.2f} kN  |  Rb = {Rb:.2f} kN")
    print(f"Mmax = {max(M):.2f} kN·m at x ≈ {x[np.argmax(M)]:.2f} m")

    return x, V, M


# --- Run it ---
L = 7.0   # meters
q = 10.0  # kN/m
P = 30.0  # kN
a = 3.0   # meters from left support

x, V, M = combined_loads(L, q, P, a)
plot_diagrams(x, V, M, title=f"Combined Loads — UDL {q} kN/m + Point Load {P} kN at {a} m")

6. Continuous Beam — Two Spans

A continuous beam over three supports is the bread-and-butter of slab and floor beam design. We'll use the Three-Moment Equation (Clapeyron's theorem) for the analytical solution.

For two equal spans with UDL:

M0 = M2 = 0  (simply supported ends)
M1 = -q × L² / 8  (hogging moment at intermediate support)

def two_span_continuous_udl(L1, L2, q):
    """
    Two-span continuous beam under UDL using Three-Moment Equation.

    Parameters:
        L1, L2 (float): Span lengths in meters
        q (float): UDL in kN/m (same on both spans)

    Returns:
        x, V, M arrays (concatenated for both spans)
    """
    # Three-Moment Equation: 2*M1*(L1+L2) = -q/4*(L1^3 + L2^3)
    M1 = -q * (L1**3 + L2**3) / (8 * (L1 + L2))

    # Reactions
    Ra = q * L1 / 2 + M1 / L1
    Rc = q * L2 / 2 + M1 / L2
    Rb = q * (L1 + L2) - Ra - Rc

    print(f"M1 (support B) = {M1:.2f} kN·m (hogging)")
    print(f"Ra = {Ra:.2f} kN  |  Rb = {Rb:.2f} kN  |  Rc = {Rc:.2f} kN")

    # Span 1: A → B
    x1 = np.linspace(0, L1, 500)
    V1 = Ra - q * x1
    M_b_left = 0  # Ma = 0
    M1_val = M1
    # Linear correction for moment at B
    M1_arr = Ra * x1 - q * x1**2 / 2 + (M1_val / L1) * x1

    # Span 2: B → C
    x2 = np.linspace(0, L2, 500)
    V2 = (q * L2 / 2 - M1 / L2) - q * x2
    M2_arr = (q * L2 / 2 - M1 / L2) * x2 - q * x2**2 / 2 + M1 * (1 - x2 / L2)

    # Concatenate
    x_full = np.concatenate([x1, x2 + L1])
    V_full = np.concatenate([V1, V2])
    M_full = np.concatenate([M1_arr, M2_arr])

    return x_full, V_full, M_full


# --- Run it ---
L1 = L2 = 5.0  # meters
q = 20.0        # kN/m

x, V, M = two_span_continuous_udl(L1, L2, q)
plot_diagrams(x, V, M, title=f"Two-Span Continuous Beam — UDL {q} kN/m, L1=L2={L1} m")

Output:

M1 (support B) = -62.50 kN·m (hogging)
Ra = 37.50 kN  |  Rb = 125.00 kN  |  Rc = 37.50 kN

Notice the hogging moment at the intermediate support — this is what drives the top reinforcement in a continuous RC slab.

7. Full Reusable Module

Here's everything packaged as a clean module you can drop into any project:

# beam_diagrams.py

import numpy as np
import matplotlib.pyplot as plt


class SimpleBeam:
    def __init__(self, L: float):
        self.L = L
        self.loads = []
        self._x = np.linspace(0, L, 2000)

    def add_point_load(self, P: float, a: float):
        """Add a point load P (kN) at position a (m) from left support."""
        self.loads.append(('point', P, a))
        return self

    def add_udl(self, q: float, x_start: float = 0, x_end: float = None):
        """Add a UDL q (kN/m) from x_start to x_end."""
        if x_end is None:
            x_end = self.L
        self.loads.append(('udl', q, x_start, x_end))
        return self

    def solve(self):
        """Compute reactions, V(x) and M(x) by superposition."""
        x = self._x
        L = self.L
        V = np.zeros_like(x)
        M = np.zeros_like(x)
        Ra_total = 0
        Rb_total = 0

        for load in self.loads:
            if load[0] == 'point':
                _, P, a = load
                b = L - a
                Ra = P * b / L
                Rb = P * a / L
                Ra_total += Ra
                Rb_total += Rb
                V_i = np.where(x < a, Ra, Ra - P)
                M_i = np.where(x < a, Ra * x, Ra * x - P * (x - a))

            elif load[0] == 'udl':
                _, q, xs, xe = load
                length = xe - xs
                resultant = q * length
                centroid = xs + length / 2
                Ra = resultant * (L - centroid) / L
                Rb = resultant * centroid / L
                Ra_total += Ra
                Rb_total += Rb

                def udl_contribution(xi, q=q, xs=xs, xe=xe, Ra=Ra):
                    loaded = np.clip(xi - xs, 0, xe - xs)
                    V_i = Ra - q * loaded
                    M_i = Ra * (xi - xs) * (xi >= xs) - q * loaded**2 / 2
                    return V_i, M_i

                V_i, M_i = udl_contribution(x)

            V += V_i
            M += M_i

        self.Ra = Ra_total
        self.Rb = Rb_total
        self.x  = x
        self.V  = V
        self.M  = M

        print(f"Ra = {Ra_total:.2f} kN  |  Rb = {Rb_total:.2f} kN")
        print(f"Mmax = {max(M):.2f} kN·m at x = {x[np.argmax(M)]:.2f} m")
        return self

    def plot(self, title="Beam Analysis"):
        plot_diagrams(self.x, self.V, self.M, title=title)
        return self


# Usage example
if __name__ == "__main__":
    beam = SimpleBeam(L=8.0)
    beam.add_udl(q=12.0)          # 12 kN/m full span
    beam.add_point_load(P=40.0, a=3.0)  # 40 kN at 3 m
    beam.solve()
    beam.plot(title="Custom Beam — UDL 12 kN/m + Point Load 40 kN at 3 m")

8. Going Further

This tutorial covers linear elastic static analysis — the foundation of any structural check. From here, you can:

Plug into Eurocode checks : use Mmax to verify the required steel area As = Med / (fyd × z) for reinforced concrete under EC2
Connect to OpenSeesPy for non-linear analysis and seismic loads
Export to PDF with ReportLab to generate automated calculation notes
Read from IFC with IfcOpenShell to extract beam geometry directly from your Revit model

If you're working on a real structural project and need expert review of your beam design — from preliminary sizing to full Eurocode justification — Structalis provides structural engineering services across France, covering reinforced concrete, steel, and timber structures.

Conclusion

In just ~100 lines of Python, we've built a tool that:

✅ Computes support reactions by superposition

✅ Generates accurate V(x) and M(x) diagrams

✅ Handles point loads, UDL, and combined loads

✅ Solves continuous beams with the Three-Moment Equation

✅ Packages everything into a reusable SimpleBeam class

The code is intentionally transparent — no black box, every formula is visible and traceable. That's exactly how structural calculation notes should work.

Found this useful? Drop a ❤️ and share with an engineer friend. Questions or improvements? Open a discussion below — I read every comment.

References:

EN 1990: Eurocode — Basis of structural design
EN 1992-1-1: Eurocode 2 — Design of concrete structures
Timoshenko, S. & Gere, J. — Mechanics of Materials
Structalis — Bureau d'études structure

推荐订阅源

DEV Community