In this paper, we propose a high-order weighted compact nonlinear scheme (WCNS) that is well-balanced and positivity-preserving for the two-dimensional Ripa system. We use a source term splitting technique and the surface reconstruction method to preserve the still-water steady states. Meanwhile, we introduce a “draining” time-step technique to ensure the positivity of the water depth and the temperature. This approach can effectively handle the wet-dry front problem, which is one of the main contributions of the present scheme. To preserve the isobaric steady states, we apply a steady-state-preserving parameter to modify the original Lax-Friedrichs flux. The moving-water steady states are preserved by reconstructing equilibrium variables rather than conservative variables, along with a source-term discretization that matches the flux discretization. Several numerical experiments are conducted on classical one-dimensional and two-dimensional problems of the Ripa system. Numerical results demonstrate the well-balanced and positivity-preserving properties, high-order accuracy, and excellent shock-capturing performance of our proposed scheme.
Fractional-order models provide a powerful framework for capturing anomalous transport, memory, and nonlocal interactions in biological microvascular networks. This study presents a hybrid 1D-3D model of cerebral drug transport under mild hyperthermia, coupling fractional intravascular dynamics with non-Fourier tissue bioheat. The vascular network is represented as a directed 1D graph, where advection, diffusion, and reactive exchange obey Caputo-Fabrizio time-fractional and Riesz space-fractional laws, linked through Robin-type mass and heat transfer to surrounding 3D tissue governed by a Cattaneo-Vernotte bioheat equation. Thermal feedback modifies permeability and perfusion, yielding a two-way thermo-chemical coupling. The explicit fractional scheme employs exponential-memory updates for Caputo-Fabrizio derivatives and symmetric Grönwald-Letnikov sums for nonlocal spatial fluxes while preserving global conservation. Simulations on an arteriole-capillary-venule network show that increasing fractional order α sharpens pulse dispersion and delays washout; higher temperature enhances vascular-tissue exchange and accelerates equilibration; and dual-phase-lag relaxation mitigates nonphysical heat spikes near vessel walls. The framework unifies anomalous vascular transport, temperature-sensitive physiology, and network geometry into a scalable model suitable for hyperthermia-optimized brain drug delivery and parameter calibration from imaging data.
The discovery of partial differential equations (PDEs) is crucial for applied science and engineering. However, data-driven discovery of PDEs is challenging due to the sensitivity of the equations to noise and the complexities involved in model selection. We propose a Bayesian sparse learning algorithm for discovering PDEs with variable coefficients, particularly when these coefficients are spatially or temporally dependent. Specifically, we apply threshold Bayesian group Lasso regression with a spike-and-slab prior (tBGL-SS) and leverage a Gibbs sampler for Bayesian posterior estimation of PDE coefficients. This approach not only enhances the robustness of point estimation with valid uncertainty quantification but also relaxes the computational burden from Bayesian inference by incorporating coefficient thresholds as an approximate MCMC method. Moreover, from the quantified uncertainties, we propose a Bayesian total error bar criteria for model selection, which outperforms classic metrics including the root mean square and the Akaike information criterion. The capability of this method is illustrated by the discovery of several classical benchmark PDEs with spatially or temporally varying coefficients from solution data obtained from the reference simulations. In these experiments, we show that the tBGL-SS method is more robust than the baseline methods under noisy environments and provides better model selection criteria along the regularization path.
Wall shear stress (WSS) is crucial in the development and rupture of abdominal aortic aneurysms (AAA). While computational fluid dynamics (CFD) enable the acquisition of hemodynamics during blood flow modeling, they are computationally expensive. To address this, we propose a deep learning (DL) framework integrating CFD, DL networks, a lightweight point cloud dataset, and a segment threshold for predicting dynamic WSS in AAA. Our workflow generates a WSS dataset from CFD simulations of 1000 AAA models, mapping temporal-spatial information to three dimensional WSS distribution using PointNet and PointNet++ networks. The DL framework accurately predicts WSS distribution with a mean relative error (MRE) below 12% and a normalized mean absolute error (NMAE) under 9%, closely matching traditional CFD results. Introducing a size ratio (SR) threshold, related to aneurysm rupture, for segmenting the test dataset improves predictive performance by 14% for SR greater than 2.7. This segmentation has significant clinical implications. Additionally, our use of lightweight point cloud data from the arterial wall enhances training efficiency five-fold compared to previous studies. This framework efficiently maps geometric information to dynamic WSS distribution, reducing computational cost and improving clinical applicability.
We present a numerical procedure for computing guaranteed two-sided bounds on the effective coefficients of elliptic partial differential operators in a three-dimensional setting. The upper bounds are obtained via a standard variational formulation, discretized using the finite element method. To derive the lower bounds, we formulate the corresponding dual variational problem and construct suitable approximation spaces within the finite element framework. To reduce the computational overhead associated with evaluating the lower bounds, we propose an efficient estimation technique based on a single FFT-based projection. Theoretical justification of the proposed procedure is provided, and its performance is demonstrated through illustrative numerical examples.
Turbulence modelling for incompressible flows remains challenging when strong anisotropy and intercomponent energy transfer are essential, as in thermal-hydraulics or atmospheric boundary layers. To address the limitations of eddy-viscosity closures, we consider the second-order turbulence-moment equations that transport the full Reynolds-stress tensor, with intercomponent energy redistribution modelled by the classical Rotta return-to-isotropy closure. A key difficulty in numerical simulations of such equations lies in maintaining the realisability of the Reynolds-stress tensor–symmetry and positive semi-definiteness–at the discrete level to avoid loss of hyperbolicity and numerical breakdown. We introduce a global two-step algorithm for multidimensional incompressible Reynolds-averaged Navier–Stokes equations: (i) an explicit Godunov-type convective predictor for stable transport of Reynolds stresses; and (ii) an implicit correction step applied to the mean velocity to enforce incompressibility and account for pressure, while the Reynolds stresses are updated using a dedicated diffusion-source operator. A key innovation is a semi-implicit, realisability-preserving strategy embedded in both the convection step and the integration of the Rotta source term, guaranteeing positive semi-definiteness in every cell and substep. Validation includes tests of the Rotta model using ordinary differential equations, where the Reynolds stress trajectories remain within the Lumley triangle, and full simulations of a plane mixing layer. Results capture turbulence onset, self-similar growth, and Reynolds stress evolution, with good agreement to experimental data even on coarse meshes.
This paper studies optimal control problems governed by the Boussinesq equations with box and L2-norm control constraints. The discretize-then-optimize strategy is utilized to address these problems. First, we derive the first-order optimality conditions for both types of control constraints. Subsequently, we discretize the problem using the L2 projection technique, the Marchuk-Yanenko method, and the P1−P1 iso P2 finite element method. Based on the optimality conditions of the discretized optimization problem, we develop algorithms for the box-constrained and L2-norm constrained cases by combining the L-BFGS method with an active set strategy, step size truncation, and projection techniques, as appropriate. The proposed algorithms are efficient and robust for solving the optimal control problems constrained by the Boussinesq equations. Numerical experiments are conducted to illustrate the effectiveness of these algorithms.
This study introduces an innovative framework that integrates radial basis functions (RBFs) into physics-informed neural networks (PINNs) for solving three-dimensional forward and inverse elastostatic problems in functionally graded materials. The proposed framework employs a fully-connected neural network that takes the angles of source points as inputs and predicts their locations along with shape parameters to approximate solutions across the computational domain. By leveraging the established RBF formulation and automatic differentiation, a loss function based on the governing equations and boundary conditions is constructed. The trainable parameters are optimized using a gradient-based constrained nonlinear optimization method, implemented by the “fmincon” function with the SQP algorithm in MATLAB to refine the source point distribution and shape parameters. The numerical results demonstrate that the proposed approach can effectively and accurately solve elastostatic and inverse Cauchy problems in three-dimensional complex geometries. Compared with the tested PINN and RBF collocation methods, the developed framework generally achieves improved accuracy and stable numerical performance.
This paper investigates the well-posedness and error analysis of a semilinear stochastic ψ˜-Caputo fractional diffusion-wave equation driven by fractionally integrated multiplicative noise. The existence, uniqueness, and regularity of mild solutions are established by combining the Banach fixed point theorem with the smoothing properties of the associated fractional solution operators. A novel fully discrete numerical scheme is developed using a spectral Galerkin method for spatial discretization together with a Mittag-Leffler Euler method for time integration. Rigorous error estimates are derived, yielding strong convergence results for the proposed scheme. Numerical experiments are presented to support the theoretical analysis and to demonstrate the effectiveness of the proposed approach.
We develop a fictitious-domain spectral method for reaction–diffusion systems defined on simply connected planar domains with curved boundaries and homogeneous Neumann conditions. By embedding the physical domain into an auxiliary disk, a polar change of variables and a Fourier expansion in the angular direction reduce each time step to a family of one-dimensional radial modified-Helmholtz problems. We split each radial problem into homogeneous and particular parts. The particular part is discretized with mode-dependent Legendre-Galerkin bases that satisfy the pole conditions at the origin, and we recover the unknown artificial boundary data from the physical Neumann condition via a boundary least-squares formulation stabilized by truncated singular value decomposition. For reaction-diffusion systems with cross-diffusion, we reduce the diffusion matrix to upper-triangular form via a real Schur decomposition. This yields sequentially decoupled scalar problems, allowing the same spatial solver to be reused without altering the radial discretization. Numerical experiments on scalar and coupled test problems confirm theoretically predicted temporal accuracy of the generalized BDF schemes, and simulations of the Schnakenberg and CIMA models on disks and flower-shaped domains demonstrate that the method captures geometry-induced mode selection and cross-diffusion-driven pattern transitions.
In this paper, we develop and analyse robust weak Galerkin methods for the steady Boussinesq model with damping. The methods respectively apply piecewise polynomials of degrees m (m ≥ 1), m−1 and m to approximate velocity, pressure and temperature variables in the interior of elements, and piecewise polynomials of degrees k(k=m−1,m), m and k respectively to their numerical traces on interfaces elements. The methods are demonstrated to yield globally divergence-free approximation for the velocity. Stability conditions, existence and uniqueness results of proposed discrete schemes are established, and optimal a priori error estimates in the energy norm and L2 norm are rigorously derived. A convergent linearized iterative method is proposed to solve the resulted nonlinear algebraic system. Numerical experiments are given to verify the theoretical analysis of presented algorithms.
The transmission model of the motor neural system is a nonlinear dynamical system wherein the action potential state is governed by the complex interplay of its intrinsic properties, axonal inputs, and prior states. This dynamical system is described by the FitzHugh-Nagumo (FHN) model, which presents severe computational challenges for direct solution owing to its strong nonlinearity and the inclusion of small parameters (e.g., ϵ) that critically influence neural impulse propagation. To balance computational accuracy, efficiency, and physical fidelity, we first construct a fully discrete finite element (FE) scheme incorporating a time-averaged linearization of the nonlinear term, effectively mitigating numerical stiffness. We then integrate this scheme with proper orthogonal decomposition (POD) to derive the ROFE method. Crucially, by leveraging the linearly-implicit structure of the parent FE scheme, the ROFE method circumvents costly hyper-reduction techniques while achieving substantial dimensionality reduction. Rigorous theoretical analysis establishes optimal error estimates for both the FE and ROFE solutions, explicitly quantifying the coupling effects of mesh size, time step, and POD basis number. Numerical experiments validate the theoretical findings and demonstrate that the ROFE scheme achieves over 80% computational time reduction (up to 89% for 2D) with only 3-5 POD modes capturing over 98% of system energy. Most importantly, both schemes strictly preserve the FHN model’s core biophysical property: normal neural impulse propagation occurs exclusively for the physiological threshold α ∈ (0, 1]. By guaranteeing numerical stability, optimal convergence, and physical consistency, this work provides a reliable and efficient framework for large-scale simulations of nonlinear reaction-diffusion systems, with direct applications in motor neural dynamics and excitable media.
In this paper, some high-order finite difference mapped unequal-sized WENO schemes with adaptive linear weights (which are termed as the ALW MUS-WENO schemes) are proposed for solving compressible flows in one, two, and three dimensions on structured meshes. The information based on three unequal-sized spatial stencils is employed to reconstruct one fourth degree polynomial, one sixth degree polynomial, and one eighth degree polynomial together with two linear polynomials for achieving high-order spatial approximations. By automatically adjusting three linear weights with only one simple condition, the fifth-order, seventh-order, and ninth-order ALW MUS-WENO schemes can attain uniformly high-order accuracies in smooth regions and maintain essentially non-oscillation properties near strong discontinuities. Then, a mapping function is applied to extremely decrease the finite difference between the linear and nonlinear weights. In comparison to the originally fifth-order finite difference US-WENO scheme, the proposed increasingly high-order ALW MUS-WENO schemes have the following advantages. First, they allow the deployment of a tiny ε=10−30 to obtain high-order accuracies even near critical points in smooth regions and suppress spurious oscillations near strong discontinuities, while the original US-WENO scheme fails to do so. Second, these new ALW MUS-WENO schemes track the limitations of the original US-WENO schemes which cannot provide the sharp shock transitions via the order increasing. Third, these high-order ALW MUS-WENO schemes are implemented to solve Euler and Navier-Stokes equations in one, two, and three dimensions, respectively, while the originally fifth-order US-WENO scheme is only confined for Euler equations in one and two dimensions. Numerical tests including the inviscid flows and viscous flows are conducted to demonstrate the attractive capability of these new ALW MUS-WENO schemes in multi-dimensions.
This work presents a nonlinear fourth-order partial differential equation (PDE) model for image denoising under strong multiplicative noise. Unlike existing second-order diffusion models, which often produce staircase artifacts and structural degradation, the proposed Laplacian-based fourth-order formulation employs gray-level and gradient-based edge information into the diffusion coefficient within a nonlinear biharmonic framework, thereby enabling effective multiplicative noise removal while preserving edges and fine structures. The existence and uniqueness of a weak solution are established using the Faedo-Galerkin framework together with compactness arguments and the Schauder-Tychonoff fixed-point theorem. Stability for both linear and nonlinear regimes is analyzed through an energy-based Lyapunov approach, demonstrating asymptotic convergence toward a nonlinear biharmonic steady state. Numerical experiments on natural and SAR images show that the proposed model consistently outperforms two recent second-order diffusion models in terms of PSNR, MSSIM, and structural preservation while maintaining strong denoising capability.
In this study, we present a comprehensive spectral analysis of the convection-dispersion equation, a class of equations that includes the Korteweg–de Vries (KdV) and modified Korteweg–de Vries (mKdV) equations as special cases, to investigate the behavior of high-order numerical schemes across a wide range of nondimensional parameters. The motivation for this analysis stems from the equation’s importance in modeling wave propagation and transport phenomena, where accurate resolution of dispersive effects is critical, and traditional numerical schemes often suffer from spurious artifacts. We analyze one sixth-order and two eighth-order compact finite difference spatial discretization schemes, encompassing both Cell-Node Compact Scheme (CNCS) and Central Compact Scheme (CCS) formulations, combined with a third-order strong stability-preserving Runge-Kutta (SSPRK3) time integrator. The analysis is performed in terms of key nondimensional parameters such as the wavenumber, Courant number (Nc), and dispersion number (Dϵ) over the full spectral plane for both one- and two-dimensional cases. Key numerical indicators, including the amplification factor, normalized phase speed, and normalized group velocity, are evaluated to characterize stability, dispersion error, errors in energy transport, and directional anisotropy. Critical dispersion thresholds and Courant numbers are identified, beyond which numerical instability and nonphysical phenomena such as spurious q-waves and reversed phase or energy transport arise. Theoretical predictions are verified through numerical experiments involving linear and nonlinear one- and two-dimensional test problems, including cases with exact solutions and established benchmark results. This comprehensive analysis uncovers subtle numerical errors and offers practical guidance for selecting reliable discretization parameters, ensuring accurate and stable simulations of convection-dispersion systems.
Nonlinear fractional-order unsteady weakly singular integral equations are encountered in many scientific fields. The combination of weakly singular kernels, nonlinear terms, and fractional-order derivatives presents significant analytical challenges. To address these, we have developed and applied two computational schemes to approximate solutions for such equations. For steady problems, we introduce a linearized fully spectral method, that represents u(x) using shifted Gegenbauer polynomials (SGPs) in the vector Λ(x), and use Picard's iterative method for the nonlinear terms. For unsteady problems, we extend a semi-discrete scheme that estimates temporal derivatives with a forward difference formula and spatial variables with the SGPs vector Λ(x). We also establish new operational matrices to estimate singular integral terms, such as:∫0x∫0y∫0zK(x,p)uγ(p,t)/(xρ1−pρ1)α1(yρ2−qρ2)α2(zρ3−rρ3)α3drdqdp;with ρ's > 1, 0 < α's < 1.These methods convert the original nonlinear problem into a system of linear algebraic equations that are straightforward to solve. We implement both schemes in Maple 2015 and validate our results against existing literature. Various test problems are used to demonstrate the accuracy, stability, and reliability of the proposed schemes. Simulations across a range of nonlinear parameters γ, α, β, ρ, M, and N confirm that our methods are accurate, stable, and effective for these challenging problems.
