Papers

Vasudeva Murthy's relaxation approach [A.S. Vasudeva Murthy, J. Comput. Appl. Math., 203(2), pp. 437-443, 2007], originally proposed for the Jin-Xin relaxation model, provides an alternative formulation with invariant properties that is consistent and retains the semilinear structure incomparison to the standard one. In this work, Vasudeva Murthy's relaxation approach for various relaxation systems are proposed such as the shallow water equations, the Broadwell model, the Euler equations with heat transfer and two-dimensional Jin-Xin model. For proposed relaxation models, the associated integral invariants are rigorously established at the theoretical level. The main advantage of the integral invariant is that it provides a conserved quantity for the relaxation system by incorporating the coupled contributions of the solution variables in vector form. To validate the analytical results, numerical simulations are carried out for each model using three second-order numerical schemes: CS-EBT2, a semi-implicit second-order central finite-volume scheme for hyperbolic systems with relaxation source terms [S. Sahu, E. Macca, and R. Samala, J. Comput. Phys., 563, 115100, 2026]; UCS2, a finite-volume central relaxation-type scheme [S. F. Liotta, V. Romano, and G. Russo, SIAM J. Numer. Anal., 38(4), 1337-1356, 2000]; and IMEX-RK2, a second-order Implicit-Explicit Runge-Kutta scheme [Pareschi and Russo, J. Sci. Comput., 25, 129-155, 2005]. Numerical results, compared to exact or finely resolved reference solutions, confirm that the models preserve integral invariants, remain stable under CFL restrictions, and exhibit robust and accurate behavior across all benchmark systems tested.

The Poisson–Nernst–Planck (PNP) equations strongly couple ionic transport and electrostatic interactions through the Poisson equation, posing substantial numerical challenges under small permittivity and complex potential boundary conditions. Underlying these equations is a natural Wasserstein gradient-flow structure, in which the Poisson equation serves as a local realization of the nonlocal electrostatic interaction energy. Exploiting this structure, we formulate each time step as a constrained convex minimization problem where the ionic continuity equations and the Poisson equation are incorporated as linear constraints, allowing the concentrations, fluxes, and electrostatic potential to be updated simultaneously. The variational structure of the scheme intrinsically guarantees the dissipation of the original free energy, mass conservation, and nonnegativity of ionic concentrations under general electrostatic boundary conditions. Moreover, the framework is structurally modular: extending from classical to modified PNP models with steric interactions and concentration-gradient corrections requires only modifying the energy functional, while all structure-preserving properties are automatically retained. To efficiently solve the resulting large-scale constrained problems, we develop preconditioned and transformed primal–dual algorithms equipped with tailored fast dual solvers, namely DCT-based direct and Schur-complement iterative methods, that exploit the coupled block structure of the PDE constraints. Numerical experiments on classical and modified PNP systems demonstrate the accuracy and structure-preserving properties of the scheme, and show that the proposed algorithms converge reliably in strongly coupled small-permittivity regimes without significant growth in computational cost.

Magnetic composites used for the conversion of electrical energy often incorporate ferromagnetic inclusions insulated from each other to mitigate eddy current losses. Numerical models for these composites must be robust enough to address potential convergence issues arising from the presence of nonlinear magnetic inclusions and presence of significant confined eddy currents in the cell. This paper introduces an $\mathbf{h}$-conforming multiscale formulation for magnetic composites in a periodic framework. The proposed method uses the Heterogeneous Multiscale Method (HMM) with two mesoscale problems: a magnetoquasistatic problem for upscaling the homogenized magnetic flux density $\mathbf{B}_M$ and a magnetostatic problem for upscaling the macroscale incremental reluctivity $(\partial \mathbf{B}_M/\partial \mathbf{H}_M)$. Additionally, the method uses relaxed Newton–Raphson schemes at both macro and mesoscale levels to mitigate the well-known NR convergence issue linked to nonlinear BH constitutive laws in $\mathbf{h}$-conforming formulations. The accuracy and performance of the formulation are evaluated using 2D and 3D idealized periodic soft magnetic composites with linear and nonlinear BH curves. Furthermore, the paper demonstrates that the magnetoquasistatic mesoscale problem can be replaced by a magnet a

We introduce and analyze a barrier-regularized symmetric Nitsche method for the scalar Signorini problem. Applying a logarithmic barrier to the nonnegative slack variable in an augmented Lagrangian and then eliminating that variable yields a smooth positive-part operator and the perturbed complementarity relation on a primal-dual central path, with barrier parameter $\mu=\gamma s$. For every $s>0$, the discrete problem is continuously differentiable, uniquely solvable, and has a symmetric positive definite Newton matrix when the Nitsche parameter is sufficiently large. Under a mild barrier-feasibility condition, the continuous logarithmic energy has a unique minimizer $u_\mu$ for every $\mu>0$, its gap is positive almost everywhere, and $\|u-u_\mu\|_{H^1(\Omega)}\lesssim\mu^{1/2}$. Under the Sobolev regularity used for finite element approximation, this minimizer satisfies the \(L^2\) central-path boundary law. If the obstacle is locally the trace of an $H^2$-function, we also prove a uniform local $H^2$ bound in the interior of each planar contact face, on neighborhoods that may contain a free-boundary point of the limiting Signorini solution. The method is exactly consistent and quasi-optimal relative to $u_\mu$, with constants independent of $s$. If $u_\mu$ is uniformly bounded in $H^r(\Omega)$, $3/2<r\leq k+1$, then $\|u-u_h\|_{H^1(\Omega)}\lesssim h^{r-1}\|u_\mu\|_{H^r(\Omega)}+\mu^{1/2}$; hence the sufficient balance $s\lesssim h^{2r-1}$ preserves the available energy-norm rate. We also derive rates for discrete penetration and the complementarity residual. Numerical experiments with $P_1$ and $P_2$ elements on structured and unstructured meshes support the predicted rates and the local Newton theory, while showing that a rate-preserving smoothing may still resolve the contact set poorly.

This paper proposes a novel frequency-time hybrid multiple scattering (FTH-MS) integral equation solver for time-dependent wave equation problems in general multi-layered media, with remarkable scalability with respect to the number of layers $N$. In light of the finite speed of wave propagation, the new methodology provides an innovative multiple-scattering idea of re-modeling the original $N$-layered problem into a sequence of $N-1$ two-layered sub-problems, for which the main advantages lie in that (i) each sub-problem enjoys much simpler wave scattering properties compared with the complicated problem in a multi-layered medium, (ii) it enables to develop high-accuracy solver utilizing Fourier transform and frequency-domain boundary integral equation (BIE) method; and (iii) numerical evaluation of the sub-problems in each multiple scattering step can be parallelized. Both multiplicative- and additive-type strategies are developed and equivalence results, which indicate that the $M$-th order multiple scattering sums can provide equivalent representations of the solutions up to a certain time $T(M)$, are rigorously derived. Owing to the existed result of exponential convergence of the perfectly-matched-layer (PML) truncation for two-layered problem, all the sub-problems is numerically resolved by means of the FTH method based on the Fourier transform and the PML-BIE method whose numerical evaluation is addressed utilizing the Chebyshev-based rectangular-polar solver with high accuracy. Numerical examples are presented to validate the efficiency and accuracy of the proposed method.

This paper introduces a family of one-step implicit methods for solving scalar ordinary differential equations. At each step, the update uses a scaled angular mean of two vector-field evaluations, and the positive scale may vary from step to step. The leading terms of the local truncation error can be expressed in terms of the derivative of the signed curvature of the scaled solution graph. We prove that the proposed method reproduces the exact solution at the mesh points when the solution graph has constant signed curvature under a fixed positive scaling and each implicit update is unique. When this special geometric condition is not satisfied, we establish second-order consistency and convergence for positive scale sequences satisfying suitable uniform conditions. Furthermore, third- and fourth-order consistency and convergence can be achieved by choosing the scale to cancel the relevant curvature terms in the local truncation error. Using only the given problem data, we classify when these improvements are possible and determine the corresponding scale choices. An example shows that the proposed fourth-order method can yield smaller errors than the classical fourth-order Runge–Kutta method at the same step size.

Independently weighted semidefinite penalties arise in augmented-Lagrangian and constrained formulations. This paper characterizes when an exact additive subspace-correction preconditioner remains uniformly effective over all nonnegative penalty weights on a fixed finite-dimensional space. Robustness holds precisely when the correction spaces decompose every joint kernel generated by a nonempty subset of penalties. If one condition fails, a computable constant determines the exact first-order decay of the smallest preconditioned eigenvalue along the associated parameter ray, and the condition number grows linearly; none of the subset conditions can be discarded in general. Filtered decompositions provide computable sufficient bounds on parameter-ordering cones, while distributive kernel lattices permit a single common splitting. Exact-additive computations confirm the characterization and predicted rates. Separate Scott-Vogelius experiments produce stable multilevel iteration counts over the tested weights and mesh levels. The analysis does not establish mesh-uniformity.

We consider the collision-transport operator for vector-valued two-relaxation-time lattice Boltzmann schemes linearized about a uniform rest state. Assume that the equilibrium blocks are positive definite and that the link-even and link-odd relaxation parameters satisfy s_+ + s_- = 2 and 0 < s_- < 2. If the homogeneous transport is unitary in the equilibrium metric and reversible under velocity exchange, then the powers of the amplification operator are bounded uniformly with respect to the number and arrangement of lattice nodes. The admissible transports include periodic transport, vector halfway bounce-back, coordinate-aligned specular reflection, tangential orthogonal involutions, and compatible multi-channel scattering. The proof reduces the population equation to a two-step macroscopic recurrence generated by a contraction. An inclusion of the numerical range in an ellipse, combined with the Crouzeix-Palencia theorem, yields a dimension-independent estimate for the companion operator and hence the population bound. For a three-coefficient off-midpoint boundary interpolation, we give an exact rational D2N5 example whose finite-domain amplification matrix has a real eigenvalue larger than one, although the interpolation coefficients and bulk parameters are admissible. Thus coefficient convexity alone does not ensure stability for this boundary family.

Checkerboard edge-midpoint geometry provides an exact planar Varignon parallelogram for every spatial quadrilateral, allowing a local tangent frame and normal to be recovered directly from nodal positions even when the raw quadrilateral is warped. Building on this property, we develop a position-only thin-shell discretization with no independent director, rotation, or strain variables. The connected edge-midpoint surface is taken as the physical midsurface: the first fundamental form is evaluated on planar B faces, while a W-centered second fundamental form is constructed from variations of neighboring B-face normals, so membrane and bending share the same geometric carrier. A variational-kernel analysis shows that smooth second-order consistency does not eliminate lattice-scale blind modes. After quotienting out the raw checkerboard gauge, the B metric has one physical membrane blind direction and the symmetric W curvature has two curvature blind directions. We introduce a quotient-minimal membrane compatibility coordinate $X_M$ and an objective reference-relative curvature coordinate $X_W^{rel}$, placed consistently in the $O(t)$ membrane and $O(t^3)$ bending sectors. The formulation admits an explicit midpoint quotient, complete flat blind-mode classification, rigid-motion objectivity, reference-state consistency, and an $O(h^2)$ near-isometry approximation result for aligned generalized cylinders. Numerical tests show second-order curvature convergence, targeted removal of the membrane defect, and a sub-percent, refinement-decaying influence of $X_W^{rel}$. Linear and nonlinear shell benchmarks further demonstrate flat bending, curved-shell membrane-bending coupling, thickness sensitivity, large rotation, nonlinear pinching, and localized ovalization within a single position-only framework.

For $N\geq 3$ and a potential phase $\vartheta\in\mathbb{R}$, we study the non-self-adjoint almost Mathieu matrix obtained by multiplying the discrete Laplacian by a complex phase with angle $\varphi\in\mathbb{R}$, $A_N(\varphi,\vartheta)=e^{i\varphi}(S+S^{-1})/2+\operatorname{diag}(\cos(2\pi j/N+\vartheta))_{j\in\mathbb{Z}/N\mathbb{Z}}$, where $S e_j=e_{j+1}$ is the periodic shift on $\mathbb{C}^N$. We derive a Chambers formula and isolate the part $Q_{N,\varphi}$ of the characteristic polynomial that depends only on $N$ and $\varphi$, but not on $\vartheta$ or on a change of boundary conditions for the shift operator. We then show, for every $N$, that the zeros of $Q_{N,\varphi}$ lie on the two perpendicular lines $e^{i\varphi/2}\mathbb{R}\cup e^{i(\varphi/2+\pi/2)}\mathbb{R}$. For even $N$, the same property holds for the matrices $A_N(\varphi,\vartheta)$ with $\vartheta\in 2\pi\mathbb{Z}/N$, and we compute their limiting eigenvalue measure explicitly. For $\varphi\in[-\pi,\pi]$, the eigenvalue distribution approximates elliptic-integral densities with masses $1-|\varphi|/\pi$ and $|\varphi|/\pi$, and maximal radii $2|\cos(\varphi/2)|$ and $2|\sin(\varphi/2)|$, respectively. At $\varphi=\pi/2$, the central polynomial $Q_{N,\varphi}$ factors into positive quartic factors. This proves that the Scottish flag matrix, after Trefethen and Chapman, has its spectrum on the two diagonal lines of the saltire.

We develop a variational–Markov construction of detailed-balance master-equation discretizations for Fokker–Planck equations directly from their energy–dissipation laws. Rather than discretizing the differential operator, we represent mass transfer between neighboring grid points by two directional jump rates. On each edge, the discrete energy–dissipation law determines the net flux, while detailed balance fixes the rate ratio; a local reconstruction of edge quantities from neighboring grid-point values then uniquely determines both rates. A logarithmic-mean reconstruction recovers the classical Scharfetter–Gummel/Wang–Peskin–Elston rates, while alternative reconstructions yield other reversible schemes. The resulting semi-discrete systems conserve mass, preserve nonnegativity, satisfy detailed balance, and dissipate a discrete free energy. The same edgewise construction extends to state-dependent and degenerate mobilities, nonlocal interaction energies and higher-dimensional problems. The corresponding implicit and semi-implicit fully discrete schemes are linear, mass-conservative, positivity-preserving, and energy-stable, with a time-step restriction required only for nonlocal interaction energies whose kernel is not negative semidefinite. Numerical experiments confirm the expected spatial convergence rates of all detailed-balance schemes and verify their equilibrium accuracy, positivity preservation, and free-energy decay with discontinuous potentials, saturation, nonlocal interactions, and two-dimensional problems.

Sparse representations are often expected to make models smaller and also reduce inference cost. For Fourier Neural Operators (FNOs), these objectives are not equivalent or do not always align: removing parts of the learned operator can leave the underlying transforms and dense computations unchanged, while changing the grid on which the model is evaluated can introduce overhead of its own. We therefore distinguish sparsity in the representation, in the stored parameters, in the theoretical operation count, and in measured runtime, and present an empirical study of several routes toward sparse FNOs that tests each transition between them separately. Coarsening the execution grid reduces the theoretical cost without reducing measured latency, and adding a correction term recovers accuracy at the cost of making the model slower. Even an 83% parameter reduction remains slower than the dense baseline under ordinary execution. These results motivate a stricter definition of useful sparsity: the deployed operator must preserve solution accuracy and map its reduced support to a genuinely cheaper execution path.

We develop an immersed interface finite difference method for a one-dimensional nonlinear parabolic interface problem with jump condition \[ [u]_\alpha=\lambda u^+u^-. \] The method combines a Crank–Nicolson immersed interface discretization with an \(s\)-parameter reduction of the nonlinear interface condition, thereby reducing the nonlinear coupling to a scalar quadratic equation. We also discuss two different viewpoints for combining Newton iteration with immersed interface discretization, namely the IIM–Newton and Newton–IIM formulations. Numerical experiments are presented to illustrate the behavior and accuracy of the method.

We construct Fortin operators for discontinuous Petrov-Galerkin (DPG) discretizations of the advection equation with a piecewise constant divergence-free advection vector $\beta$, on simplicial meshes of any spatial dimension $N$ and for any polynomial degree $k \ge 1$ of the trial space. A minimal test space is built on each element from facet and interior bubbles and later augmented. The Fortin operator is shown to be bounded, uniformly over shape-regular mesh families, in the natural $\beta$-weighted broken test graph norm built on $L_q$ for every $1 < q < \infty$, where $q$ is the exponent conjugate to the trial exponent $p$. A non-characteristic facet condition is assumed when $k \ge 2$, while the lowest-order case requires no such condition and admits characteristic facets. As applications we prove that the practical fully discrete residual minimization method is quasioptimal in the DPG energy norm for every $1 < p < \infty$, with a quasioptimality constant governed solely by the Fortin operator, that its computable residual is a globally reliable and efficient a posteriori error estimator, and that augmenting the test space and changing the test norm improves the results in the lowest-order case.

Finite difference methods for interface problems often exhibit large local truncation errors near the interface, particularly when discontinuities in coefficients or solution derivatives are present. Nevertheless, many such schemes achieve second-order accuracy, a phenomenon not fully explained by standard pointwise consistency arguments. In this work, we provide a precise explanation of this behavior through the structure of the discrete Green's function associated with the underlying conservative difference operator. An explicit representation of the Green's function reveals a two-plateau structure in its weighted increments. By expressing the numerical error in terms of this Green kernel, we show that the dominant interface truncation errors, although individually of order \(O(1)\), possess a cancellation structure that reduces their effective contribution to the global error. The exact Green's-function analysis leads to a class of conservative interface flux–balance schemes for which the cancellation mechanism yields second-order accuracy in the maximum norm. The weighted harmonic discretization is included as a particular member, and its cancellation property is established directly. For Cartesian grids with a flat interface parallel to a coordinate direction, the same conservative cancellation mechanism persists along grid lines crossing the interface. Numerical experiments in one and two dimensions directly illustrate the predicted cancellation behavior and the resulting second-order accuracy. These results demonstrate that second-order accuracy in interface problems can arise from the interaction between localized truncation errors and the global structure of the discrete operator, rather than from pointwise consistency alone.

In this paper, we propose an efficient preconditioner for solving indefinite complex symmetric linear systems within a block preconditioning framework. We analyze the convergence of the corresponding iterative method and investigate several spectral properties of the preconditioned matrix, including eigenvalue distributions and eigenvector structures. The new preconditioner is used to accelerate the convergence of the flexible version of GMRES. Numerical experiments are presented to illustrate the effectiveness of the proposed preconditioner, and comparisons with existing block preconditioners demonstrate its superior performance.

This paper focuses on computing spectral solutions for high-order eigenvalue problems using an efficient discretization method based on Chebfun spectral discretization algorithms and domain truncation. We solve several numerical eigenvalue problems, demonstrating both the accuracy and computational efficiency of the proposed approach.

We construct high-order Runge–Kutta extensions of Barrett–Garcke–Nürnberg (BGN) parametric finite element methods for curve-shortening flow and curve diffusion of planar curves and for mean curvature flow and surface diffusion of closed genus-$0$ surfaces. On each time slab $I_m=(t_m,t_{m+1}]$, the continuous equations are posed on the left-endpoint surface $\Gamma^m=\Gamma^{t_m}$ through the map $\X^{m,t}:\Gamma^m\to\Gamma^t$, whose target is the evolving surface at time $t$. For curves, this formulation follows by harmonic pullback from $\Gamma^0$ to $\Gamma^m$. For surfaces, an orientation-preserving harmonic diffeomorphism is posed separately on each time slab, and conformality yields the weight $\frac12|\nabla_{\Gamma^m}\X^{m,t}|^2$. Evaluating the time-slab equations at the Runge–Kutta internal times and applying mass-lumped parametric finite elements yields systems on the common domain $\Gamma^m$, with $\Gamma^{t_m+c_i\tau_m}$ as the intermediate target geometry. This internal-time discretization constructs high-order BGN-structured extensions for the four flows. For an algebraically stable tableau with nonnegative weights, every exact stage solution with nondegenerate intermediate configurations satisfies monotone decay of the discrete curve length or surface area for every positive time step for which that solution exists. Radau IIA experiments exhibit energy decay for all four flows and BGN-type mesh redistribution in the curve tests. The Hausdorff self-convergence results exhibit high-order behavior consistent with the corresponding design orders.

We extend the Wigner-Eckart factorization of the spectral Boltzmann collision operator to polyatomic gases with continuous internal energy. Because internal energies are invariant under spatial rotations, the SO(3) reduction survives the Borgnakke-Larsen energy exchange, and the twelve-dimensional collision integral collapses onto a nine-dimensional kinematic core. The core splits into a sparse geometric tensor, evaluated exactly, and a dense physical tensor, integrated by singularity-resolving Gauss rules with an auxiliary Laplace representation of the fractional energy couplings. The quadrature attains near machine precision at the fractional exponents of real gases. The collision invariants are embedded exactly, preserving the translational-internal energy exchange. The factorization compresses the operator by three to nearly four orders of magnitude and accelerates its evaluation 40-fold over dense formulations. The method is validated against the exact monatomic limit, Landau-Teller relaxation, and an analytic frozen-channel Prandtl number, and it matches a published calibration of the same kernel for N2, CO, and H2.

This study introduces and analyzes a three-step sequential decoupling algorithm designed to address nonlinear, fully coupled quasi-static thermo-poroelasticity systems incorporating convective transport. The finite element method is employed for spatial discretization and the backward Euler method for temporal discretization. The proposed sequential method has a higher computing efficiency than the fully implicit nonlinear numerical scheme, since it does not require any internal iterations. The well-posedness of the numerical solution is discussed by introducing a cut-off operator and the stability analysis of the algorithm is performed. Rigorous analysis yields optimal convergence order estimates for both spatial and temporal discretizations. In order to confirm the theoretical results and the effectiveness of the suggested approach, numerical experiments are finally carried out.

Special relativistic hydrodynamics (SRHD) equations arise in the modeling of high-speed fluid flows encountered in astrophysical phenomena such as jets, supernova explosions, and gamma-ray bursts. Owing to their highly nonlinear hyperbolic nature, solutions often develop strong discontinuities, making the design of stable and accurate numerical schemes challenging. Although Weighted Essentially Non-Oscillatory (WENO) schemes are widely used for such problems, component-wise WENO reconstruction may produce spurious oscillations near discontinuities. On the other hand, characteristic-wise WENO reconstruction provides accurate non-oscillatory solutions for systems of conservation laws, but it involves the computation of eigenvectors in each cell, which leads to high computational cost. In this work, we intend to develop hybrid schemes which maintain the non-oscillatory feature of characteristic-wise WENO while being less costly. We propose three hybrid schemes, namely the H1-WENO, H2-WENO, and H3-WENO schemes, based on a new troubled-cell indicator constructed from the smoothness indicators of the WENO scheme. The proposed troubled-cell indicator effectively distinguishes smooth and discontinuous regions, allowing the hybrid schemes to employ inexpensive reconstructions in smooth regions and the characteristic-wise WENO reconstruction only near discontinuities. Numerical experiments demonstrate that the proposed schemes retain the accuracy and robustness of characteristic-wise WENO methods while significantly reducing the computational cost. In particular, the H1-WENO scheme achieves an approximately 30–40% improvement in computational efficiency compared to the standard WENO scheme in 2D test cases.

We develop a high-order physical-constraint-preserving (PCP) alternative finite difference weighted essentially non-oscillatory (AFD-WENO) scheme for the special relativistic hydrodynamics equations with general equations of state. The proposed scheme comprises two key limiters: a state limiter, which acts after the WENO state interpolation step, and a flux limiter, which acts on the final high-order fluxes. The state limiter ensures that the interpolated states are physically admissible, while the flux limiter ensures that the numerical fluxes are physically admissible. The resulting scheme is rigorously proved to satisfy the physical constraints. Incorporating multiple WENO interpolation techniques, including an improved adaptive-order formulation (WENO-AOI), the method is validated through extensive one- and two-dimensional numerical benchmarks with various equations of state. The numerical results demonstrate high-order accuracy, sharp resolution of discontinuities, and robust stability in extreme relativistic regimes.

Tensor representations have emerged as a fundamental paradigm for modeling multidimensional data by preserving intrinsic correlations across multiple modes. This paper proposes a novel theoretical framework, termed Tensor Orthogonal Subspace Split (TOSS), which explicitly splits a tensor, along a prescribed mode, into two orthogonal components: a dominant component lying in a prescribed subspace and a residual component lying in the corresponding orthogonal complement. We first present the general formulation of TOSS and systematically investigate its fundamental properties. As an important and practically meaningful special case, we further introduce the rank-one TOSS, which imposes a separable rank-one structure along the splitting mode and admits a clear geometric interpretation. This formulation naturally captures dominant consistent patterns while effectively isolating orthogonal residual component. The proposed framework establishes a unified theoretical foundation for tensor-domain orthogonal split and opens new avenues for structured tensor modeling across diverse applications. Building upon the developed TOSS theory, hyperspectral image restoration and color video background modeling are considered as two representative tasks, for which corresponding optimization models are formulated. Efficient algorithms are developed to solve the resulting problems. Extensive experimental results validate the effectiveness and superiority of the proposed approaches.

This paper investigates a stochastic multiscale subdiffusion model driven by fractional Gaussian noise, where the multiscale Abel kernel with variable exponent $\alpha(t)\in(0,1)$ is used to capture multiscale and crossover behavior in anomalous diffusion. The main difficulties of this model lie in the complexity of the multiscale Abel kernel (e.g. non-monotonicity and non-coercivity) and the low regularity caused by the noise. Concerning these issues, we prove the well-posedness and regularity of the mild solutions by means of solution operator approach and a perturbation technique for multiscale Abel kernel. Then both the semidiscrete-in-time and fully-discrete numerical schemes are proposed and analyzed under the low-regularity numerical analysis framework, with proved temporal and spatial convergence rates. Numerical experiments are presented to substantiate the theoretical results.

Physics-informed neural networks (PINNs) have emerged as a meshless alternative to conventional numerical methods for solving partial differential equations (PDEs). However, their limited ability to capture sharp gradients can lead to substantial errors when resolving boundary and interior layers. Here, we introduce an energy-norm-enhanced PINN (ENPINN) that incorporates gradient information and variational structure into the loss function to improve the resolution of layer-dominated solutions. We first examine two related formulations: weak-loss PINNs (WLPINNs), which incorporate test functions into the conventional PINN residual, and gradient-enhanced PINNs (gPINNs), which augment the loss with spatial derivatives of the PDE residual. By analyzing these formulations, we identify their limitations in resolving steep solution gradients and motivate the systematic construction of ENPINN. We establish theoretically how the energy-norm error depends on the ENPINN loss and show that a suitably modified residual-derivative term is essential for accurately capturing boundary layers. We further establish the existence of neural-network approximations with arbitrarily small energy error and derive corresponding derivative bounds, providing a theoretical foundation for the proposed framework. The performance of ENPINN is assessed through systematic comparisons with existing PINN variants for convection-diffusion-reaction problems exhibiting steep gradients. Numerical experiments include a combustion model, a coupled multi-scale system, a two-dimensional Burgers equation with an interior layer, and a three-dimensional time-dependent problem.

We analyse a non-residual variational multiscale finite element formulation with dynamic subscales for incompressible generalised Newtonian Navier–Stokes flows. The apparent viscosity is assumed to be bounded and Lipschitz continuous with respect to the shear rate, covering several regularised rheological laws. The method separates, through orthogonal projections, the unresolved pressure-divergence, convective, and pressure-gradient contributions. For the linearised semi-discrete problem, we prove well-posedness and an unconditional stability estimate in an anisotropic VMS norm. The estimate controls viscous dissipation, subscale energies, discrete divergence, and a coupled acceleration-convection-pressure balance. For the non-linear formulation, a fixed-point argument yields existence and optimal-order a priori bounds under suitable regularity and smallness assumptions. The analysis also provides weak-in-time control of the discrete pressure-convection balance. A key point is that the dynamic pressure-related subscale provides the time-derivative contribution needed to control the projected pressure-gradient error, a mechanism unavailable in the corresponding quasi-static setting.

The uncontrolled equations of motion in the Nechvile frame for the restricted three-body problem have been well-known since at least the 1960s. It would seem that adding an external force to these equations is quite trivial: simply add an external force per mass term to the acceleration equations. Here we show that the last statement is not true. In fact, we show that the additive generic external force must be multiplied by the inverse of $(1 + e \cos\theta)^3$ where $e$ is the relative eccentricity of the primaries and $\theta$ is the true anomaly of the rotating frame located at the barycenter. Furthermore, when this result is combined with the mass flow rate equation, it generates several surprising results due to the mismatch between the resulting quadratic term and the cubic term in the equations of motion. This leads to a corresponding modification of the rocket equation itself. A Birkhoff-theoretic solution to an illustrative cislunar space mission problem shows propellent savings of 80% with the use of the correct cost functional. The popular quadratic cost utilizes more than $2X$ the minimum propellant consumption.

We develop an entropy-stable and physical-constraint-preserving discontinuous Galerkin spectral element method for symmetry-reduced general-relativistic hydrodynamics on prescribed stationary spacetimes. Using a local orthonormal transformation, the fluid variables are expressed in a form for which the relativistic hydrodynamic algebra and the admissible set are independent of the spatial metric, while the spacetime geometry enters through stationary coefficients. This separation allows entropy-conservative special-relativistic fluxes to be combined with a compatible discretization of the geometric source terms. On affine tensor-product meshes, the resulting DGSEM is conservative and satisfies a semidiscrete entropy inequality, while the transformed variables provide a convex framework for physical-constraint preservation. For practical stabilization, we use a geometry-only causal speed that is sufficient for both classical local Lax–Friedrichs entropy dissipation and the physical-constraint-preserving Lax–Friedrichs splitting. The fully discrete method combines this stabilization with SSP Runge–Kutta time stepping, oscillation elimination, and conservative local-orthonormal-state scaling. Numerical experiments cover smooth and strongly shocked special-relativistic flows, an axisymmetric jet, stationary Michel accretion, Schwarzschild Bondi–Hoyle flow, and four Kerr accretion cases. The results demonstrate the designed high-order accuracy in smooth regimes and robust performance for demanding relativistic flows on curved stationary backgrounds.

We develop a high-order tensor neural network (TNN) method for locally periodic elliptic multiscale problems of the form $-\nabla\cdot(A(x,x/\varepsilon)\nabla u_\varepsilon)=f$. Because the coefficient depends on both the slow variable $x$ and the fast periodic variable $y=x/\varepsilon$, the high-order cell problems and macroscopic corrector equations are more involved than in the classical case $A=A(y)$, and the correctors depend parametrically on $x$. We derive a computable high-order two-scale expansion and prove an $H^1$ convergence estimate for the partial expansion in boundary-layer-free settings, including periodic domains and ideal boundary-matching configurations. The proof uses the recursive compatibility structure of the corrector hierarchy and a zero-mean oscillation estimate in $H^{-1}$. We then construct a TNN framework for the high-dimensional corrector problems. Its tensor-product structure permits deterministic one-dimensional quadrature for the cell problems, homogenized coefficients, macroscopic source terms, and loss functions, avoiding Monte Carlo integration error. The numerical realization assumes that the coefficient entries and assembled data admit finite or controlled tensor-product representations; this computational assumption is separate from the general matrix-valued coefficient class used in the analysis. Experiments with scalar locally periodic coefficients show accurate high-order correctors. The $H^1$ semi-norm errors are consistent with the proved estimate, while point-normalized $L^2$ errors display the nominal high-order behavior predicted by the formal expansion.

For a stochastic autoparametric block-and-pendulum system, the long-time dynamics exhibit two fundamental features: the almost-sure stability of the single mode solution, characterized by its Lyapunov exponent, and the global asymptotic dynamics when this single mode solution loses stability. This naturally raises the question of whether these dynamical features are preserved under discretization, since such preservation is essential for the resulting discrete system to faithfully capture the qualitative behavior of the continuous system. To address this question, we first establish the existence of a random attractor for the continuous system subject to multiplicative stochastic excitation, providing a rigorous characterization of the global asymptotic dynamics. We then propose a numerical discretization that induces a discrete random dynamical system and prove the convergence of its random attractor to the continuous one as the step size tends to zero. In addition, we show that the numerical Lyapunov exponent of the single mode solution has the same sign as its continuous counterpart for sufficiently small step sizes, thus preserving the corresponding almost-sure stability or instability classification. These results demonstrate that the proposed discretization captures both the global asymptotic dynamics and the stability characteristics of the underlying stochastic autoparametric system.