Python for Engineers

Python for Engineers: Code Repository

Explore practical Python scripts for engineering calculations, numerical analysis, differential equations, signal processing, and engineering data processing using NumPy, SciPy, Pandas, and Matplotlib.

Spring-Mass-Damper ODE Solver

SciPyODEDynamics

Solves the second-order ordinary differential equation of a 1D spring-mass-damper system using SciPy.

Spring-Mass-Damper ODE Solver

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# System parameters
m = 10.0      # Mass (kg)
k = 250.0     # Spring constant (N/m)
c = 15.0      # Damping coefficient (N·s/m)

# Free vibration: F(t) = 0
def system_dynamics(t, y):
    x, v = y

    dxdt = v
    dvdt = -(c / m) * v - (k / m) * x

    return [dxdt, dvdt]

# Initial conditions
y0 = [0.5, 0.0]

t_span = (0, 10)
t_eval = np.linspace(t_span[0], t_span[1], 500)

solution = solve_ivp(
    system_dynamics,
    t_span,
    y0,
    t_eval=t_eval
)

plt.plot(solution.t, solution.y[0], label='Displacement x(t)')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Damped Free Vibration')
plt.grid(True)
plt.legend()
plt.show()

Forced Vibration Response

NumPyDynamicsVibration

Calculates and plots the steady-state response of a damped single-degree-of-freedom system under harmonic excitation.

Forced Vibration Response

import numpy as np
import matplotlib.pyplot as plt

# System parameters
m = 10.0       # kg
k = 1000.0     # N/m
zeta = 0.05    # Damping ratio
F0 = 100.0     # Force amplitude (N)

wn = np.sqrt(k / m)
c = 2 * zeta * m * wn

frequency = np.linspace(0.1, 30, 1000)
omega = 2 * np.pi * frequency

# Frequency response
H = 1 / np.sqrt(
    (k - m * omega**2)**2 +
    (c * omega)**2
)

amplitude = F0 * H

plt.plot(frequency, amplitude)
plt.xlabel('Excitation Frequency (Hz)')
plt.ylabel('Displacement Amplitude (m)')
plt.title('Forced Vibration Frequency Response')
plt.grid(True)
plt.show()

Natural Frequency Calculator

NumPyDynamicsFrequency

Calculates natural frequency and period for a simple mass-spring system.

Natural Frequency Calculator

import numpy as np

mass = 10.0       # kg
stiffness = 2500  # N/m

omega_n = np.sqrt(stiffness / mass)
frequency = omega_n / (2 * np.pi)
period = 1 / frequency

print(f"Natural angular frequency: {omega_n:.3f} rad/s")
print(f"Natural frequency: {frequency:.3f} Hz")
print(f"Period: {period:.4f} s")

Vibration Signal FFT Analysis

NumPyFFTVibration

Applies Fast Fourier Transform to a time-domain vibration signal to identify dominant frequency components.

Vibration Signal FFT Analysis

import numpy as np
import matplotlib.pyplot as plt

# Sampling parameters
fs = 1000
t = np.arange(0, 1.0, 1 / fs)

# Signal components
freq1 = 50
freq2 = 120

noise = np.random.normal(0, 0.5, t.shape)

signal = (
    np.sin(2 * np.pi * freq1 * t)
    + 0.5 * np.sin(2 * np.pi * freq2 * t)
    + noise
)

# FFT
n = len(signal)

fft_values = np.fft.fft(signal)
fft_freqs = np.fft.fftfreq(n, 1 / fs)

positive = fft_freqs > 0

freqs = fft_freqs[positive]
magnitudes = np.abs(fft_values[positive]) / (n / 2)

plt.plot(freqs, magnitudes)
plt.xlim(0, 200)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.title('Vibration Signal FFT')
plt.grid(True)
plt.show()

Vibration Low-Pass Filter

SciPyFilterSignal Processing

Applies a Butterworth low-pass filter to remove high-frequency noise from a measured engineering signal.

Vibration Low-Pass Filter

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt

fs = 1000
t = np.arange(0, 2, 1 / fs)

# Example signal
signal = (
    np.sin(2 * np.pi * 20 * t)
    + 0.4 * np.sin(2 * np.pi * 250 * t)
    + 0.2 * np.random.randn(len(t))
)

# Butterworth filter
cutoff = 50
order = 4

b, a = butter(
    order,
    cutoff / (fs / 2),
    btype='low'
)

filtered = filtfilt(b, a, signal)

plt.plot(t, signal, alpha=0.4, label='Raw')
plt.plot(t, filtered, label='Filtered')

plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('Low-Pass Filter')
plt.legend()
plt.grid(True)
plt.show()

Simply Supported Beam Bending

NumPyBeamStressDeflection

Calculates bending moment, shear force and deflection for a simply supported beam with a central point load.

Simply Supported Beam Bending

import numpy as np
import matplotlib.pyplot as plt

# Beam parameters
L = 2.0          # Beam length (m)
P = 1000.0       # Point load (N)
E = 210e9        # Young's modulus (Pa)
I = 8.0e-6       # Second moment of area (m^4)

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

# Reactions
RA = P / 2
RB = P / 2

# Shear force
V = np.where(x < L / 2, RA, -RB)

# Bending moment
M = np.where(
    x <= L / 2,
    RA * x,
    RB * (L - x)
)

# Deflection
delta = np.where(
    x <= L / 2,
    P * x * (3 * L**2 - 4 * x**2) / (48 * E * I),
    P * (L - x) * (3 * L**2 - 4 * (L - x)**2)
    / (48 * E * I)
)

plt.plot(x, M)
plt.xlabel('Position (m)')
plt.ylabel('Bending Moment (N·m)')
plt.title('Simply Supported Beam')
plt.grid(True)
plt.show()

print(f"Maximum bending moment: {P * L / 4:.2f} N·m")
print(f"Maximum deflection: {np.max(delta) * 1000:.3f} mm")

Circular Shaft Torsion Analysis

TorsionShaftStress

Calculates shear stress and angle of twist for a circular shaft subjected to torque.

Circular Shaft Torsion Analysis

import numpy as np

# Shaft parameters
T = 1500.0       # Torque (N·m)
L = 1.2          # Shaft length (m)
d = 0.04         # Diameter (m)
G = 80e9         # Shear modulus (Pa)

# Polar moment of inertia
J = np.pi * d**4 / 32

# Maximum shear stress
tau_max = T * (d / 2) / J

# Angle of twist
theta = T * L / (G * J)

print(f"Polar moment J: {J:.4e} m^4")
print(f"Maximum shear stress: {tau_max / 1e6:.2f} MPa")
print(f"Angle of twist: {np.degrees(theta):.3f} degrees")

Euler Column Buckling

BucklingEulerStructures

Calculates the critical Euler buckling load for a slender compression member.

Euler Column Buckling

import numpy as np

# Column parameters
E = 210e9       # Pa
I = 8.0e-6      # m^4
L = 2.0         # m

# Effective length factor
K = 1.0         # Pinned-pinned

# Euler critical load
Pcr = (
    np.pi**2 * E * I
    / (K * L)**2
)

print(f"Critical buckling load: {Pcr / 1000:.2f} kN")

Von Mises Stress Calculator

StressVon MisesFEA

Calculates equivalent Von Mises stress from a 3D stress state.

Von Mises Stress Calculator

import numpy as np

# Normal stresses
sigma_x = 120.0  # MPa
sigma_y = 60.0   # MPa
sigma_z = 0.0    # MPa

# Shear stresses
tau_xy = 20.0    # MPa
tau_yz = 10.0    # MPa
tau_zx = 15.0    # MPa

von_mises = np.sqrt(
    0.5 * (
        (sigma_x - sigma_y)**2
        + (sigma_y - sigma_z)**2
        + (sigma_z - sigma_x)**2
    )
    + 3 * (
        tau_xy**2
        + tau_yz**2
        + tau_zx**2
    )
)

print(f"Von Mises stress: {von_mises:.2f} MPa")

Linear Thermal Expansion

ThermalMaterialsExpansion

Calculates dimensional change caused by temperature variation using the linear thermal expansion equation.

Linear Thermal Expansion

# Linear thermal expansion

L0 = 1.0             # Initial length (m)
alpha = 12e-6        # Thermal expansion coefficient (1/K)
delta_T = 80.0       # Temperature change (°C)

delta_L = alpha * L0 * delta_T
L_final = L0 + delta_L

print(f"Length change: {delta_L * 1000:.3f} mm")
print(f"Final length: {L_final:.6f} m")

1D Heat Conduction

Heat TransferFourierThermal

Calculates heat transfer through a flat wall using Fourier conduction.

1D Heat Conduction

# Fourier's law of conduction

k = 45.0       # Thermal conductivity (W/m·K)
A = 0.02       # Area (m²)
L = 0.01       # Wall thickness (m)

T_hot = 120.0
T_cold = 30.0

Q = k * A * (T_hot - T_cold) / L

print(f"Heat transfer rate: {Q:.2f} W")

Reynolds Number Analysis

Fluid MechanicsReynoldsPipe Flow

Calculates Reynolds number and classifies internal pipe flow as laminar or turbulent.

Reynolds Number Analysis

# Reynolds number

rho = 998.0       # Density (kg/m³)
V = 2.0           # Velocity (m/s)
D = 0.05          # Pipe diameter (m)
mu = 0.001        # Dynamic viscosity (Pa·s)

Re = rho * V * D / mu

print(f"Reynolds number: {Re:.0f}")

if Re < 2300:
    print("Flow regime: Laminar")
elif Re < 4000:
    print("Flow regime: Transitional")
else:
    print("Flow regime: Turbulent")

Darcy-Weisbach Pressure Drop

Pipe FlowDarcy-WeisbachPressure Loss

Calculates pressure loss in a straight pipe using the Darcy-Weisbach equation.

Darcy-Weisbach Pressure Drop

import numpy as np

rho = 1000.0      # kg/m³
V = 2.0           # m/s
D = 0.05          # m
L = 20.0          # m
f = 0.025         # Darcy friction factor

delta_p = f * (L / D) * (rho * V**2 / 2)

print(f"Pressure drop: {delta_p / 1000:.2f} kPa")

Ohm's Law & Power Calculator

ElectricalOhm LawPower

Calculates voltage, current, resistance and electrical power using Ohm’s law.

Ohm's Law & Power Calculator

# Ohm's Law

V = 24.0       # Voltage (V)
R = 12.0       # Resistance (Ohm)

I = V / R
P = V * I

print(f"Current: {I:.2f} A")
print(f"Power: {P:.2f} W")

Voltage Divider Calculator

ElectricalVoltage DividerCircuits

Calculates output voltage from a two-resistor voltage divider.

Voltage Divider Calculator

# Voltage divider

Vin = 12.0
R1 = 10000.0
R2 = 5000.0

Vout = Vin * R2 / (R1 + R2)

print(f"Output voltage: {Vout:.3f} V")

RC Circuit Transient Response

RC CircuitTransientNumPy

Simulates the charging response of a first-order RC circuit.

RC Circuit Transient Response

import numpy as np
import matplotlib.pyplot as plt

V = 5.0
R = 10000.0
C = 100e-6

tau = R * C

t = np.linspace(0, 5 * tau, 500)

voltage = V * (
    1 - np.exp(-t / tau)
)

plt.plot(t, voltage)
plt.xlabel('Time (s)')
plt.ylabel('Capacitor Voltage (V)')
plt.title('RC Charging Curve')
plt.grid(True)
plt.show()

print(f"Time constant: {tau:.3f} s")

Nonlinear Equation Solver

SciPyNumericalRoot Finding

Finds the root of a nonlinear engineering equation using SciPy’s root-finding algorithms.

Nonlinear Equation Solver

import numpy as np
from scipy.optimize import brentq

# Example nonlinear equation
def equation(x):
    return x**3 - 5*x - 2

root = brentq(
    equation,
    0,
    3
)

print(f"Root: {root:.6f}")
print(f"Equation value: {equation(root):.3e}")

Engineering Curve Fitting

SciPyCurve FitExperimental Data

Fits experimental engineering data to a nonlinear model using SciPy.

Engineering Curve Fitting

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# Experimental data
x = np.array([1, 2, 3, 4, 5, 6])
y = np.array([2.1, 4.3, 6.2, 8.4, 10.1, 12.5])

# Linear model
def model(x, a, b):
    return a * x + b

params, covariance = curve_fit(
    model,
    x,
    y
)

a, b = params

x_fit = np.linspace(1, 6, 200)
y_fit = model(x_fit, a, b)

plt.scatter(x, y, label='Experimental')
plt.plot(x_fit, y_fit, label='Fitted model')

plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

print(f"Slope: {a:.4f}")
print(f"Intercept: {b:.4f}")

Engineering Measurement Statistics

NumPyStatisticsMeasurements

Calculates mean, standard deviation, variance and basic statistics from experimental measurements.

Engineering Measurement Statistics

import numpy as np

measurements = np.array([
    10.12,
    10.08,
    10.15,
    10.10,
    10.11,
    10.09,
    10.14
])

mean = np.mean(measurements)
std = np.std(measurements, ddof=1)
variance = np.var(measurements, ddof=1)

print(f"Mean: {mean:.4f}")
print(f"Standard deviation: {std:.4f}")
print(f"Variance: {variance:.6f}")

Gear Ratio & Torque Calculator

GearsTorqueMechanical

Calculates gear ratio, output speed and output torque for a simple gear train.

Gear Ratio & Torque Calculator

# Gear ratio calculation

input_speed = 3000.0     # rpm
input_torque = 20.0      # N·m

driver_teeth = 20
driven_teeth = 60

ratio = driven_teeth / driver_teeth

output_speed = input_speed / ratio
output_torque = input_torque * ratio

print(f"Gear ratio: {ratio:.2f}:1")
print(f"Output speed: {output_speed:.1f} rpm")
print(f"Output torque: {output_torque:.1f} N·m")

Vehicle Braking Distance

Vehicle DynamicsBrakingAutomotive

Estimates theoretical braking distance from vehicle speed and tire-road friction coefficient.

Vehicle Braking Distance

# Theoretical braking distance

import numpy as np

speed_kmh = 100.0
mu = 0.8
g = 9.81

speed_ms = speed_kmh / 3.6

braking_distance = speed_ms**2 / (2 * mu * g)

print(f"Initial speed: {speed_kmh:.1f} km/h")
print(f"Friction coefficient: {mu:.2f}")
print(f"Braking distance: {braking_distance:.2f} m")

Wheel Torque & Power

AutomotivePowerTorque

Calculates wheel force, wheel torque and power from vehicle speed and motor output.

Wheel Torque & Power

import numpy as np

motor_power = 50000.0     # W
wheel_radius = 0.32       # m
vehicle_speed = 80.0      # km/h

speed_ms = vehicle_speed / 3.6

wheel_force = motor_power / speed_ms
wheel_torque = wheel_force * wheel_radius

print(f"Vehicle speed: {vehicle_speed:.1f} km/h")
print(f"Wheel force: {wheel_force:.1f} N")
print(f"Wheel torque: {wheel_torque:.1f} N·m")

Python for Engineering Analysis

Python has become a widely used programming language for engineering analysis, automation, numerical computing, and technical data processing. Its clear syntax and extensive scientific computing ecosystem make it useful for automating calculations, analyzing experimental data, and developing numerical models.

NumPy: The Foundation of Numerical Computing

For engineers working with structural matrices, kinematics, numerical calculations, or large sensor datasets, NumPy provides an efficient foundation for numerical computing. Its multidimensional arrays and mathematical functions are useful for matrix operations, vector calculations, data transformation, and engineering algorithms.

SciPy: Advanced Engineering Mathematics

SciPy extends Python with algorithms for optimization, integration, interpolation, eigenvalue problems, and ordinary differential equations (ODEs). Engineers can use functions such as solve_ivp to model dynamic systems including mass-spring-damper systems and other time-dependent engineering problems.

Pandas and Engineering Data Processing

Engineering workflows often involve data from simulations, experiments, measurements, and sensors. Pandas provides tools for importing, cleaning, filtering, transforming, and analyzing structured datasets such as CSV and Excel files. This makes it useful for preparing engineering data before statistical analysis, visualization, or numerical modeling.

Signal Processing and FFT Analysis

Python can also be used for engineering signal processing, including noise filtering, frequency-domain analysis, and Fast Fourier Transform (FFT) calculations. These techniques are commonly applied to vibration measurements, sensor signals, rotating machinery data, and other time-series engineering datasets.

Continue Calculating

reference

Engineering Formulas Cheat Sheet

A quick-reference library of fundamental physics and mechanical engineering equations.

Open

reference

Engineering Materials Database

Interactive database to compare yield strength, tensile strength, and elastic modulus.

Open

calculator

Von Mises & Tresca Yield Calculator

Evaluate a 2D plane stress state against material failure theories.

Open

Frequently Asked Questions

Do I need to install any external libraries to run these scripts?

Most of these engineering scripts rely on the scientific Python stack. Depending on the script, you may need NumPy, SciPy, Pandas, and Matplotlib. These packages can typically be installed using a Python package manager such as pip.

Can I use these snippets in Jupyter Notebooks?

Yes. The scripts can be adapted for Jupyter Notebook environments, making them useful for interactive engineering analysis, documentation, experimentation, and data visualization.