A Lattice Boltzmann BGK Model with an Amending Function for Two-Dimensional Second-Order Nonlinear Partial Differential Equations

TL;DR

This paper introduces a new lattice Boltzmann model for solving complex nonlinear equations, showing it is efficient and accurate.

Contribution

The novel contribution is a BGK lattice Boltzmann model with an amending function for second-order nonlinear PDEs.

Findings

The model successfully recovers macroscopic equations via Chapman–Enskog analysis.

Numerical experiments show excellent agreement with exact solutions.

The model is robust for strongly nonlinear systems.

Abstract

A mesoscopic lattice Boltzmann method based on the BGK model is proposed to solve a class of two-dimensional second-order nonlinear partial differential equations by incorporating an amending function. The model provides an efficient and stable framework for simulating initial value problems of second-order nonlinear partial differential equations and is adaptable to various nonlinear systems, including strongly nonlinear cases. The numerical characteristics and evolution patterns of these nonlinear equations are systematically investigated. A D2Q4 lattice model is employed, and the kinetic moment constraints for both local equilibrium and correction distribution functions are derived in the four velocity directions. Explicit analytical expressions for these distribution functions are presented. The model is verified to recover the target macroscopic equations in the continuous limit via Chapman–Enskog analysis. Numerical experiments using exact solutions are performed to assess the model’s accuracy and stability. The results show excellent agreement with exact solutions and demonstrate the model’s robustness in capturing nonlinear dynamics.

Figures

Full text

1. Introduction

The solution of nonlinear partial differential Equations (PDEs) constitutes a core area of research in applied mathematics and computational science, as such equations frequently arise in the modeling of complex physical, biological, and engineering systems [1,2]. Nonlinear PDEs are widely used to describe phenomena such as fluid dynamics, heat transfer, wave propagation, chemical reactions, biological pattern formation, and financial modeling [3,4,5,6]. Due to the intrinsic nonlinearity, these equations often exhibit rich and diverse behaviors, including shock formation, wave interactions, and pattern emergence, which render exact solutions intractable in most cases [2,3]. Consequently, considerable efforts have been devoted to the development of robust and efficient numerical methods, such as finite difference, finite element, spectral, and kinetic-based schemes, for approximating solutions and exploring the qualitative dynamics of nonlinear PDEs under various initial and boundary conditions [7,8]. Due to the fact that nonlinear PDEs only have exact solutions under special initial and boundary conditions and lack a unified expression, it is therefore necessary to vigorously develop efficient and high-precision numerical simulation methods.

In recent years, kinetic theory-based methods have witnessed rapid development. Originally devised to simulate fluid flows, approaches such as the lattice Boltzmann method (LBM) [9,10,11] and the discrete Boltzmann method (DBM) [12,13,14,15,16,17,18,19,20] have also proven effective in solving nonlinear problems. Among them, the LBM has undergone significant advancements, leading to the development of various models tailored to different classes of nonlinear PDEs. By modeling the evolution of particle distribution functions on a discrete lattice, LBM captures nonlinear macroscopic behaviors through local interactions and collision rules. Its inherent simplicity, parallelizability, and flexibility make it well-suited for simulating nonlinear systems with complex geometries and boundary conditions [21]. Compared with the conventional numerical methods, the LBM enjoys many notable advantages, such as geometrical flexibility, numerical efficiency and ease in incorporating complex boundary conditions. Furthermore, it can be naturally adapted to parallel processes computing [22,23].

Many researchers attempted to extend the LBM to solve nonlinear PDEs. Wei et al. (2018) [24,25] introduced a two-dimensional lattice Boltzmann model for thermal incompressible flows. By modifying the equilibrium distribution function and coupling a discrete force, the model improves stability and accuracy. Wang (2020) developed a lattice Boltzmann model to simulate $[eqn]$ D solitary waves in dusty plasma governed by the extended Zakharov–Kuznetsov equation, confirming the model’s validity through numerical experiments [26]. Gorakifard et al. (2022) [27] proposed the meshless local Petrov–Galerkin cumulant LBM to enhance aeroacoustic analysis by improving stability at low viscosities and computational efficiency. Their method integrates cumulant collision modeling and mesh-free streaming, demonstrating superior accuracy and reduced runtime compared to conventional approaches. Boghosian et al. (2024) [28] explore numerical approximations of LBM for inhomogeneous advection in one dimension. Their work innovatively derives high-order equivalent partial differential equations, demonstrating the critical role of initialization in accuracy. Spectral methods validate asymptotic error behavior in steady and unsteady states. Chen et al. (2024) [29] presented a novel fourth-order multiple-relaxation-time LBM for two-dimensional diffusion equations, bridging mesoscopic and macroscopic approaches. By deriving a corresponding finite-difference scheme, they rigorously establish accuracy and unconditional stability, offering valuable theoretical insights and computational advancements in modeling diffusion phenomena. Chai et al. (2016) [30] developed a multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection–diffusion equations. Their MRT model recovers the correct macroscopic equations and achieves second-order spatial convergence. Compared to the BGK model, the MRT approach reduces numerical slip on boundaries and improves stability by tuning relaxation parameters. Wang et al. (2024) [31] applied the LBM to solve the generalized Zakharov equation, enhancing numerical accuracy and efficiency in plasma wave simulations. Their model demonstrates second-order convergence and superior computational performance compared to traditional methods, advancing nonlinear wave analysis in plasma physics. Wang et al. (2024) [32] applied the lattice Boltzmann method to study wave propagation governed by variable coefficient nonlinear Schrödinger equations, demonstrating advantages in accuracy and efficiency. Liang et al. (2024) [33] proposed a novel thermal lattice Boltzmann model for accurately simulating liquid-vapor phase changes. By avoiding complex gradient and Laplacian calculations, their approach enhances numerical precision while simplifying the modeling process. Chen et al. (2025) [34] proposed a fourth-order multiple-relaxation-time LBM for d-dimensional coupled Burgers’ equations using the Cole–Hopf transformation. This approach eliminates nonlinear convection terms, enhancing stability and accuracy. Chai et al. (2025) [35] proposed a novel mesoscopic regularized LBM-based macroscopic numerical scheme for the Allen–Cahn equation. Li et al. (2025) [36] proposed a novel LBM for a class of viscous wave equations by transforming the original equation to eliminate mixed third-order derivatives. Their method improves numerical stability and achieves second-order accuracy, validated through simulations. Du et al. (2025) constructed a lattice Boltzmann model with BGK operator for solving the fractional Laplacian Reynolds-averaged Navier–Stokes equations by leveraging the fractional centered difference scheme [37]. They also proposed a lattice Boltzmann model with a suitable equilibrium distribution function to solve the surface quasi-geostrophic equations with fractional Laplacian [38].

Apart from successful verification in general problems of fluid mechanics, it has elegant application in relevant domains, ranging from fluid turbulence, multi-phase flows, chemical reactions and hydromagnetics. Recent research has mainly focused on expanding its application areas and improving computational methods for LBM. In applications, Bocanegra et al. (2024) applied LBM to simulate sound wave propagation, demonstrating its excellent performance in handling complex geometric boundaries and the coupling of sound fields with mean flow [39]. Ginzburg (2025) used three-dimensional LBM to simulate the full-volume flow of Newtonian fluids in the colon, contributing to a better understanding of gastrointestinal physiology and pathology [40]. In computational methods, Maquart et al. (2022) developed an LBM version with deformable boundaries, using a Dirichlet velocity scheme to precisely describe surface deformation on regular computational grids without requiring traditional LBM adaptive boundary grids or surface remeshing [41]. This approach is crucial for simulating dynamic shape changes in biological tissues. Hosseini et al. (2024) innovatively proposed a new LBM model and numerical simulation strategy for combustion simulations, making it more stable and efficient in handling compressible flow and multi-component transport problems in combustion simulations [42]. Wawrzyniak et al. (2025) introduced an LBM-based quantum algorithm for solving convection–diffusion equations, providing a new method for solving Navier–Stokes equations on quantum computers in the future [43].

In this paper, we consider two-dimensional second-order nonlinear partial differential Equation (SONPDE) with initial condition. The basic form is shown as follows

$$[eqn]$$

where $[eqn]$ is a bounded region, $[eqn]$   $[eqn]$ , $[eqn]$ represents the wave displacement at position $[eqn]$ and time t. $[eqn]$ is the known differential functions of $[eqn]$ , used to model nonlinear convective effects. $[eqn]$ is the molecular diffusivity characterizing the Brownian motion, and $[eqn]$ are convection coefficients. $[eqn]$ and $[eqn]$ are known functions.

This equation can be considered as the linearized version of the two-dimensional Navier–Stokes equations in vorticity formulation with viscosity $[eqn]$ (Re is the Reynolds number). When $[eqn]$ , the equation becomes the classical convection–diffusion equation. The convection–diffusion equation is a fundamental model in fluid dynamics, heat transfer, and mass transport phenomena. It describes the combined effects of convection (transport due to fluid motion) and diffusion (spreading due to concentration gradients) on the evolution of a scalar quantity, such as temperature or concentration. The equation is widely used in various fields, including environmental science, chemical engineering, and materials science. This equation can be considered to be the linearized version of the two-dimensional Navier–Stokes equations in vorticity formulation with viscosity $[eqn]$ (Re is the Reynolds number). Although multiple-relaxation-time (MRT) and entropic collision operators have been shown to improve numerical stability and robustness for the full Navier–Stokes equations, particularly at high Reynolds numbers, our current work focuses on the BGK model applied to this linearized system as a foundational step. The BGK framework allows for a clear demonstration and validation of the proposed lattice Boltzmann method with an amending function. We acknowledge that incorporating MRT and entropic collision operators could further enhance stability and adaptability, and these extensions are planned for future work once the current model’s effectiveness is fully established. This staged approach ensures the methodological soundness and clarity before advancing to more complex collision models and nonlinear systems.

In order to improve computational efficiency and extend the application of LBM in solving nonlinear partial differential equations, we adopt a D2Q4 discrete velocity model to construct a lattice BGK Boltzmann model with an amending function for Equation (1). By using the Chapman–Enskog expansion [44], The governing evolution equation can be accurately derived from the continuous BGK Boltzmann equation. Comprehensive comparisons between numerical results and analytical solutions are conducted. The simulations demonstrate excellent agreement with the analytical solutions.

The content of this paper is arranged as follows. In Section 2, we introduce the lattice Boltzmann model and derive the governing equation. In Section 3, we present numerical experiments to validate the accuracy and reliability of the proposed model. Finally, we summarize our findings in Section 4.

2. Lattice Boltzmann Model

The lattice Boltzmann model employed in this study is the standard lattice Bhatnagar–Gross–Krook (LBGK) model based on the D2Q4 velocity set. The discrete velocity model is shown in Figure 1, and the D2Q4 lattice consists of four discrete velocities, each aligned with the Cartesian axes. The corresponding velocity vectors $[eqn]$ for $[eqn]$ are defined as follows:

$$[eqn]$$

The evolution of the particle distribution function $[eqn]$ with an amending source term $[eqn]$ is governed by the following LBGK equation:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are defined as distribution function and equilibrium distribution function, respectively, $[eqn]$ is an amending function. Here, c is the lattice characteristic speed relating spatial and temporal scales via the lattice CFL condition $[eqn]$ . The parameter $[eqn]$ is a small expansion parameter related to the Knudsen number, h is the spatial step size, and T is a characteristic timescale satisfying the diffusive scaling $[eqn]$ . The dimensionless relaxation time $[eqn]$ controls the rate of relaxation towards equilibrium and must satisfy the stability condition $[eqn]$ [45].

The macroscopic variable $[eqn]$ can be defined as the zeroth moment of the distribution function:

$$[eqn]$$

and the equilibrium distribution function $[eqn]$ should meet the following conservation laws:

$$[eqn]$$

Then, through choosing appropriate local equilibrium distribution, we can retrieve the corresponding macroscopic equation correctly.

Indeed, applying Taylor expansion to the left hand of Equation (2) and retaining terms up to $[eqn]$ , we can get

$$[eqn]$$

To derive the macroscopic equation Equation (1), the Chapman–Enskog expansion is applied:

$$[eqn]$$

The Knudsen number $[eqn]$ is defined as follows:

$$[eqn]$$

where ℓ is the mean free path, and L is the characteristic length, which can be taken as the time step $[eqn]$ . And $[eqn]$   $[eqn]$ represent higher-order corrections associated with nonequilibrium effects at different time scales, which satisfy the solvability conditions:

$$[eqn]$$

Substituting Equation (6) into Equation (5), we obtain a series of lattice Boltzmann equations in different timescales:

$[eqn]$ :

$$[eqn]$$

$[eqn]$ :

$$[eqn]$$

Substituting Equation (9) into Equation (10), we get

$$[eqn]$$

Summing the two sides of Equation (11) with respect to i and using Equations (3) and (8), we can obtain

$$[eqn]$$

We select the local equilibrium distribution functions $[eqn]$ to satisfy the following:

$$[eqn]$$

and define the amending functions $[eqn]$ as follows:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are some constants to be determined. The structure of the amending function $[eqn]$ ensures isotropy and consistency with the target convection term, while allowing for flexible tuning via the parameters $[eqn]$ and $[eqn]$ .

Then, we get

$$[eqn]$$

and

$$[eqn]$$

Substituting Equations (15) and (16) into Equation (12), we can obtain

$$[eqn]$$

To recover Equation (1), just let

$$[eqn]$$

In the computational process, we can let

$$[eqn]$$

then

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are known parameters of the problem.

With the above choice of parameters, the proposed lattice Boltzmann BGK model, equipped with an amending function, successfully recovers the general form of a two-dimensional nonlinear convection–diffusion equation:

$$[eqn]$$

This formulation offers a flexible and consistent framework for simulating a broad class of second-order nonlinear partial differential equations, with tunable convective and diffusive components. The incorporation of the amending function $[eqn]$ enables accurate recovery of macroscopic dynamics beyond the scope of standard equilibrium-based LBM approaches.

In summary, the proposed mesoscopic lattice Boltzmann BGK model, built upon the D2Q4 velocity lattice and augmented with analytically derived equilibrium and correction distribution functions, satisfies the required kinetic moment constraints and recovers the target macroscopic equation through Chapman–Enskog analysis in the diffusive limit. This theoretical foundation paves the way for the numerical experiments in the next section, where the model’s accuracy and robustness are systematically validated using benchmark problems with known exact solutions.

3. Numerical Simulation

In this section, to test the above model, numerical simulations of SONPDE are performed. The distribution function $[eqn]$ is initialized by setting to equal $[eqn]$ for all nodes at $[eqn]$ . And the macroscopic variable $[eqn]$ in Equation (1) is initialized by initial condition. The non-equilibrium extrapolation scheme proposed by Guo et al. (2002) is used for implementing boundary conditions [46]. To quantitatively assess the numerical accuracy of our proposed model, three different error metrics (MAXE, MAE, GRE) are defined, where $[eqn]$ and $[eqn]$ represent the numerical and exact solutions, respectively, with $[eqn]$ and $[eqn]$ denoting the number of grid points in the x and y directions. All summations are performed over the entire computational domain. In the following numerical simulation, unless otherwise specified, the evolution time t is defined as $[eqn]$ .

The maximum absolute error (MAXE) is defined as follows:

$$[eqn]$$

The mean absolute error (MAE) is defined as follows:

$$[eqn]$$

The global relative error (GRE) is defined as follows:

$$[eqn]$$

Example 1. Consider the following convection–diffusion equation in the region $[eqn]$ , $[eqn]$ :

$$[eqn]$$

the initial and boundary conditions are given below:

$$[eqn]$$

The corresponding exact solution for this instance is given as

$$[eqn]$$

where

$$[eqn]$$

In the simulation, we use $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ . The simulation results of Figure 2 demonstrate that the numerical solution behaves in a similar pattern to the analytical one across the domain builds confidence in the correctness and reliability of the numerical approach. It reveals the global behavior of the system, showing its structure and characteristics after long-term evolution. This visualization aids in understanding the dynamic behavior of the complex system.

Additionally, we present the two-dimensional visual comparison of profiles at $[eqn]$ at specific times in Figure 3. The two-dimensional profile comparison at $[eqn]$ at specific times in Figure 3 provides detailed local evolution information, further validating the accuracy of the numerical results against the exact solutions. The numerical results exhibit a good match with the exact solutions, confirming the accuracy of the method. The comparison shows that the numerical solution captures the essential features of the exact solution, including its spatial distribution and temporal evolution.

Figure 4 presents the contour plot of the numerical solution to the two-dimensional convection–diffusion equation. It is observed that the solution exhibits a smooth spatial distribution, with values ranging approximately from 0.005 to 0.025. The solution reaches relatively higher values near the central region, indicating the presence of a source term or boundary-driven high concentration area. A clear gradient along a principal direction, possibly from the lower left to the upper right, suggests a dominant convective flow in that direction. This behavior is typical in convection–diffusion problems, where the solution is advected along the flow direction while being smoothed by diffusion. Overall, the contour lines are well-behaved and uniformly distributed without numerical oscillations or non-physical artifacts. This indicates that the adopted numerical scheme possesses good stability and conservativeness, effectively capturing the combined effects of convection and diffusion.

Furthermore, the MAXE, MAE, and GRE for solutions of Example 1 at different times are summarized in Table 1. The table demonstrates the decreasing trend of errors as time progresses, indicating the stability and accuracy of the numerical method. Specifically, the MAXE reduces from $[eqn]$ at $[eqn]$ to $[eqn]$ at $[eqn]$ , showcasing a significant improvement in precision over time. Similarly, the MAE and GRE also exhibit consistent reductions, confirming the reliability of the proposed lattice Boltzmann model for solving the convection–diffusion equation. These results highlight the effectiveness of the numerical scheme in capturing the dynamics of the system with high accuracy.

To verify the spatial accuracy of the proposed lattice Boltzmann model, we use this example for simulation. The simulation is conducted using four different grid resolutions: $[eqn]$ , $[eqn]$ , and $[eqn]$ . The parameters are set as $[eqn]$ , $[eqn]$ , and the simulation is run until $[eqn]$ . The GRE and the corresponding convergence rate are presented in Table 2. The convergence rate is calculated using the following formula:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the MAXE corresponding to grid spacings $[eqn]$ and $[eqn]$ , respectively. As shown in Table 2, the GRE decreases as the grid is refined, and the convergence rate approaches 2, which is consistent with the theoretical second-order accuracy of the LBM.

Example 2. Consider the following the convection–diffusion equation in the region $[eqn]$ , $[eqn]$ with the initial and boundary conditions [47]:

$$[eqn]$$

The initial conditions are given below:

$$[eqn]$$

The corresponding exact solution for this instance is given that

$$[eqn]$$

The boundary conditions can be extracted from the exact solution.

In the simulation, we use $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ . The simulation results in Figure 5 show a strong agreement between the numerical and analytical solutions over the entire computational domain. The solution surfaces exhibit similar structural patterns, demonstrating that the numerical scheme accurately captures the essential features of the exact solution. This close match not only verifies the correctness of the proposed method but also reflects the method’s capability in resolving the long-term dynamic behavior of the system. The three-dimensional visualization further enhances our understanding of the system’s global evolution.

In addition, we present the two-dimensional visual comparison of profiles at $[eqn]$ at specific times in Figure 6. The comparison provides detailed insights into the local evolution of the solution, allowing for a more precise evaluation of the numerical method’s accuracy. The numerical profiles are in good agreement with the corresponding exact solutions, clearly reproducing both the spatial structure and the temporal development of the exact solution. This consistency further confirms that the proposed numerical scheme effectively captures the key dynamic features of the system.

Figure 7 presents the contour plot of the numerical solution to the two-dimensional convection–diffusion equation. The solution exhibits a smooth and continuous spatial distribution, with values ranging approximately from 0.02 to 0.16. A prominent high-concentration region is observed near the upper central part of the domain, likely corresponding to a localized source term or boundary-induced effect. The contour lines reveal a distinct gradient oriented from the lower left toward the upper right, indicative of a dominant convective flow in that direction. This directional behavior is consistent with the expected dynamics in convection–diffusion problems, where convection transports the quantity along a flow path while diffusion smooths out local variations. The uniformly spaced and well-behaved contour lines suggest that the employed numerical scheme is stable and conservative, effectively capturing the interplay between convective and diffusive processes without introducing non-physical artifacts.

Furthermore, the MAXE, MAE, and GRE for solutions of Example 2 at different times are summarized in Table 3. The table shows that the errors decrease as time progresses, indicating the stability and accuracy of the numerical method. Specifically, the MAXE reduces from $[eqn]$ at $[eqn]$ to $[eqn]$ at $[eqn]$ , showcasing a significant improvement in precision over time. Similarly, the MAE and GRE also exhibit consistent reductions, confirming the reliability of the proposed lattice Boltzmann model for solving the convection–diffusion equation. These results highlight the effectiveness of the numerical scheme in capturing the dynamics of the system with high accuracy.

To investigate the influence of the relaxation time $[eqn]$ on the numerical accuracy of the proposed LBM model, we conduct an additional test based on the setup of Example 2. The spatial discretization is $[eqn]$ , and the time step is fixed at $[eqn]$ . The relaxation time $[eqn]$ is varied by selecting different viscosity coefficients $[eqn]$ . All simulations are performed up to a final time of $[eqn]$ . The errors (MAXE, MAE, GRE) are computed and compared for these different $[eqn]$ values (and corresponding $[eqn]$ values), as shown in Table 4.

The results presented in Table 4 (with illustrative error values) and Figure 8 show how varying the viscosity coefficient $[eqn]$ , and consequently the relaxation time $[eqn]$ , affects numerical accuracy when other parameters like $[eqn]$ and h are fixed. In this setup, as $[eqn]$ increases, $[eqn]$ also increases (from 0.52 to 0.60 in the table). The illustrative error values reveal a clear trend that larger $[eqn]$ (and thus larger $[eqn]$ ) correspond to decreased errors for the fixed $[eqn]$ and h considered here. This implies that, for a given discretization, higher physical viscosity $[eqn]$ (leading to larger $[eqn]$ values) may result in improved numerical accuracy, while still maintaining $[eqn]$ to ensure stability. It is important to emphasize that $[eqn]$ is an inherent physical property of the system. The numerical parameters $[eqn]$ and h must be selected in coordination with the physical viscosity $[eqn]$ to guarantee that the relaxation time $[eqn]$ remains within a stable and accurate range, balancing computational cost and precision. When the physical $[eqn]$ is very small, causing $[eqn]$ to approach the lower stability limit of 0.5, adjustments to $[eqn]$ or h may be necessary to maintain both stability and accuracy, which often entails a trade-off with computational efficiency.

Example 3. Consider the following the diffusion equation in the region $[eqn]$ with initial and boundary conditions [48]:

$$[eqn]$$

The initial and boundary conditions are given below:

$$[eqn]$$

The corresponding exact solution for this instance is given as

$$[eqn]$$

In this example, we perform a basic sensitivity and stability analysis by simulating the model with two different sets of spatial and temporal resolutions: $[eqn]$ and $[eqn]$ for the simulation. To satisfy the Courant–Friedrichs–Lewy (CFL) stability condition in numerical simulations, we set $[eqn]$ for the case of $[eqn]$ , and $[eqn]$ for the case of $[eqn]$ . For both spatial step sizes, we set $[eqn]$ . This test is designed to evaluate the robustness of the proposed scheme under different discretization parameters. The simulation results indicate that, regardless of the selected spatial and temporal step sizes, the numerical solutions remain in excellent agreement with the exact solutions throughout the time evolution. The simulation results demonstrate that regardless of the values chosen for $[eqn]$ and $[eqn]$ , the numerical solutions maintain excellent agreement with the exact solutions during long-term evolution. Similar to previous examples, we employ smaller time steps and spatial steps to ensure high computational accuracy. Figure 9a,c display three-dimensional images of the exact solution at t = 0.5, showing a shape with a high center and low periphery, which presents certain computational challenges due to significant numerical variations. However, comparison between the three-dimensional numerical solution illustrations in Figure 9b,d and the exact solution illustrations reveals a high degree of concordance between numerical and exact solutions, independent of the selected values for $[eqn]$ and $[eqn]$ . These results provide preliminary evidence of the reliability and robustness of the proposed model. Further tests involving a wider range of parameters, including relaxation times, will be explored in future work.

Figure 10a,b present two-dimensional profiles at $[eqn]$ at different time intervals, providing detailed information about local function variations during the evolution process, further confirming the consistency between numerical and exact solutions. For error analysis, Table 5 display the MAXE, MAE, and GRE for this example. The numerical values in Table 5 indicate that errors remain minimal over extended time periods, validating the effectiveness and reliability of this model.

Example 4. Consider the following unsteady convection–diffusion equation with initial and boundary conditions in the domain $[eqn]$ [49]:

$$[eqn]$$

The initial condition is given by the following:

$$[eqn]$$

where

$$[eqn]$$

It can be verified that the exact solution is the following

$$[eqn]$$

The boundary conditions can be easily derived from the exact solution.

In the simulation, $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ . The three-dimensional surface plots of the two-dimensional numerical solutions at $[eqn]$ are illustrated in Figure 11, which provide a clearer visualization of spatial variations over the 2D domain. A close agreement between the numerical solution (left) and the exact solution (right) is observed, demonstrating the high accuracy and reliability of the proposed lattice Boltzmann model. Both solutions exhibit smooth, symmetric surfaces characterized by a pronounced peak near the center and a gradual decay towards the domain boundaries along both the x- and y-axes. The numerical method successfully reproduces the essential dynamics of the convection–diffusion system, capturing both the spatial smoothness and the outward decay in solution values. These results confirm the model’s capability to accurately resolve the combined effects of convection and diffusion in unsteady transport processes.

Furthermore, Figure 12 provides a two-dimensional comparison of the numerical and exact solutions along the line $[eqn]$ at multiple time instances. The numerical profiles align well with the exact solutions, confirming the model’s ability to accurately simulate the temporal evolution of the scalar field $[eqn]$ .

The contour plot in Figure 13 further illustrates the spatial distribution of the numerical solution at $[eqn]$ . The contours are smooth and well-behaved, indicating the stability of the numerical scheme. The results highlight the model’s capability to handle complex boundary conditions and accurately resolve the interplay between convection and diffusion.

The MAXE, MAE, and GRE values in Table 6 demonstrate the model’s high accuracy over time. Due to the system’s unsteadiness, the MAXE increases slightly from $[eqn]$ at $[eqn]$ to $[eqn]$ at $[eqn]$ , indicating that the numerical solution still remains close to the exact solution throughout the simulation. Similarly, the MAE values reflect the overall consistency of the numerical results. The GRE values remain stable around $[eqn]$ , further confirming the robustness and reliability of the proposed method for solving unsteady convection–diffusion equations. These results highlight the effectiveness of the lattice Boltzmann model in accurately capturing the dynamics of the system over time.

4. Conclusions

In this study, a lattice BGK model incorporating an amending function was developed to solve two-dimensional second-order nonlinear partial differential Equation (SONPDE). The Chapman–Enskog expansion rigorously demonstrated the model’s ability to recover the target macroscopic equations in the continuous limit, ensuring theoretical consistency between the mesoscopic framework and the macroscopic physics. A D2Q4 discrete velocity set was employed to construct the model, enabling efficient computation and clear derivation of the moment constraints for both equilibrium and correction terms. Extensive numerical simulations were conducted on multiple benchmark problems. The proposed model exhibited excellent agreement with exact solutions, confirming its accuracy, stability, and robustness across a variety of parameter regimes and initial conditions. Furthermore, the model demonstrated strong performance over long-term evolutions, showcasing its capability to capture essential convection–diffusion dynamics with high precision. The simplicity of the D2Q4 lattice and the modularity of the correction scheme make the approach computationally attractive for large-scale simulations and parallel implementations. This study establishes a solid foundation for further kinetic-based modeling of nonlinear PDEs in two dimensions.

Despite these promising results, the current model also faces some limitations. First, the use of the D2Q4 lattice restricts the method to relatively simple geometries and structured grids; its application to irregular or complex domains may require non-trivial modifications or hybrid approaches. Second, the extension to three-dimensional problems is not straightforward, as it demands careful design of the discrete velocity set and appropriate correction mechanisms to maintain accuracy and stability. Furthermore, the current formulation is tailored for scalar equations; extending it to coupled systems or vector-valued fields introduces additional modeling challenges that merit future investigation.

Future research will aim to (1) extend the current framework to nonlinear, coupled, higher-order or more general partial differential equations, such as reaction–diffusion systems, heat and mass transfer equations with variable coefficients, and fluid–structure interaction problems; (2) integrate multiple relaxation time (MRT) and entropic collision operators to enhance numerical stability and adaptability to high Reynolds number flows; (3) incorporate thermodynamic and hydrodynamic nonequilibrium effects to simulate strongly coupled or multiphysics processes such as reactive flows and multiphase transport; and (4) explore coupling with machine learning techniques for adaptive parameter tuning and model acceleration. These directions will not only expand the applicability of the proposed model but also facilitate its use in more realistic and complex scientific and engineering problems.