Presentation on Siam Applied book club(ABC)

Abstract

This session, we will focus on a general question regarding the finite element method for handling the special differential equation, giving some basic ideas and definitions. We also may look at some codes. The topic will cover the first 3 sections of this slideshow.

Mathematical Formulation

The 1D second-order elliptic equation we solve is:

For Example 1: $c(x) = 1$, $f(x) = 2$, $a = 0$, $b = 1$, $g_a = 0$, $g_b = 1$

The weak formulation leads to the linear system $Au = b$ where:

• Stiffness matrix: $A_{ij} = \int_a^b c(x) \phi_i'(x) \phi_j'(x) dx$
• Load vector: $b_i = \int_a^b f(x) \phi_i(x) dx$

Implementation

Stiffness Matrix Generation

def generate_stiff_matrix_by_linear(node_sum, mesh_size, c_func=None):
    """
    Generate stiffness matrix for 1D finite element method
    
    For the equation -d/dx(c(x) * du/dx) = f(x), the stiffness matrix is:
    A[i,j] = ∫[0,1] c(x) * φ'_i(x) * φ'_j(x) dx
    """
    N = node_sum
    h = mesh_size
    
    # Default c(x) = 1 if no function provided
    if c_func is None:
        c_func = lambda x: 1.0
    
    # Initialize the matrix
    A = np.zeros((N + 1, N + 1))
    
    # Process each element [x_i, x_{i+1}] for i = 0, 1, ..., N-1
    for i in range(N):  # i = 0 to N-1 (N elements)
        x_i = i * h
        x_i_plus_1 = (i + 1) * h
        
        # Compute integral of c(x) over element [x_i, x_{i+1}]
        c_integral = compute_integral(c_func, x_i, x_i_plus_1)
        element_contribution = c_integral / h
        
        # Add contribution to diagonal terms
        A[i, i] += element_contribution
        A[i + 1, i + 1] += element_contribution
        
        # Add contribution to off-diagonal terms
        A[i, i + 1] -= element_contribution
        A[i + 1, i] -= element_contribution
    
    # Scale the matrix by 1/h for the 1D Poisson equation
    A = A / h
    
    return A

Load Vector Generation

def generate_load_vector_by_linear(node_sum, mesh_size, f_func=None):
    """
    Generate load vector for 1D finite element method
    
    For the equation -d/dx(c(x) * du/dx) = f(x), the load vector is:
    b[i] = ∫[0,1] f(x) * φ_i(x) dx
    """
    N = node_sum
    h = mesh_size
    
    # Default f(x) = 1 if no function provided
    if f_func is None:
        f_func = lambda x: 1.0
    
    # Initialize the load vector
    b = np.zeros((N + 1, 1))
    
    # Process each element [x_i, x_{i+1}] for i = 0, 1, ..., N-1
    for i in range(N):  # i = 0 to N-1 (N elements)
        x_i = i * h
        x_i_plus_1 = (i + 1) * h
        
        # For element [x_i, x_{i+1}], the basis functions are:
        # φ_i(x) = (x_{i+1} - x) / h
        # φ_{i+1}(x) = (x - x_i) / h
        
        # Contribution to b[i] from element [x_i, x_{i+1}]
        integrand_i = lambda x: f_func(x) * (x_i_plus_1 - x) / h
        b[i] += compute_integral(integrand_i, x_i, x_i_plus_1)
        
        # Contribution to b[i+1] from element [x_i, x_{i+1}]
        integrand_i_plus_1 = lambda x: f_func(x) * (x - x_i) / h
        b[i + 1] += compute_integral(integrand_i_plus_1, x_i, x_i_plus_1)
    
    return b

Example 1: Simple Constant Coefficient Case

We solve a simple case with constant coefficients:

For this example: $c(x) = 1$, $f(x) = 2$, $a = 0$, $b = 1$, $g_a = 0$, $g_b = 1$

The analytical solution of this problem is $u(x) = -x^2 + 2x$, which can be verified by direct substitution.

Results

Here are the results from our FEM implementation:

Finite Element Method results showing solution comparison, error distribution, convergence study, and stiffness matrix visualization

Performance Summary

• L2 Error: 3.79 × 10⁻¹⁶ (machine precision)
• L∞ Error: 6.66 × 10⁻¹⁶ (machine precision)
• Mesh size: h = 0.1 (N = 10 internal nodes)
• Analytical solution: u(x) = -x² + 2x
• Result: Essentially exact solution due to quadratic exactness of linear FEM

This example demonstrates that for problems where the analytical solution is a quadratic polynomial, the linear finite element method achieves machine precision accuracy, as the method is exact for polynomials of degree up to 1, and the integration is exact for quadratic functions.

Example 2: Numerical Example with Variable Coefficient

We solve the following 1D second-order elliptic equation using the linear finite element method:

The analytical solution of this problem is $u(x) = x\cos(x)$, which can be used to compute the error of the numerical solution.

Implementation Steps

1. Construct P Matrix (Node Coordinates)

def construct_P_matrix_linear(a, b, N):
    """
    Construct P matrix (node coordinates matrix) for linear finite elements
    The j-th column stores the coordinates of the j-th finite element node.
    """
    h = (b - a) / N
    x_nodes = np.array([a + j * h for j in range(N + 1)])
    
    # For 1D linear elements, P is (2, N+1) where each column is [x_j, x_j]
    P = np.zeros((2, N + 1))
    P[0, :] = x_nodes
    P[1, :] = x_nodes
    
    return P

2. Construct T Matrix (Element Connectivity)

def construct_T_matrix_linear(N):
    """
    Construct T matrix (element connectivity matrix) for linear finite elements
    The n-th column stores the global node indices of the n-th mesh element.
    """
    # For 1D linear elements, T is (2, N) where n-th column is [n, n+1]^T
    T = np.zeros((2, N), dtype=int)
    for n in range(N):
        T[0, n] = n      # First node of element n
        T[1, n] = n + 1  # Second node of element n
    
    return T

3. Construct Stiffness Matrix A (Algorithm IV)

def generate_stiff_matrix_by_assembly(P, T, c_func=None):
    """
    Generate stiffness matrix using Algorithm IV (local assembly)
    
    Algorithm IV:
    FOR n = 1, ..., N:
        FOR α = 1, ..., Nlb:
            FOR β = 1, ..., Nlb:
                Compute r = ∫_{x_n}^{x_{n+1}} c ψ'_{nα} ψ'_{nβ} dx
                Add r to A(Tb(β, n), Tb(α, n))
    """
    Nb = P.shape[1]  # Number of global basis functions (nodes)
    N = T.shape[1]   # Number of elements
    Nlb = 2          # Number of local basis functions for linear elements
    
    if c_func is None:
        c_func = lambda x: 1.0
    
    A = np.zeros((Nb, Nb))
    
    for n in range(N):
        i_global_1 = T[0, n]
        i_global_2 = T[1, n]
        x_n = P[0, i_global_1]
        x_n_plus_1 = P[0, i_global_2]
        h = x_n_plus_1 - x_n
        
        # Local basis derivatives: ψ'_{n1} = -1/h, ψ'_{n2} = 1/h
        psi_prime_1 = -1.0 / h
        psi_prime_2 = 1.0 / h
        
        for alpha in range(1, Nlb + 1):
            for beta in range(1, Nlb + 1):
                psi_prime_alpha = psi_prime_1 if alpha == 1 else psi_prime_2
                psi_prime_beta = psi_prime_1 if beta == 1 else psi_prime_2
                
                # Compute r = ∫_{x_n}^{x_{n+1}} c(x) * ψ'_{nα} * ψ'_{nβ} dx
                c_integral = compute_integral(c_func, x_n, x_n_plus_1)
                r = psi_prime_alpha * psi_prime_beta * c_integral
                
                # Map to global indices
                i_global_beta = T[beta - 1, n]
                i_global_alpha = T[alpha - 1, n]
                
                # Add r to A(Tb(β, n), Tb(α, n))
                A[i_global_beta, i_global_alpha] += r
    
    return A

4. Construct Load Vector b (Algorithm V)

def generate_load_vector_by_assembly(P, T, f_func=None):
    """
    Generate load vector using Algorithm V (local assembly)
    
    Algorithm V:
    b = zeros(N_b, 1)
    FOR n = 1, ..., N:
        FOR β = 1, ..., N_lb:
            Compute r = ∫_{x_n}^{x_{n+1}} f * ψ_{nβ} dx
            b(T_b(β, n), 1) = b(T_b(β, n), 1) + r
    """
    Nb = P.shape[1]
    N = T.shape[1]
    Nlb = 2
    
    if f_func is None:
        f_func = lambda x: 1.0
    
    b = np.zeros((Nb, 1))
    
    for n in range(N):
        i_global_1 = T[0, n]
        i_global_2 = T[1, n]
        x_n = P[0, i_global_1]
        x_n_plus_1 = P[0, i_global_2]
        h = x_n_plus_1 - x_n
        
        for beta in range(1, Nlb + 1):
            # Compute r = ∫_{x_n}^{x_{n+1}} f(x) * ψ_{nβ}(x) dx
            def integrand(x):
                if beta == 1:
                    psi_beta = (x_n_plus_1 - x) / h
                else:
                    psi_beta = (x - x_n) / h
                return f_func(x) * psi_beta
            
            r = compute_integral(integrand, x_n, x_n_plus_1)
            
            # Map to global index
            i_global_beta = T[beta - 1, n]
            b[i_global_beta, 0] += r
    
    return b

5. Apply Dirichlet Boundary Conditions (Algorithm VI)

def apply_dirichlet_boundary_conditions(A, b, P, boundary_nodes, g_func):
    """
    Apply Dirichlet boundary conditions using Algorithm VI
    
    Algorithm VI:
    FOR k = 1, ..., nbn:
        IF boundarynodes(1, k) shows Dirichlet condition, THEN
            i = boundarynodes(2, k)
            A(i, :) = 0
            A(i, i) = 1
            b(i) = g(Pb(i))
        ENDIF
    END
    """
    A_mod = A.copy()
    b_mod = b.copy()
    nbn = boundary_nodes.shape[1]
    
    for k in range(nbn):
        condition_type = boundary_nodes[0, k]
        if condition_type == 1:  # Dirichlet condition
            i = int(boundary_nodes[1, k])
            A_mod[i, :] = 0
            A_mod[i, i] = 1
            x_i = P[0, i]
            b_mod[i, 0] = g_func(x_i)
    
    return A_mod, b_mod

Results

Using the assembly-based implementation (Algorithms IV, V, and VI), we obtain the following results:

FEM assembly results showing solution comparison, error distribution, convergence study, stiffness matrix, node coordinates, and element connectivity

Maximum Numerical Errors at All Mesh Nodes

The following table shows the maximum absolute error at all mesh nodes for different mesh sizes $h$:

Convergence Analysis

The convergence rate check shows that when the mesh size $h$ is halved, the error decreases by approximately a factor of 4, confirming $O(h^2)$ convergence:

• Error($h=1/4$) / Error($h=1/8$) = 4.0022
• Error($h=1/8$) / Error($h=1/16$) = 3.9822
• Error($h=1/16$) / Error($h=1/32$) = 3.9930
• Error($h=1/32$) / Error($h=1/64$) = 3.9995
• Error($h=1/64$) / Error($h=1/128$) = 3.9993

This demonstrates the expected quadratic convergence rate for linear finite elements, where the error is proportional to $h^2$.