error_analysis.py

# Imports
import numpy as np
import matplotlib.pyplot as plt
from solve_fem import solve_poisson, exact_solution, f0, f1, f2, f3, f4
from typing import Callable, Tuple


def compute_error(u_h: np.ndarray, x: np.ndarray, u_exact: Callable, du_exact: Callable=None, norm: str='L2') -> float:
    """
    Computes the L2 or H1 error between numerical solution u_h and exact solution u_exact using Gauss-Legendre 3-point quadrature.
    """

    N = (len(x) - 1) // 2 # Number of elements

    # Initialize squared errors
    L2_sq = 0.0
    H1_semi_sq = 0.0

    # 3-point Gauss-Legendre quadrature on reference interval [-1, 1]
    xi_ref, w_ref = np.polynomial.legendre.leggauss(3)
    # Values are:
    # xi_ref = np.array([-np.sqrt(3/5), 0.0, np.sqrt(3/5)])
    # w_ref  = np.array([5/9, 8/9, 5/9])

    # Shape functions
    def phi1(xi): return xi * (xi - 1) / 2
    def phi2(xi): return 1 - xi**2
    def phi3(xi): return xi * (xi + 1) / 2

    # Derivatives of shape functions
    def dphi1(xi): return xi - 0.5
    def dphi2(xi): return -2 * xi
    def dphi3(xi): return xi + 0.5

    for i in range(N):

        # Left, middle, and right nodes for i-th element
        i0, i1, i2 = 2*i, 2*i + 1, 2*i + 2

        # Left and right physical coordinates
        xL, xR = x[i0], x[i2]

        # Element size
        h = xR - xL

        # Jacobian for the transformation
        jac = h / 2

        # Extracting FEM solution values at the nodes
        u_vals = u_h[[i0, i1, i2]]

        # Loop over quadrature points
        for j in range(3):

            xi_j = xi_ref[j] # Reference Gauss point
            w_j = w_ref[j] # Corresponding weight

            # Map to physical coordinates
            xj = ((1 - xi_j) * xL + (1 + xi_j) * xR) / 2

            # Evaluate shape functions at xi_j
            phi_vals = np.array([phi1(xi_j), phi2(xi_j), phi3(xi_j)])

            # Numerical solution at xj
            uh_j = np.dot(u_vals, phi_vals)

            # Exact solution at xj
            uex_j = u_exact(xj)

            # Compute L2-norm error
            L2_sq += w_j * (uh_j - uex_j)**2 * jac

            if norm == 'H1':

                if du_exact is None:
                    raise ValueError("du_exact must be provided for H1 norm computation.")
                
                # Evaluate shape function derivatives at xi_j
                dphi_vals = np.array([dphi1(xi_j), dphi2(xi_j), dphi3(xi_j)])

                # Numerical derivative at xj (chain rule with dxi/dx = 2/h)
                duh_j = np.dot(u_vals, dphi_vals) * (2 / h)

                # Exact derivative at xj
                duex_j = du_exact(xj)

                # Compute H1-semi-norm error
                H1_semi_sq += w_j * (duh_j - duex_j)**2 * jac

    if norm == 'L2':
        return np.sqrt(L2_sq)
    elif norm == 'H1':
        return np.sqrt(L2_sq + H1_semi_sq)
    else:
        raise ValueError("Norm must be either 'L2' or 'H1'")


def get_exact_and_derivative(f: Callable) -> Tuple[Callable, Callable]:
    """
    Returns exact solution and its derivative corresponding to the given RHS function f.
    """

    # Get exact solution
    u = lambda x: exact_solution(x, f)

    # Get derivative of the exact solution
    if f == f0:
        du = lambda x: 0.5 - x
    elif f == f1:
        du = lambda x: np.pi * np.cos(np.pi * x)
    elif f == f2:
        du = lambda x: np.cos(np.pi * x) / np.pi
    elif f == f3:
        du = lambda x: -(x**2 / 2 - x**3 / 3) + 1 / 12
    elif f == f4:
        du = lambda x: -np.exp(x) + (np.e - 1)
    else:
        raise ValueError("Exact solution and derivative not implemented for this function")
    
    return u, du

def convergence_analysis(f: Callable) -> None:
    """
    Perform convergence analysis and plot the results.
    """

    # Get exact solution and its derivative
    u_exact, du_exact = get_exact_and_derivative(f)

    # Array with different number of elements
    N_vals = np.array([10, 20, 40, 80, 160, 320, 640], dtype=int)

    # Initialize error arrays
    L2_errors = np.zeros_like(N_vals, dtype=np.float64)
    H1_errors = np.zeros_like(N_vals, dtype=np.float64)

    # Compute errors for different step sizes
    for i, N in enumerate(N_vals):
        x, u_num = solve_poisson(N, f)
        L2_errors[i] = compute_error(u_num, x, u_exact, norm='L2')
        H1_errors[i] = compute_error(u_num, x, u_exact, du_exact, norm='H1')

    # Compute convergence rates
    order_L2 = np.polyfit(np.log10(1/N_vals), np.log10(L2_errors), 1)[0]
    order_H1 = np.polyfit(np.log10(1/N_vals), np.log10(H1_errors), 1)[0]

    # Plotting
    fig, axs = plt.subplots(1, 2, figsize=(14, 7))
    
    axs[0].loglog(1/N_vals, L2_errors, 'o-', label=f"$L^2$ Order ≈ {order_L2:.2f}")
    axs[0].set_xlabel("h (step size)", fontsize=12)
    axs[0].set_ylabel("$L^2$ Error", fontsize=12)
    axs[0].set_title("Convergence in $L^2$-norm", fontsize=14)
    axs[0].grid(True, which="both", linestyle="--")
    axs[0].legend()

    axs[1].loglog(1/N_vals, H1_errors, 'o-', color='crimson', label=f"$H^1$ Order ≈ {order_H1:.2f}")
    axs[1].set_xlabel("h (step size)", fontsize=12)
    axs[1].set_ylabel("$H^1$ Error", fontsize=12)
    axs[1].set_title("Convergence in $H^1$-norm", fontsize=14)
    axs[1].grid(True, which="both", linestyle="--")
    axs[1].legend()

    plt.tight_layout()
    plt.savefig("convergence.svg", format="svg")
    plt.show()

    # Print convergence rates
    print(f"Estimated $L^2$ order of convergence: {order_L2:.2f}")
    print(f"Estimated $H^1$ order of convergence: {order_H1:.2f}")

#convergence_analysis(f0)
convergence_analysis(f1)
#convergence_analysis(f2)
#convergence_analysis(f3)
#convergence_analysis(f4)