Structural Analysis of a Beam with Python: From Equations to FEM

Some links on this page are affiliate links: if you buy through them we may earn a commission, at no extra cost to you.

Python can calculate beam reactions, shear force, bending moment, deflection, and elastic bending stress—provided the beam model, units, supports, loads, and assumptions are defined correctly. For a single slender beam, NumPy and Matplotlib are often enough. SymPy adds symbolic equations, while PyNite or OpenSeesPy becomes useful for larger, three-dimensional, nonlinear, or dynamic models.

This guide starts with a fully checkable simply supported beam example, then shows how the same problem extends into matrix stiffness analysis and finite-element modeling.

#	Preview	Product	Price
1		Structural Analysis	$171.56	Buy on Amazon
2		Structural Analysis: Understanding Behavior	$76.21	Buy on Amazon
3		Fundamentals of Structural Analysis	$56.12	Buy on Amazon
4		Structural Analysis	$19.31	Buy on Amazon
5		Structural Analysis Made Easy: A practice book for calculating statically determined systems	$27.36	Buy on Amazon

As an Amazon Associate I earn from qualifying purchases.

What “beam analysis” means

Beam analysis can refer to several different levels of work:

1. Strength-of-materials analysis: reactions, shear force, bending moment, slope, deflection, and idealized stress.
2. Matrix structural analysis: several beam elements assembled into a global stiffness system.
3. Finite-element analysis: meshes, element formulations, nonlinear behavior, stability, dynamics, and more complex geometry.
4. Design verification: code-specific load combinations, strength, serviceability, buckling, connections, detailing, and professional review.

The Python example below is primarily a linear elastic strength-of-materials calculation. It is not a complete structural design check.

Popular on Amazon

Structural Analysis

$171.56

Buy on Amazon

Model assumptions

The worked example uses the Euler–Bernoulli beam model. It assumes:

• Static loading and small deflection.
• Linear elastic material behavior.
• Constant Young’s modulus, E, and second moment of area, I.
• A slender, prismatic beam.
• Plane sections remain plane and shear deformation is negligible.
• No cracking, yielding, local buckling, contact, large rotations, or material nonlinearity.
• Supports have been idealized correctly.
• All calculations use a consistent unit system.

The governing differential equation is:

EI d⁴v/dx⁴ = q(x)

Here, v(x) is transverse deflection and q(x) is the distributed load. Euler–Bernoulli theory may be inadequate for short or deep beams, sandwich construction, thick composites, or any case where shear deformation is significant. Such problems may require Timoshenko beam theory. Dynamic analysis also requires mass and damping assumptions.

Section properties and units

At minimum, a beam model needs its length, cross-sectional area, second moment of area, Young’s modulus, load values and positions, and support restraints.

For a rectangular section:

A = bh

I = bh³/12

b is the width and h is the bending depth—the dimension cubed. Increasing depth can therefore improve bending stiffness much more effectively than increasing width, where other constraints are unchanged.

Use N, m, and Pa internally. Convert only for display:

• kN for force,
• kN·m for moment,
• mm for deflection,
• MPa for stress.

Mixing kN with N, mm with m, GPa with Pa, or mm⁴ with m⁴ is one of the most common causes of incorrect results.

Worked example: simply supported beam with a uniform load

Consider a simply supported beam carrying a uniform load over its full span:

Also popular on Amazon

Sale

Structural Analysis: Understanding Behavior

$76.21

Buy on Amazon

• Span: L = 6 m
• Uniform load: w = 10 kN/m = 10,000 N/m
• Young’s modulus: E = 200 GPa = 200 × 10⁹ Pa
• Rectangular section: b = 0.15 m, h = 0.30 m

The section properties are:

A = 0.15 × 0.30 = 0.045 m²

I = 0.15 × 0.30³ / 12 = 0.0003375 m⁴

The reactions are equal by symmetry:

R_A = R_B = wL/2 = 30 kN

Taking upward reaction as positive, the shear and bending moment functions are:

V(x) = R_A − wx

M(x) = R_Ax − wx²/2

The maximum moment occurs at midspan, where shear changes sign:

M_max = wL²/8 = 45 kN·m

The maximum elastic deflection is:

v_max = 5wL⁴/(384EI) ≈ 0.0119 m = 11.9 mm

At the extreme fiber, with y = h/2:

σ = My/I

At maximum moment, the idealized elastic bending stress is approximately 20.0 MPa. This is a verification value for the stated model, not a conclusion that the beam is safe. A design decision requires material properties, factored and service load combinations, allowable limits, stability checks, connections, applicable codes, and appropriate review.

Complete NumPy and Matplotlib calculation

The following script evaluates the closed-form equations at 501 positions along the beam and plots the response. It is not yet a finite-element solver: it samples known analytical equations.

import numpy as np
import matplotlib.pyplot as plt

# Geometry, material, and load
L = 6.0                 # beam length, m
w = 10_000.0            # UDL, N/m
E = 200e9               # Young's modulus, Pa
b = 0.15                # rectangular width, m
h = 0.30                # rectangular depth, m

A = b * h
I = b * h**3 / 12
y_max = h / 2

# Static response
RA = w * L / 2
RB = w * L / 2

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

V = RA - w * x
M = RA * x - w * x**2 / 2

# Downward deflection magnitude
v = w * x * (L**3 - 2 * L * x**2 + x**3) / (24 * E * I)

# Elastic bending stress at the extreme fiber
sigma = M * y_max / I

# Checks
M_max = np.max(M)
x_M_max = x[np.argmax(M)]
v_max = np.max(v)
x_v_max = x[np.argmax(v)]
sigma_max = np.max(np.abs(sigma))

print(f"RA = {RA / 1000:.3f} kN")
print(f"RB = {RB / 1000:.3f} kN")
print(f"Maximum moment = {M_max / 1000:.3f} kN·m at x = {x_M_max:.3f} m")
print(f"Maximum deflection = {v_max * 1000:.3f} mm at x = {x_v_max:.3f} m")
print(f"Maximum elastic bending stress = {sigma_max / 1e6:.3f} MPa")

# Closed-form verification
M_expected = w * L**2 / 8
v_expected = 5 * w * L**4 / (384 * E * I)

print(f"Closed-form maximum moment = {M_expected / 1000:.3f} kN·m")
print(f"Closed-form maximum deflection = {v_expected * 1000:.3f} mm")

# Plots
fig, axes = plt.subplots(3, 1, figsize=(9, 10), sharex=True)

axes[0].plot(x, V / 1000, color="tab:blue")
axes[0].axhline(0, color="black", linewidth=0.8)
axes[0].set_ylabel("V (kN)")
axes[0].set_title("Shear-force diagram")
axes[0].grid(True, alpha=0.3)

axes[1].plot(x, M / 1000, color="tab:red")
axes[1].axhline(0, color="black", linewidth=0.8)
axes[1].set_ylabel("M (kN·m)")
axes[1].set_title("Bending-moment diagram")
axes[1].grid(True, alpha=0.3)

axes[2].plot(x, v * 1000, color="tab:green")
axes[2].set_xlabel("Position x (m)")
axes[2].set_ylabel("Deflection (mm)")
axes[2].set_title("Elastic deflection")
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

The expected results are approximately:

Independent checks

Never rely on a plot that merely looks plausible. Check the following:

Equilibrium

The reactions must balance the applied load:

R_A + R_B = wL

Taking moments about either support must also produce zero net moment.

Boundary conditions

For ideal simple supports:

v(0) = 0 and v(L) = 0

With no applied end moments, the bending moment is also zero at both supports.

Also popular on Amazon

Sale

Fundamentals of Structural Analysis

$56.12

Buy on Amazon

Symmetry

A symmetric beam, section, and load should have equal reactions, with maximum moment and deflection at midspan. Shear should change sign there.

Scaling behavior

• Doubling w doubles force, moment, and deflection in a linear model.
• Doubling E halves deflection.
• Doubling I halves deflection.
• Reversing the load reverses the response signs.
• Zero load produces zero response.

Other standard load cases

Once the workflow is understood, the same checks can be applied to other classic cases:

For off-centre point loads, partial distributed loads, internal hinges, springs, or multiple spans, piecewise equations or a stiffness method usually becomes more practical than manually extending one formula.

Symbolic beam analysis with SymPy

SymPy’s Beam class supports beam loads, supports, reactions, boundary conditions, shear, bending moment, slope, and deflection. It is useful when the equations themselves matter—for example, when deriving a piecewise result or retaining symbolic parameters.

import sympy as sp
from sympy.physics.continuum_mechanics.beam import Beam

L = 6
E = 200e9
I = 0.0003375
w = 10_000

beam = Beam(L, E, I)
R1, R2 = sp.symbols("R1 R2")

beam.apply_load(R1, 0, -1)
beam.apply_load(w, 0, 0, end=L)
beam.apply_load(R2, L, -1)

beam.bc_deflection = [(0, 0), (L, 0)]
beam.solve_for_reaction_loads(R1, R2)

print(beam.reaction_loads)
print(beam.shear_force())
print(beam.bending_moment())
print(beam.slope())
print(beam.deflection())

Check the SymPy documentation for the Beam API and its sign convention. Force, shear, moment, slope, and deflection signs depend on the convention you choose and apply consistently; no single diagram orientation is universal.

From formulas to the direct stiffness method

A general beam solver can be built by dividing the member into elements and solving the global system:

Ku = F

1. Divide the beam into elements.
2. Assign two nodes to every element.
3. Give each two-dimensional Euler–Bernoulli node a transverse displacement and a rotation degree of freedom.
4. Form each element stiffness matrix.
5. Assemble the element matrices into the global stiffness matrix.
6. Apply support restraints and loads.
7. Solve for nodal displacements.
8. Recover element shears, moments, and deflections.
9. Refine the mesh and compare the result with an analytical solution.

For a prismatic two-dimensional bending element of length L_e, the local stiffness matrix is:

Also popular on Amazon

Structural Analysis

$19.31

Buy on Amazon

k_e = (EI/L_e³) [12, 6L_e, −12, 6L_e; 6L_e, 4L_e², −6L_e, 2L_e²; −12, −6L_e, 12, −6L_e; 6L_e, 2L_e², −6L_e, 4L_e²]

In an actual implementation, distributed loads must be converted consistently into equivalent nodal loads. An unconstrained model produces a singular stiffness matrix because it can move as a mechanism. That usually indicates missing restraints or an unstable structure, not a Python defect.

Mesh convergence

For a finite-element implementation, solve the same beam with 1, 2, 4, 8, and 16 elements. Track midspan deflection and compare it with 11.9 mm for this example. Convergence of global displacement does not automatically guarantee accurate local stresses near point loads, supports, releases, or other discontinuities. Point loads also create jumps in shear, while applied point moments create jumps in bending moment; plotting should not smooth these features away.

Choosing a Python tool

PyNite documents elastic three-dimensional structural analysis, load combinations, member force and deflection diagrams, modal analysis, P–Δ analysis, pushover analysis, and tension-only or compression-only elements. Its analysis documentation also describes linear and more advanced analysis options, including sparse and dense solvers.

OpenSeesPy documents elastic and nonlinear beam-column elements. Ansys documents Python interfaces for Mechanical and MAPDL; exact API support, product edition, version, and licensing requirements vary. These tools do not make a model automatically correct or code-compliant.

Common failure modes

Inconsistent units

Convert inputs before calculation and display a summary of the converted values. A mistake such as using 10 instead of 10,000 for a load in N/m, or using a value in mm⁴ where m⁴ is required, can change the answer by orders of magnitude.

Incorrect second moment of area

For I = bh³/12, the bending depth must be cubed. Reversing width and depth can substantially change predicted stiffness and stress.

Also popular on Amazon

Sale

Structural Analysis Made Easy: A practice book for calculating statically determined systems

$27.36

Buy on Amazon

Wrong support idealization

A pin, roller, fixed support, internal hinge, spring, and rigid link impose different constraints. Too few restraints create a mechanism; too many can produce incorrect reactions.

Sign errors

Define positive force, shear, sagging or hogging moment, and deflection before plotting. A result can be mathematically correct but appear reversed relative to a textbook convention.

Point-load placement

In an FEM model, a point load applied at a node is not the same as a distributed load over adjacent elements. Interpret local shear and moment results around the load carefully.

Linear analysis used beyond its limits

A first-order elastic script does not include geometric stiffness, P–Δ or P–δ effects, yielding, cracking, plastic hinges, composite action, or local buckling unless those behaviors are explicitly modeled. PyNite’s analysis documentation, for example, treats P–Δ analysis as a distinct capability rather than something implied by ordinary linear calculations.

When a simple Python script is not enough

Escalate the analysis when the structure involves:

• Short or deep members where shear deformation matters.
• Large deflection, second-order effects, or instability.
• Cracking, yielding, plasticity, composite action, or other nonlinear material behavior.
• Dynamic, seismic, fatigue, impact, or vibration loading.
• Three-dimensional behavior, torsion, warping, or complex connections.
• Irregular geometry, contact, local buckling, or detailed stress concentrations.
• Code-based design decisions requiring formal load combinations, resistance factors, detailing, and professional review.

Python remains valuable in these cases for preprocessing, parameter studies, automation, post-processing, and independent checks—but the underlying solver and engineering model must match the problem.

A repeatable workflow

1. Choose N–m–Pa or another fully consistent unit system.
2. Define geometry, material, section orientation, A, and I.
3. Define supports, degrees of freedom, load locations, and sign conventions.
4. Calculate reactions.
5. Calculate shear, moment, slope, and deflection.
6. Plot results with units and zero lines.
7. Check equilibrium and boundary conditions.
8. Compare at least one result with a closed-form or independent calculation.
9. For FEM, perform a mesh-convergence study.
10. Escalate to a suitable finite-element formulation and professional design review when the assumptions no longer hold.

For a basic beam, the most reliable progression is therefore theory → transparent NumPy calculation → symbolic verification with SymPy → stiffness or finite-element modeling only when the problem requires it.

Quick Recap

Bestseller No. 1

Structural Analysis

$171.56

SaleBestseller No. 2

Structural Analysis: Understanding Behavior

$76.21

SaleBestseller No. 3

Fundamentals of Structural Analysis

$56.12

Bestseller No. 4

Structural Analysis

$19.31

SaleBestseller No. 5

Structural Analysis Made Easy: A practice book for calculating statically determined systems

$27.36

Product prices and availability are accurate as of the date/time indicated and are subject to change. Any price and availability information displayed on Amazon at the time of purchase will apply.