Work through these hands-on notebooks alongside this chapter:
```{python}
#| code-fold: true
fem.display_strong_form()
```Strong form of the equation:
\(\displaystyle A{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)} + B{\left(x \right)} \frac{d}{d x} y{\left(x \right)} + C{\left(x \right)} y{\left(x \right)} = F{\left(x \right)}\)
With boundary conditions:
\(\displaystyle y{\left(0 \right)} = 0\)
\(\displaystyle y{\left(L \right)} = 0\)
where x ∈ [0, L]
```{python}
#| code-fold: true
weak_eq = fem.derive_weak_form()
```Step 1: Multiply the strong form by a test function v(x):
EXPLANATION: We start with the strong form of our differential equation and multiply both sides by a test function v(x)
This is the first step in converting to the weak form.
The test function v(x) will eventually vanish at the boundaries, helping us incorporate boundary conditions naturally.
\(\displaystyle \left(A{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)} + B{\left(x \right)} \frac{d}{d x} y{\left(x \right)} + C{\left(x \right)} y{\left(x \right)}\right) v{\left(x \right)} = F{\left(x \right)} v{\left(x \right)}\)
Step 2: Integrate over the domain [0, L]:
EXPLANATION: We integrate both sides of the equation over the entire domain [0, L].
This transforms our pointwise equation into an integral equation that will hold over the whole domain.
This step moves us from a local formulation to a global one, which is a key aspect of the finite element method.
\(\displaystyle \int\limits_{0}^{L} \left(A{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)} + B{\left(x \right)} \frac{d}{d x} y{\left(x \right)} + C{\left(x \right)} y{\left(x \right)}\right) v{\left(x \right)}\, dx = \int\limits_{0}^{L} F{\left(x \right)} v{\left(x \right)}\, dx\)
Step 3: Expand the left-hand side:
EXPLANATION:
We expand the left-hand side to separate the terms with different derivatives of y. This makes it easier to apply integration by parts to the second-derivative term.
Breaking down the differential operator allows us to handle each term appropriately based on the order of derivatives.
\(\displaystyle \int\limits_{0}^{L} A{\left(x \right)} v{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)}\, dx + \int\limits_{0}^{L} B{\left(x \right)} v{\left(x \right)} \frac{d}{d x} y{\left(x \right)}\, dx + \int\limits_{0}^{L} C{\left(x \right)} v{\left(x \right)} y{\left(x \right)}\, dx = \int\limits_{0}^{L} F{\left(x \right)} v{\left(x \right)}\, dx\)
Step 4: Integration by parts for the term with second derivative:
EXPLANATION:
This is the critical step that reduces the order of derivatives in our equation. We apply integration by parts to the term containing the second derivative of y.
The formula for integration by parts is:
∫u·dv = [u·v] - ∫v·du
For our case:
u = v(x)·A(x) and dv = d²y/dx² dx
du = d(v(x)·A(x))/dx dx and v = dy/dx
Using the product rule:
d/dx(v·A·dy/dx) = v·A·d²y/dx² + (v·A)’·dy/dx
Rearranging:
v·A·d²y/dx² = d/dx(v·A·dy/dx) - (v·A)’·dy/dx
This substitution allows us to replace the second derivative with first derivatives, reducing the continuity requirements on our approximation functions.
After integration by parts:
\(\displaystyle \int\limits_{0}^{L} A{\left(x \right)} v{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)}\, dx = - A{\left(0 \right)} v{\left(0 \right)} \left. \frac{d}{d x} y{\left(x \right)} \right|_{\substack{ x=0 }} + A{\left(L \right)} v{\left(L \right)} \frac{d}{d L} y{\left(L \right)} - \int\limits_{0}^{L} \left(A{\left(x \right)} \frac{d}{d x} v{\left(x \right)} + v{\left(x \right)} \frac{d}{d x} A{\left(x \right)}\right) \frac{d}{d x} y{\left(x \right)}\, dx\)
Step 5: Apply boundary conditions v(0) = v(L) = 0:
EXPLANATION:
We choose our test functions v(x) to vanish at the domain boundaries (v(0) = v(L) = 0).
This is a crucial step that eliminates the boundary terms from integration by parts. By carefully selecting test functions that satisfy these conditions, we naturally incorporate the essential (Dirichlet) boundary conditions into our formulation.
This is one of the elegant aspects of the finite element method.
The boundary term [v·A·dy/dx]₀ᴸ vanishes because v(0) = v(L) = 0.
Step 6: Assemble the final weak form:
EXPLANATION: Now we collect all terms to form the complete weak formulation.
This form requires lower continuity requirements on our solution (only first derivatives appear).
The weak form is equivalent to the strong form for sufficiently smooth solutions, but it allows us to find approximate solutions with less smoothness.
This is the foundation for the finite element approximation where we’ll represent the solution as a combination of piecewise polynomial functions.
\(\displaystyle - \int\limits_{0}^{L} \left(A{\left(x \right)} \frac{d}{d x} v{\left(x \right)} + v{\left(x \right)} \frac{d}{d x} A{\left(x \right)}\right) \frac{d}{d x} y{\left(x \right)}\, dx + \int\limits_{0}^{L} B{\left(x \right)} v{\left(x \right)} \frac{d}{d x} y{\left(x \right)}\, dx + \int\limits_{0}^{L} C{\left(x \right)} v{\left(x \right)} y{\left(x \right)}\, dx = \int\limits_{0}^{L} F{\left(x \right)} v{\left(x \right)}\, dx\)
```{python}
#| code-fold: true
N = 5
phi, a, b, y_approx, v_approx = fem.define_shape_functions(N)
```Step 7: Define shape functions for N elements:
EXPLANATION:
Now we discretize the problem by representing both the solution y(x) and test function v(x) as linear combinations of shape functions.
This transforms our continuous problem into a discrete one with a finite number of unknowns (the coefficients).
The shape functions are typically chosen to be simple piecewise polynomials (often linear) that are non-zero only over a small portion of the domain .
This ‘local support’ property often leads to sparse matrices in the final system.
\(\displaystyle \left[ \phi_{1}{\left(x \right)}, \ \phi_{2}{\left(x \right)}, \ \phi_{3}{\left(x \right)}, \ \phi_{4}{\left(x \right)}, \ \phi_{5}{\left(x \right)}\right]\)
Step 8: Express solution and test functions using shape functions:
EXPLANATION:
We approximate both the solution y(x) and test function v(x) using the same set of shape functions.
The solution is written as a sum of shape functions with unknown coefficients a_j. Similarly, the test function is written as a sum with coefficients b_i.
This is the Galerkin approach, where we use the same functions for both approximation and testing.
\(\displaystyle y{\left(x \right)} = a_{1} \phi_{1}{\left(x \right)} + a_{2} \phi_{2}{\left(x \right)} + a_{3} \phi_{3}{\left(x \right)} + a_{4} \phi_{4}{\left(x \right)} + a_{5} \phi_{5}{\left(x \right)}\)
\(\displaystyle v{\left(x \right)} = b_{1} \phi_{1}{\left(x \right)} + b_{2} \phi_{2}{\left(x \right)} + b_{3} \phi_{3}{\left(x \right)} + b_{4} \phi_{4}{\left(x \right)} + b_{5} \phi_{5}{\left(x \right)}\)
Where:
NOTE: In the finite element method, we often choose shape functions that are equal to 1 at their corresponding node and 0 at all other nodes.
This makes it easy to enforce boundary conditions and interpret the coefficients aⱼ as the solution values at the nodes.
```{python}
#| code-fold: true
K, F_vec = fem.assemble_fem_system(weak_eq, phi, a, b, y_approx, v_approx, N)
```** Step 9: Substitute discretized functions into the weak form:**
EXPLANATION: We now substitute our approximations for y(x) and v(x) into the weak form.
This transforms the continuous weak form into a discrete algebraic system of equations.
Since the test function coefficients bᵢ are arbitrary, we can collect terms and set up a system of equations for the unknown coefficients aⱼ.
Step 10: Assemble the system matrix and right-hand side vector:
EXPLANATION:
When we substitute the discretized functions and collect terms, we obtain a linear system of equations K·a = F, where:
Each entry K[i,j] is computed by evaluating an integral that involves the shape functions phi_i and phi_j and their derivatives.
Similarly, each entry F[i] comes from an integral involving phi_i and the force function F(x).
For a simplified case where A(x) = B(x) = C(x) = 1, the entries would be:
\(\displaystyle K_{ij} = \int_0^L \frac{d\phi_i}{dx} \frac{d\phi_j}{dx} dx + \int_0^L \phi_i(x) \frac{d\phi_j}{dx} dx + \int_0^L \phi_i(x) \phi_j(x) dx\)
\(\displaystyle F_i = \int_0^L \phi_i(x) F(x) dx\)
NOTE: In practice, the assembly process is typically done by looping over elements,computing the element contributions, and adding them to the global system.
This element-by-element approach is more efficient and aligns with the local support property of shape functions.
```{python}
#| code-fold: true
fem.finite_element_method_1d_demo()
```======================================================================
INTRODUCTION TO THE FINITE ELEMENT METHOD
Lecture #1: Introduction + 1D Setup
======================================================================
Strong form of the equation:
\(\displaystyle A{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)} + B{\left(x \right)} \frac{d}{d x} y{\left(x \right)} + C{\left(x \right)} y{\left(x \right)} = F{\left(x \right)}\)
With boundary conditions:
\(\displaystyle y{\left(0 \right)} = 0\)
\(\displaystyle y{\left(L \right)} = 0\)
where x ∈ [0, L]
======================================================================
DERIVING THE WEAK FORM
======================================================================
Step 1: Multiply the strong form by a test function v(x):
EXPLANATION: We start with the strong form of our differential equation and multiply both sides by a test function v(x)
This is the first step in converting to the weak form.
The test function v(x) will eventually vanish at the boundaries, helping us incorporate boundary conditions naturally.
\(\displaystyle \left(A{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)} + B{\left(x \right)} \frac{d}{d x} y{\left(x \right)} + C{\left(x \right)} y{\left(x \right)}\right) v{\left(x \right)} = F{\left(x \right)} v{\left(x \right)}\)
Step 2: Integrate over the domain [0, L]:
EXPLANATION: We integrate both sides of the equation over the entire domain [0, L].
This transforms our pointwise equation into an integral equation that will hold over the whole domain.
This step moves us from a local formulation to a global one, which is a key aspect of the finite element method.
\(\displaystyle \int\limits_{0}^{L} \left(A{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)} + B{\left(x \right)} \frac{d}{d x} y{\left(x \right)} + C{\left(x \right)} y{\left(x \right)}\right) v{\left(x \right)}\, dx = \int\limits_{0}^{L} F{\left(x \right)} v{\left(x \right)}\, dx\)
Step 3: Expand the left-hand side:
EXPLANATION:
We expand the left-hand side to separate the terms with different derivatives of y. This makes it easier to apply integration by parts to the second-derivative term.
Breaking down the differential operator allows us to handle each term appropriately based on the order of derivatives.
\(\displaystyle \int\limits_{0}^{L} A{\left(x \right)} v{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)}\, dx + \int\limits_{0}^{L} B{\left(x \right)} v{\left(x \right)} \frac{d}{d x} y{\left(x \right)}\, dx + \int\limits_{0}^{L} C{\left(x \right)} v{\left(x \right)} y{\left(x \right)}\, dx = \int\limits_{0}^{L} F{\left(x \right)} v{\left(x \right)}\, dx\)
Step 4: Integration by parts for the term with second derivative:
EXPLANATION:
This is the critical step that reduces the order of derivatives in our equation. We apply integration by parts to the term containing the second derivative of y.
The formula for integration by parts is:
∫u·dv = [u·v] - ∫v·du
For our case:
u = v(x)·A(x) and dv = d²y/dx² dx
du = d(v(x)·A(x))/dx dx and v = dy/dx
Using the product rule:
d/dx(v·A·dy/dx) = v·A·d²y/dx² + (v·A)’·dy/dx
Rearranging:
v·A·d²y/dx² = d/dx(v·A·dy/dx) - (v·A)’·dy/dx
This substitution allows us to replace the second derivative with first derivatives, reducing the continuity requirements on our approximation functions.
After integration by parts:
\(\displaystyle \int\limits_{0}^{L} A{\left(x \right)} v{\left(x \right)} \frac{d^{2}}{d x^{2}} y{\left(x \right)}\, dx = - A{\left(0 \right)} v{\left(0 \right)} \left. \frac{d}{d x} y{\left(x \right)} \right|_{\substack{ x=0 }} + A{\left(L \right)} v{\left(L \right)} \frac{d}{d L} y{\left(L \right)} - \int\limits_{0}^{L} \left(A{\left(x \right)} \frac{d}{d x} v{\left(x \right)} + v{\left(x \right)} \frac{d}{d x} A{\left(x \right)}\right) \frac{d}{d x} y{\left(x \right)}\, dx\)
Step 5: Apply boundary conditions v(0) = v(L) = 0:
EXPLANATION:
We choose our test functions v(x) to vanish at the domain boundaries (v(0) = v(L) = 0).
This is a crucial step that eliminates the boundary terms from integration by parts. By carefully selecting test functions that satisfy these conditions, we naturally incorporate the essential (Dirichlet) boundary conditions into our formulation.
This is one of the elegant aspects of the finite element method.
The boundary term [v·A·dy/dx]₀ᴸ vanishes because v(0) = v(L) = 0.
Step 6: Assemble the final weak form:
EXPLANATION: Now we collect all terms to form the complete weak formulation.
This form requires lower continuity requirements on our solution (only first derivatives appear).
The weak form is equivalent to the strong form for sufficiently smooth solutions, but it allows us to find approximate solutions with less smoothness.
This is the foundation for the finite element approximation where we’ll represent the solution as a combination of piecewise polynomial functions.
\(\displaystyle - \int\limits_{0}^{L} \left(A{\left(x \right)} \frac{d}{d x} v{\left(x \right)} + v{\left(x \right)} \frac{d}{d x} A{\left(x \right)}\right) \frac{d}{d x} y{\left(x \right)}\, dx + \int\limits_{0}^{L} B{\left(x \right)} v{\left(x \right)} \frac{d}{d x} y{\left(x \right)}\, dx + \int\limits_{0}^{L} C{\left(x \right)} v{\left(x \right)} y{\left(x \right)}\, dx = \int\limits_{0}^{L} F{\left(x \right)} v{\left(x \right)}\, dx\)
======================================================================
DISCRETIZATION USING SHAPE FUNCTIONS
======================================================================
Step 7: Define shape functions for N elements:
EXPLANATION:
Now we discretize the problem by representing both the solution y(x) and test function v(x) as linear combinations of shape functions.
This transforms our continuous problem into a discrete one with a finite number of unknowns (the coefficients).
The shape functions are typically chosen to be simple piecewise polynomials (often linear) that are non-zero only over a small portion of the domain .
This ‘local support’ property often leads to sparse matrices in the final system.
\(\displaystyle \left[ \phi_{1}{\left(x \right)}, \ \phi_{2}{\left(x \right)}, \ \phi_{3}{\left(x \right)}, \ \phi_{4}{\left(x \right)}, \ \phi_{5}{\left(x \right)}\right]\)
Step 8: Express solution and test functions using shape functions:
EXPLANATION:
We approximate both the solution y(x) and test function v(x) using the same set of shape functions.
The solution is written as a sum of shape functions with unknown coefficients a_j. Similarly, the test function is written as a sum with coefficients b_i.
This is the Galerkin approach, where we use the same functions for both approximation and testing.
\(\displaystyle y{\left(x \right)} = a_{1} \phi_{1}{\left(x \right)} + a_{2} \phi_{2}{\left(x \right)} + a_{3} \phi_{3}{\left(x \right)} + a_{4} \phi_{4}{\left(x \right)} + a_{5} \phi_{5}{\left(x \right)}\)
\(\displaystyle v{\left(x \right)} = b_{1} \phi_{1}{\left(x \right)} + b_{2} \phi_{2}{\left(x \right)} + b_{3} \phi_{3}{\left(x \right)} + b_{4} \phi_{4}{\left(x \right)} + b_{5} \phi_{5}{\left(x \right)}\)
Where:
NOTE: In the finite element method, we often choose shape functions that are equal to 1 at their corresponding node and 0 at all other nodes.
This makes it easy to enforce boundary conditions and interpret the coefficients aⱼ as the solution values at the nodes.
======================================================================
ASSEMBLING THE FEM SYSTEM
======================================================================
** Step 9: Substitute discretized functions into the weak form:**
EXPLANATION: We now substitute our approximations for y(x) and v(x) into the weak form.
This transforms the continuous weak form into a discrete algebraic system of equations.
Since the test function coefficients bᵢ are arbitrary, we can collect terms and set up a system of equations for the unknown coefficients aⱼ.
Step 10: Assemble the system matrix and right-hand side vector:
EXPLANATION:
When we substitute the discretized functions and collect terms, we obtain a linear system of equations K·a = F, where:
Each entry K[i,j] is computed by evaluating an integral that involves the shape functions phi_i and phi_j and their derivatives.
Similarly, each entry F[i] comes from an integral involving phi_i and the force function F(x).
For a simplified case where A(x) = B(x) = C(x) = 1, the entries would be:
\(\displaystyle K_{ij} = \int_0^L \frac{d\phi_i}{dx} \frac{d\phi_j}{dx} dx + \int_0^L \phi_i(x) \frac{d\phi_j}{dx} dx + \int_0^L \phi_i(x) \phi_j(x) dx\)
\(\displaystyle F_i = \int_0^L \phi_i(x) F(x) dx\)
NOTE: In practice, the assembly process is typically done by looping over elements,computing the element contributions, and adding them to the global system.
This element-by-element approach is more efficient and aligns with the local support property of shape functions.
======================================================================
SOLVING THE FEM SYSTEM
======================================================================
The final step would be solving the linear system:
\(\displaystyle \text{True}\)
After solving for the unknown coefficients a_j, we can reconstruct the solution y(x)
The following plots show the FEM solution for increasing numbers of elements, demonstrating how the numerical solution converges to the exact solution as the mesh is refined.
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 2, figsize=(10, 8))
x_exact = np.linspace(0, 1, 100)
y_exact = 0.5 * x_exact * (1 - x_exact)
for idx, n_elem in enumerate([2, 5, 10, 20]):
ax = axes[idx // 2, idx % 2]
# --- Compute FEM solution ---
n_nodes = n_elem + 1
h = 1.0 / n_elem
K = np.zeros((n_nodes, n_nodes))
F_vec = np.zeros(n_nodes)
for e in range(n_elem):
k_e = np.array([[1, -1], [-1, 1]]) / h
f_e = np.array([h / 2, h / 2])
K[e:e+2, e:e+2] += k_e
F_vec[e:e+2] += f_e
a = np.linalg.solve(K[1:-1, 1:-1], F_vec[1:-1])
y_fem = np.zeros(n_nodes)
y_fem[1:-1] = a
x_fem = np.linspace(0, 1, n_nodes)
# --- Plot on the subplot ---
ax.plot(x_exact, y_exact, 'b-', label='Exact: y = 0.5x(1-x)')
ax.plot(x_fem, y_fem, 'ro-', label=f'FEM ({n_elem} elements)')
ax.set_title(f'{n_elem} elements')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid(True)
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()