import sympy as sp
import numpy as np
from IPython.display import Math, display
from symbolic_fem_workbench import build_bar_1d_local_problem4 1D Element Stiffness: From Symbolic to Numerical
This notebook demonstrates how to convert symbolic element matrices into efficient numerical kernels using SymPy’s lambdify. We’ll examine the numerical behavior of the 1D bar element and perform parameter studies.
4.1 Step 1: Symbolic derivation
We start with the symbolic element stiffness matrix and load vector derived in the previous notebooks. These are expressed in terms of material parameters (\(E\), \(A\)) and geometric parameters (\(L\)).
problem = build_bar_1d_local_problem()
Ke = problem["Ke"]
fe = problem["fe"]
L_sym = problem["L"]
E_sym, A_sym = problem["E"], problem["A"]
q_sym, P_sym = problem["q"], problem["P"]
print("Symbolic Ke:")
display(Math(r"K_e = " + sp.latex(Ke)))
print("\nSymbolic fe:")
display(Math(r"f_e = " + sp.latex(fe)))Symbolic Ke:\(\displaystyle K_e = \left[\begin{matrix}\frac{A E}{L} & - \frac{A E}{L}\\- \frac{A E}{L} & \frac{A E}{L}\end{matrix}\right]\)
Symbolic fe:\(\displaystyle f_e = \left[\begin{matrix}\frac{L q}{2}\\\frac{L q}{2}\end{matrix}\right]\)
4.2 Step 2: Convert to numerical functions using lambdify
SymPy’s lambdify converts symbolic expressions into NumPy-compatible functions. This allows us to rapidly evaluate matrices for different parameter values without recompiling.
ke_fn = sp.lambdify((L_sym, E_sym, A_sym), Ke, "numpy")
fe_fn = sp.lambdify((L_sym, q_sym, P_sym), fe, "numpy")
# Example: steel bar element, L=0.5m, E=200e9 Pa, A=1e-4 m^2
L_val, E_val, A_val = 0.5, 200e9, 1e-4
Ke_num = np.array(ke_fn(L_val, E_val, A_val), dtype=float)
print(f"Numerical Ke for L={L_val}, E={E_val:.0e}, A={A_val}:")
print(Ke_num)
q_val, P_val = 1000.0, 500.0
fe_num = np.array(fe_fn(L_val, q_val, P_val), dtype=float).flatten()
print(f"\nNumerical fe for q={q_val}, P={P_val}:")
print(fe_num)Numerical Ke for L=0.5, E=2e+11, A=0.0001:
[[ 40000000. -40000000.]
[-40000000. 40000000.]]
Numerical fe for q=1000.0, P=500.0:
[250. 250.]4.3 Step 3: Verify properties
The element stiffness matrix should satisfy important mathematical properties: - Symmetry: \(K_e = K_e^T\) (from the symmetry of the bilinear form) - Positive semi-definiteness: All eigenvalues \(\geq 0\) - Singular: One null eigenvector (rigid body translation)
print("Symmetry check: Ke == Ke.T?", np.allclose(Ke_num, Ke_num.T))
print("Row sums (should be zero for homogeneous problem):", Ke_num.sum(axis=1))
print("Positive diagonal entries?", np.all(np.diag(Ke_num) > 0))
# Eigenvalue analysis
eigvals = np.linalg.eigvalsh(Ke_num)
print(f"\nEigenvalues: {eigvals}")
print(f"Smallest eigenvalue (should be ~0): {eigvals[0]:.2e}")
print(f"Positive semi-definite: {np.all(eigvals >= -1e-10)}")Symmetry check: Ke == Ke.T? True
Row sums (should be zero for homogeneous problem): [0. 0.]
Positive diagonal entries? True
Eigenvalues: [ 0. 80000000.]
Smallest eigenvalue (should be ~0): 0.00e+00
Positive semi-definite: True4.4 Step 4: Parameter study
The stiffness scales with the ratio \(EA/L\). Let’s verify this behavior numerically.
print("Effect of element length on stiffness:")
for L_test in [0.1, 0.5, 1.0, 2.0]:
K = np.array(ke_fn(L_test, E_val, A_val), dtype=float)
print(f" L = {L_test:4.1f} -> K[0,0] = {K[0,0]:.2e} (EA/L = {E_val*A_val/L_test:.2e})")
print("\nEffect of cross-section on stiffness:")
for A_test in [1e-5, 1e-4, 5e-4, 1e-3]:
K = np.array(ke_fn(L_val, E_val, A_test), dtype=float)
print(f" A = {A_test:.0e} -> K[0,0] = {K[0,0]:.2e}")Effect of element length on stiffness:
L = 0.1 -> K[0,0] = 2.00e+08 (EA/L = 2.00e+08)
L = 0.5 -> K[0,0] = 4.00e+07 (EA/L = 4.00e+07)
L = 1.0 -> K[0,0] = 2.00e+07 (EA/L = 2.00e+07)
L = 2.0 -> K[0,0] = 1.00e+07 (EA/L = 1.00e+07)
Effect of cross-section on stiffness:
A = 1e-05 -> K[0,0] = 4.00e+06
A = 1e-04 -> K[0,0] = 4.00e+07
A = 5e-04 -> K[0,0] = 2.00e+08
A = 1e-03 -> K[0,0] = 4.00e+084.5 Exercise
Question: Assemble a 3-element bar manually using the numerical kernels: 1. Create 3 elements of equal length \(L_{\text{total}}/3\) 2. Assemble the global \(4 \times 4\) stiffness matrix 3. Apply boundary condition \(u(0) = 0\) and a prescribed force \(P\) at \(x=L_{\text{total}}\) 4. Solve for displacements and plot the solution
Hint: Use np.zeros to initialize the global matrix and an assembly loop with connectivity information.
# Exercise: Assemble a 3-element bar manually.
# Create 3 elements of equal length L_total/3, assemble K_global (4x4),
# apply u(0) = 0 and P at x=L_total, solve, and plot the solution.