import sympy as sp
from IPython.display import Math, display
from symbolic_fem_workbench import (
make_field_1d, weighted_residual, integrate_divergence_1d,
drop_dirichlet_boundary, apply_neumann_flux, split_linear_weak_form,
)2 From Strong Form to Weak Form
This notebook demonstrates the fundamental transformation from a strong form PDE into a weak form, which is the mathematical foundation of the Finite Element Method. We will convert a 1D bar equilibrium equation into its weak form using the method of weighted residuals and integration by parts.
2.1 Step 1: Define the problem
We consider a 1D bar under axial loading. The strong form of the equilibrium equation is: \[\frac{d}{dx}\left(EA\frac{du}{dx}\right) + q = 0\]
where: - \(u(x)\) is the axial displacement - \(E\) is Young’s modulus - \(A\) is the cross-sectional area - \(q\) is a distributed load - Domain: \(x \in [0, L]\)
x, L = sp.symbols("x L", positive=True)
E, A = sp.symbols("E A", positive=True)
q, P = sp.symbols("q P", real=True)
u = make_field_1d("u", x)
v = make_field_1d("v", x)
domain = (x, 0, L)
print("Trial field u(x):", u)
print("Test field v(x):", v)
print("Domain:", domain)Trial field u(x): u(x)
Test field v(x): v(x)
Domain: (x, 0, L)2.2 Step 2: Build the weighted residual
We multiply the strong form by a test function \(v(x)\) (also called a weight function) and integrate over the domain. This is the method of weighted residuals: \[\int_0^L \left[\frac{d}{dx}\left(EA\frac{du}{dx}\right) + q\right] v \, dx = 0\]
flux = E * A * sp.diff(u, x)
strong_form_lhs = sp.diff(flux, x) + q
display(Math(r"\text{Strong form: } " + sp.latex(strong_form_lhs) + " = 0"))
wr = weighted_residual(strong_form_lhs, 0, v, domain)
display(Math(r"\text{Weighted residual: } " + sp.latex(wr.as_expression()) + " = 0"))\(\displaystyle \text{Strong form: } A E \frac{d^{2}}{d x^{2}} u{\left(x \right)} + q = 0\)
\(\displaystyle \text{Weighted residual: } \int\limits_{0}^{L} \left(A E v{\left(x \right)} \frac{d^{2}}{d x^{2}} u{\left(x \right)} + q v{\left(x \right)}\right)\, dx = 0\)
2.3 Step 3: Integration by parts
We now apply integration by parts to shift one derivative from \(u\) to \(v\). This lowers the regularity requirement on \(u\) and introduces natural boundary conditions: \[\int_0^L \frac{d}{dx}\left(EA\frac{du}{dx}\right) v \, dx = \left[EA\frac{du}{dx}v\right]_0^L - \int_0^L EA\frac{du}{dx}\frac{dv}{dx} \, dx\]
domain_term, (left_bc, right_bc) = integrate_divergence_1d(flux, v, x, domain)
print("Domain integral (after IBP):")
display(Math(sp.latex(domain_term.as_integral())))
print("\nLeft boundary term:")
display(Math(sp.latex(left_bc.evaluate())))
print("\nRight boundary term:")
display(Math(sp.latex(right_bc.evaluate())))Domain integral (after IBP):\(\displaystyle \int\limits_{0}^{L} A E \frac{d}{d x} u{\left(x \right)} \frac{d}{d x} v{\left(x \right)}\, dx\)
Left boundary term:\(\displaystyle - A E v{\left(0 \right)} \left. \frac{d}{d x} u{\left(x \right)} \right|_{\substack{ x=0 }}\)
Right boundary term:\(\displaystyle - A E v{\left(L \right)} \frac{d}{d L} u{\left(L \right)}\)
2.4 Step 4: Apply boundary conditions
Dirichlet BC (Essential): The test function must vanish where displacement is prescribed. At \(x=0\): \(u(0) = 0\) (homogeneous Dirichlet) \(\Rightarrow\) \(v(0) = 0\)
Neumann BC (Natural): The prescribed flux/traction appears as a source term. At \(x=L\): Prescribed force \(P\) \(\Rightarrow\) \(EA\frac{du}{dx}\big|_x=L = P\)
# Homogeneous Dirichlet at x=0: test function vanishes
left_val = drop_dirichlet_boundary(left_bc, v)
print("Left BC (Dirichlet, v=0 at x=0):", left_val)
# Neumann at x=L: prescribed traction P
right_val = apply_neumann_flux(right_bc, flux_expr=flux, prescribed_flux=P)
print("Right BC (Neumann, flux=P at x=L):", right_val)Left BC (Dirichlet, v=0 at x=0): 0
Right BC (Neumann, flux=P at x=L): -P*v(L)2.5 Step 5: Assemble and split the weak form
We now assemble all terms into a single weak form and decompose it into bilinear and linear parts: \[a(u, v) + F(v) = 0\] where \(a(u,v)\) is bilinear in \(u\) and \(v\), and \(F(v)\) is linear in \(v\).
source_integral = sp.Integral(q * v, domain)
weak_expr = domain_term.as_integral() - source_integral + left_val + right_val
display(Math(r"\text{Full weak form expression: } " + sp.latex(weak_expr) + " = 0"))
weak_form = split_linear_weak_form(weak_expr, u, v)
print("\nBilinear form a(u,v):")
display(Math(sp.latex(weak_form.bilinear)))
print("\nLinear form F(v):")
display(Math(sp.latex(weak_form.linear)))\(\displaystyle \text{Full weak form expression: } - P v{\left(L \right)} - \int\limits_{0}^{L} q v{\left(x \right)}\, dx + \int\limits_{0}^{L} A E \frac{d}{d x} u{\left(x \right)} \frac{d}{d x} v{\left(x \right)}\, dx = 0\)
Bilinear form a(u,v):\(\displaystyle A E \int\limits_{0}^{L} \frac{d}{d x} u{\left(x \right)} \frac{d}{d x} v{\left(x \right)}\, dx\)
Linear form F(v):\(\displaystyle q \int\limits_{0}^{L} v{\left(x \right)}\, dx\)
2.6 Exercise
Question: Modify the problem to have Dirichlet conditions on both ends (\(u(0) = 0\), \(u(L) = 0\)) with a distributed source term \(q\). How does the boundary condition handling change? What happens to the natural boundary term contributions?
# Exercise: Modify the problem to have Dirichlet conditions on both ends:
# u(0) = 0, u(L) = 0, with source term q.
# What changes in the boundary condition handling?