B-spline Finite Element Method (FEM)

The B-spline Finite Element Method (FEM) is a numerical technique that integrates B-spline basis functions—known for their smoothness and local support—into the finite element framework for the solution of differential equations. By employing exponential, polynomial, or trigonometric variants, the method achieves high-order accuracy, stability, and computational efficiency across a range of applications including nonlinear PDEs, geometric modeling, electronic structure, and atmospheric simulations.

1. Mathematical Foundations and B-spline Basis Construction

In B-spline FEM, the solution field is represented as a linear combination of B-spline basis functions. The generic form for a 1D function is

$U_N(x, t) = \sum_{j=-1}^{N+1} \delta_j(t) B_j(x),$

where $B_j(x)$ are B-spline basis functions (e.g., exponential cubic, polynomial, or trigonometric cubic), and $\delta_j(t)$ are the time-dependent coefficients.

B-spline functions are defined recursively and possess the following key properties:

• Local support: Each B-spline is nonzero only in a limited region, which yields sparse system matrices.
• High continuity: For degree $p$, B-splines are typically $C^{p-1}$-continuous except at repeated knots; this enables the accurate approximation of higher derivatives essential in, for example, biharmonic or bending problems.
• Partition of unity and flexibility: The family can be adapted via knot sequences or shape parameters (such as exponential or extension parameters) for optimally fitting function spaces or boundary layers.

The exponential cubic B-spline introduces a shape parameter $p$, enhancing flexibility. Trigonometric cubic B-splines provide spectral-like properties for oscillatory problems.

2. Temporal and Nonlinear Discretization: Crank-Nicolson and Linearization

Time-dependent PDEs are typically discretized using the Crank-Nicolson method, which is unconditionally stable and second-order accurate: $$\frac{U^{n+1}_i - U^n_i}{\Delta t} = \frac12\left[\mathcal{L}(U^{n+1}_i) + \mathcal{L}(U^n_i)\right],$$ where $\mathcal{L}$ is the spatial differential operator (including nonlinearities). Nonlinear terms are linearized at each time step using semi-implicit Rubin-Graves type schemes. For example, for a quadratic nonlinearity: $$(UU_x)^{n+1} \approx U^{n+1} U^n_x + U^{n} U^{n+1}_x - U^{n} U^n_x.$$ This ensures a linear system to solve at every time step while maintaining second-order accuracy.

3. Assembly and Solution of the Discrete System

Upon spatial and temporal discretization, the system at each time step reduces to a linear algebraic system due to the linearization of nonlinearities. For the coupled Burgers’ equation, the discrete system is written as: $A x^{n+1} = B x^n + b,$ where $x^{n+1}$ concatenates all the unknown spline coefficients at the new time, and $A, B$ are sparse, banded matrices constructed from the collocation of the basis and its derivatives. Boundary conditions (typically Dirichlet types) are enforced by modifying these matrices.

Because of the local support of B-splines, the resulting system is sparse and often banded, enabling efficient solution with algorithms such as the Thomas algorithm or other banded solvers.

4. Practical Advantages and Comparison with Other Methods

The B-spline FEM exhibits several key advantages:

• Higher accuracy with fewer nodes: Optimal parameter selection (e.g., exponential parameter $p$) allows exponential B-splines to accurately fit boundary and internal layers, outperforming traditional cubic B-splines, finite differences, and even some spectral methods, as demonstrated by error norms ($L_\infty$) that are an order of magnitude smaller for the same mesh.
• Flexibility and adaptivity: Free parameters in the spline basis and possible extension to higher-order splines accommodate problem-specific features, including steep gradients and oscillatory behavior.
• Sparse computational structures: Local support yields banded systems, reducing computational complexity.
• Unconditional stability: Crank-Nicolson discretization ensures stability for arbitrary time step sizes, as confirmed by von Neumann analysis.
• Simplicity and extensibility: The collocation approach, in contrast to Galerkin, requires only the evaluation of the basis and its derivatives at nodes, avoiding element-wise integrations.
• Robustness for nonlinear and coupled systems: The linearization procedure transforms even nonlinear, coupled PDEs into efficiently solvable systems.

Comparison with other spline-based and spectral methods:

5. Initialization and Boundary Conditions

Initial coefficients are determined by projecting initial data onto the B-spline basis using collocation at mesh nodes, potentially including derivative matching to conform with the PDE's order and boundary conditions.

Boundary values are incorporated directly into the system—either by modification of the system matrix and right-hand side vector or by substituting known values for coefficients associated with boundary basis functions.

6. Applications and Numerical Performance

The exponential cubic B-spline FEM has been validated on a variety of benchmark problems for the nonlinear coupled Burgers’ equation:

• Test Problems: Three classes of initial-boundary value problems, including sharp-layer and interacting front scenarios.
• Performance: Demonstrates second-order convergence in both space and time, with errors in the $L_\infty$ norm much smaller than those of classical methods for equivalent mesh sizes.
• Special features: Particularly accurate for challenging cases with complex initial data, steep layers, or boundary/interior sharp fronts.

7. Integration of the Crank-Nicolson Scheme

The synergy of B-spline spatial discretization with the Crank-Nicolson temporal scheme ensures that the stability properties are preserved while exploiting the high spatial accuracy, leading to globally accurate and robust time evolution. The method supports large time steps without loss of accuracy or the need for special stabilization mechanisms.

Key Formulae for Reference

• Solution expansion: $U_N(x, t) = \sum_{j=-1}^{N+1} \delta_j(t) B_j(x)$
• Crank-Nicolson time-stepping: $$\frac{U^{n+1}_i - U^n_i}{\Delta t} = \frac12 [LU^{n+1}_i + LU^n_i]$$
• Linear system (each time): $A x^{n+1} = B x^n + b$
• Error norm: $$L_\infty = \max_i | U_{\text{exact}}(x_i, t) - U_{\text{num}}(x_i, t) |$$

The exponential cubic B-spline finite element method, especially when combined with Crank-Nicolson time discretization, provides an accurate, efficient, and flexible approach for nonlinear coupled PDEs such as the coupled Burgers’ equation. It achieves superior accuracy for a wide range of physical and mathematical scenarios by leveraging the analytical flexibility of the B-spline basis and stable, high-order temporal integration.