Search

The approximation of solutions of partial differential equations (PDEs) with numerical algorithms is a central topic in applied mathematics. For many decades, various types of methods for this purpose have been developed and extensively studied. One class of methods which has received a lot of attention in recent years are machine learning-based methods, which typically involve the training of artificial neural networks (ANNs) by means of stochastic gradient descent type optimization methods. While approximation methods for PDEs using ANNs have first been proposed in the 1990s they have only gained wide popularity in the last decade with the rise of deep learning. This article aims to provide an introduction to some of these methods and the mathematical theory on which they are based. We discuss methods such as physics-informed neural networks (PINNs) and deep BSDE methods and consider several operator learning approaches.

We introduce a new level-set shape optimization approach based on polytopic (i.e., polygonal in two and polyhedral in three spatial dimensions) discontinuous Galerkin methods. The approach benefits from the geometric mesh flexibility of polytopic discontinuous Galerkin methods to resolve the zero-level set accurately and efficiently. Additionally, we employ suitable Runge-Kutta discontinuous Galerkin methods to update the level-set function on a fine underlying simplicial mesh. We discuss the construction and implementation of the approach, explaining how to modify shape derivate formulas to compute consistent shape gradient approximations using discontinuous Galerkin methods, and how to recover dG functions into smoother ones. Numerical experiments on unconstrained and PDE-constrained test cases evidence the good properties of the proposed methodology., Comment: 20 pages, 8 figures

This paper is concerned with developing and analyzing two novel implicit temporal discretization methods for the stochastic semilinear wave equations with multiplicative noise. The proposed methods are natural extensions of well-known time-discrete schemes for deterministic wave equations, hence, they are easy to implement. It is proved that both methods are energy-stable. Moreover, the first method is shown to converge with the linear order in the energy norm, while the second method converges with the $\mathcal{O}(\tau^{\frac32})$ order in the $L^2$-norm, which is optimal with respect to the time regularity of the solution to the underlying stochastic PDE. The convergence analyses of both methods, which are different and quite involved, require some novel numerical techniques to overcome difficulties caused by the nonlinear noise term and the interplay between nonlinear drift and diffusion. Numerical experiments are provided to validate the sharpness of the theoretical error estimate results., Comment: 28 pages, 0 figure, 3 tables

Differential equations posed on quadratic matrix Lie groups arise in the context of classical mechanics and quantum dynamical systems. Lie group numerical integrators preserve the constants of motions defining the Lie group. Thus, they respect important physical laws of the dynamical system, such as unitarity and energy conservation in the context of quantum dynamical systems, for instance. In this article we develop a high-order commutator free Lie group integrator for non-autonomous differential equations evolving on quadratic Lie groups. Instead of matrix exponentials, which are expensive to evaluate and need to be approximated by appropriate rational functions in order to preserve the Lie group structure, the proposed method is obtained as a composition of Cayley transforms which naturally respect the structure of quadratic Lie groups while being computationally efficient to evaluate. Unlike Cayley-Magnus methods the method is also free from nested matrix commutators., Comment: 18 pages, 3 figures, submitted for a review in SINUM (SIAM Journal on Numerical Analysis)

We give a complete expansion, at any accuracy order, for the iterated convolution of a complex valued integrable sequence in one space dimension. The remainders are estimated sharply with generalized Gaussian bounds. The result applies in probability theory for random walks as well as in numerical analysis for studying the large time behavior of numerical schemes.

Surface-subsurface flow models for hydrological applications solve a coupled multiphysics problem. This usually consists of some form of the Richards and shallow water equations. A typical setup couples these two nonlinear partial differential equations in a partitioned approach via boundary conditions. Full interaction between the subsolvers is ensured by an iterative coupling procedure. This can be accelerated using relaxation. In this paper, we apply continuous and fully discrete linear analysis techniques to study an idealized, linear, 1D-0D version of a surface-subsurface model. These result in explicit expressions for the convergence factor and an optimal relaxation parameter, depending on material and discretization parameters. We test our analysis results numerically for fully nonlinear 2D-1D experiments based on existing benchmark problems. The linear analysis can explain fast convergence of iterations observed in practice for different materials and test cases, even though we are not able to capture various nonlinear effects., Comment: 29 pages, 10 figures

We present here the classical Schwarz method with a time domain decomposition applied to unconstrained parabolic optimal control problems. Unlike Dirichlet-Neumann and Neumann-Neumann algorithms, we find different properties based on the forward-backward structure of the optimality system. Variants can be found using only Dirichlet and Neumann transmission conditions. Some of these variants are only good smoothers, while others could lead to efficient solvers., Comment: 28th International Conference on Domain Decomposition Methods (DD28)

This work presents a matrix-free finite element solver for finite-strain elasticity adopting an hp-multigrid preconditioner. Compared to classical algorithms relying on a global sparse matrix, matrix-free solution strategies significantly reduce memory traffic by on-the-fly evaluation of the finite element integrals. Following this approach in the context of finite-strain elasticity, the precise statement of the final weak form is crucial for performance, and it is not clear a priori whether to choose problem formulations in the material or spatial domain. With a focus on hyperelastic solids in biomechanics, the arithmetic costs to evaluate the material law at each quadrature point might favor an evaluation strategy where some quantities are precomputed in each Newton iteration and reused in the Krylov solver for the linearized problem. Hence, we discuss storage strategies to balance the compute load against memory access in compressible and incompressible neo-Hookean models and an anisotropic tissue model. Additionally, numerical stability becomes increasingly important using lower/mixed-precision ingredients and approximate preconditioners to better utilize modern hardware architectures. Application of the presented method to a patient-specific geometry of an iliac bifurcation shows significant speed-ups, especially for higher polynomial degrees, when compared to alternative approaches with matrix-based geometric or black-box algebraic multigrid preconditioners.

The modeling of electric machines and power transformers typically involves systems of nonlinear magnetostatics or -quasistatics, and their efficient and accurate simulation is required for the reliable design, control, and optimization of such devices. We study the numerical solution of the vector potential formulation of nonlinear magnetostatics by means of higher-order finite element methods. Numerical quadrature is used for the efficient handling of the nonlinearities and domain mappings are employed for the consideration of curved boundaries. The existence of a unique solution is proven on the continuous and discrete level and a full convergence analysis of the resulting finite element schemes is presented indicating order optimal convergence rates under appropriate smoothness assumptions. For the solution of the nonlinear discretized problems, we consider a Newton method with line search for which we establish global linear convergence with convergence rates that are independent of the discretization parameters. We further prove local quadratic convergence in a mesh-dependent neighborhood of the solution which becomes effective when high accuracy of the nonlinear solver is demanded. The assumptions required for our analysis cover inhomogeneous, nonlinear, and anisotropic materials, which may arise in typical applications, including the presence of permanent magnets. The theoretical results are illustrated by numerical tests for some typical benchmark problems.

We present results of numerical experiments for neural networks with stochastic gradient-based optimization with adaptive momentum. This widely applied optimization has proved convergence and practical efficiency, but for long-run training becomes numerically unstable. We show that numerical artifacts are observable not only for large-scale models and finally lead to divergence also for case of shallow narrow networks. We argue this theory by experiments with more than 1600 neural networks trained for 50000 epochs. Local observations show presence of the same behavior of network parameters in both stable and unstable training segments. Geometrical behavior of parameters forms double twisted spirals in the parameter space and is caused by alternating of numerical perturbations with next relaxation oscillations in values for 1st and 2nd momentum.

Ordinary differential equations (ODE) are a popular tool to model the spread of infectious diseases, yet they implicitly assume an exponential distribution to describe the flow from one infection state to another. However, scientific experience yields more plausible distributions where the likelihood of disease progression or recovery changes accordingly with the duration spent in a particular state of the disease. Furthermore, transmission dynamics depend heavily on the infectiousness of individuals. The corresponding nonlinear variation with the time individuals have already spent in an infectious state requires more realistic models. The previously mentioned items are particularly crucial when modeling dynamics at change points such as the implementation of nonpharmaceutical interventions. In order to capture these aspects and to enhance the accuracy of simulations, integro-differential equations (IDE) can be used. In this paper, we propose a generalized model based on integro-differential equations with eight infection states. The model allows for variable stay time distributions and generalizes the concept of ODE-based models as well as IDE-based age-of-infection models. In this, we include particular infection states for severe and critical cases to allow for surveillance of the clinical sector, avoiding bottlenecks and overloads in critical epidemic situations. We will extend a recently introduced nonstandard numerical scheme to solve a simpler IDE-based model. This scheme is adapted to our more advanced model and we prove important mathematical and biological properties for the numerical solution. Furthermore, we validate our approach numerically by demonstrating the convergence rate. Eventually, we also show that our novel model is intrinsically capable of better assessing disease dynamics upon the introduction of nonpharmaceutical interventions.

Reconstructing cardiac electrical activity from body surface electric potential measurements results in the severely ill-posed inverse problem in electrocardiography. Many different regularization approaches have been proposed to improve numerical results and provide unique results. This work presents a novel approach for reconstructing the epicardial potential from body surface potential maps based on a space-time total variation-type regularization using finite elements, where a first-order primal-dual algorithm solves the underlying convex optimization problem. In several numerical experiments, the superior performance of this method and the benefit of space-time regularization for the reconstruction of epicardial potential on two-dimensional torso data and a three-dimensional rabbit heart compared to state-of-the-art methods are demonstrated.

A rigorous Bayesian formulation of the inverse doping profile problem in infinite dimensions for a stationary linearized unipolar drift-diffusion model for semiconductor devices is given. The goal is to estimate the posterior probability distribution of the doping profile and to compute its posterior mean. This allows for the reconstruction of the doping profile from voltage-current measurements. The well-posedness of the Bayesian inverse problem is shown by proving boundedness and continuity properties of the semiconductor model with respect to the unknown parameter. A preconditioned Crank-Nicolson Markov chain Monte-Carlo method for the Bayesian estimation of the doping profile, using a physics-informed prior model, is proposed. The numerical results for a two-dimensional diode illustrate the efficiency of the proposed approach.

In scattered data approximation, the span of a finite number of translates of a chosen radial basis function is used as approximation space and the basis of translates is used for representing the approximate. However, this natural choice is by no means mandatory and different choices, like, for example, the Lagrange basis, are possible and might offer additional features. In this article, we discuss different alternatives together with their canonical duals. We study a localized version of the Lagrange basis, localized orthogonal bases, such as the Newton basis, and multiresolution versions thereof, constructed by means of samplets. We argue that the choice of orthogonal bases is particularly useful as they lead to symmetric preconditioners. All bases under consideration are compared numerically to illustrate their feasibility for scattered data approximation. We provide benchmark experiments in two spatial dimensions and consider the reconstruction of an implicit surface as a relevant application from computer graphics.

In nature, many complex multi-physics coupling problems exhibit significant diffusivity inhomogeneity, where one process occurs several orders of magnitude faster than others in temporal. Simulating rapid diffusion alongside slower processes demands intensive computational resources due to the necessity for small time steps. To address these computational challenges, we have developed an efficient numerical solver named Finite Difference informed Random Walker (FDiRW). In this study, we propose a GPU-accelerated, mixed-precision configuration for the FDiRW solver to maximize efficiency through GPU multi-threaded parallel computation and lower precision computation. Numerical evaluation results reveal that the proposed GPU-accelerated mixed-precision FDiRW solver can achieve a 117X speedup over the CPU baseline, while an additional 1.75X speedup by employing lower precision GPU computation. Notably, for large model sizes, the GPU-accelerated mixed-precision FDiRW solver demonstrates strong scaling with the number of nodes used in simulation. When simulating radionuclide absorption processes by porous wasteform particles with a medium-sized model of 192x192x192, this approach reduces the total computational time to 10 minutes, enabling the simulation of larger systems with strongly inhomogeneous diffusivity., Comment: 24 pages, 16 figures, 2 tables

We consider the use of multipreconditioning, which allows for multiple preconditioners to be applied in parallel, on high-frequency Helmholtz problems. Typical applications present challenging sparse linear systems which are complex non-Hermitian and, due to the pollution effect, either very large or else still large but under-resolved in terms of the physics. These factors make finding general purpose, efficient and scalable solvers difficult and no one approach has become the clear method of choice. In this work we take inspiration from domain decomposition strategies known as sweeping methods, which have gained notable interest for their ability to yield nearly-linear asymptotic complexity and which can also be favourable for high-frequency problems. While successful approaches exist, such as those based on higher-order interface conditions, perfectly matched layers (PMLs), or complex tracking of wave fronts, they can often be quite involved or tedious to implement. We investigate here the use of simple sweeping techniques applied in different directions which can then be incorporated in parallel into a multipreconditioned GMRES strategy. Preliminary numerical results on a two-dimensional benchmark problem will demonstrate the potential of this approach.

This paper introduces an auto-stabilized weak Galerkin (WG) finite element method with a built-in stabilizer for Poisson equations. By utilizing bubble functions as a key analytical tool, our method extends to both convex and non-convex elements in finite element partitions, marking a significant advancement over existing stabilizer-free WG methods. It overcomes the restrictive conditions of previous approaches and is applicable in any dimension $d$, offering substantial advantages. The proposed method maintains a simple, symmetric, and positive definite structure. These benefits are evidenced by optimal order error estimates in both discrete $H^1$ and $L^2$ norms, highlighting the effectiveness and accuracy of our WG method for practical applications., Comment: 19 pages

In this paper, we propose a parameterized shift-splitting (PSS) preconditioner, along with its two relaxed variants to solve the regularized three-by-three block saddle point problem (RTBSPP). The convergence of the associated PSS iteration method is analyzed, and sufficient conditions on the parameters for convergence are established. Spectral analyses are performed to derive sharp bounds for the eigenvalues of the preconditioned matrices. Numerical experiments based on examples arising from the PDE-constrained optimization problems are performed to demonstrate the effectiveness and robustness of the proposed preconditioners compared with existing state-of-the-art preconditioners.

Computations in high-dimensional spaces can often be realized only approximately, using a certain number of projections onto lower dimensional subspaces or sampling from distributions. In this paper, we are interested in pairs of real-valued functions $(F,f)$ on $[0,\infty)$ that are related by the projection/slicing formula $F (\| x \|) = \mathbb E_{\xi} \big[ f \big(|\langle x,\xi \rangle| \big) \big]$ for $x\in\mathbb R^d$, where the expectation value is taken over uniformly distributed direction in $\mathbb R^d$. While it is known that $F$ can be obtained from $f$ by an Abel-like integral formula, we construct conversely $f$ from given $F$ using their Fourier transforms. First, we consider the relation between $F$ and $f$ for radial functions $F(\| \cdot\| )$ that are Fourier transforms of $L^1$ functions. Besides $d$- and one-dimensional Fourier transforms, it relies on a rotation operator, an averaging operator and a multiplication operator to manage the walk from $d$ to one dimension in the Fourier space. Then, we generalize the results to tempered distributions, where we are mainly interested in radial regular tempered distributions. Based on Bochner's theorem, this includes positive definite functions $F(\| \cdot\| )$., Comment: 33 pages, 2 figures

Significant research efforts have been dedicated recently to explore the structured backward error (BE) for saddle point problems (SPPs). However, these investigations overlook the inherent sparsity pattern of the coefficient matrix of the SPP. Moreover, the existing techniques are not applicable when the block matrices have Circulant, Toeplitz, or symmetric-Toeplitz structures and do not even provide structure preserving minimal perturbation matrices for which the BE is attained. To overcome these limitations, we investigate the structured BEs of SPPs when the perturbation matrices exploit the sparsity pattern as well as Circulant, Toeplitz, and symmetric-Toeplitz structures. Furthermore, we construct minimal perturbation matrices that preserve the sparsity pattern and the aforementioned structures. Two applications of the obtained results are discussed: (i) deriving structured BEs for the weighted regularized least squares problem and (ii) finding an effective and reliable stopping criterion for numerical algorithms. Finally, numerical experiments validate our findings, showcasing the utility of the obtained structured BEs in assessing the strong stability of numerical algorithms

Deep learning has become a popular tool across many scientific fields, including the study of differential equations, particularly partial differential equations. This work introduces the basic principles of deep learning and the Deep Galerkin method, which uses deep neural networks to solve differential equations. This primer aims to provide technical and practical insights into the Deep Galerkin method and its implementation. We demonstrate how to solve the one-dimensional heat equation step-by-step. We also show how to apply the Deep Galerkin method to solve systems of ordinary differential equations and integral equations, such as the Fredholm of the second kind. Additionally, we provide code snippets within the text and the complete source code on Github. The examples are designed so that one can run them on a simple computer without needing a GPU., Comment: 32 pages, 12 figures, primer (tutorial)

In [1], we calibrated a one-factor Cheyette SLV model with a local volatility that is linear in the benchmark forward rate and an uncorrelated CIR stochastic variance to 3M caplets of various maturities. While caplet smiles for many maturities could be reasonably well calibrated across the range of strikes, for instance the 1Y maturity could not be calibrated well across that entire range of strikes. Here, we study whether models with alternative local volatility terms and/or alternative stochastic volatility or variance models can calibrate the 1Y caplet smile better across the strike range better than the model studied in [1]. This is made possible and feasible by the generic simulation, pricing, and calibration frameworks introduced in [1] and some new frameworks presented in this paper. We find that some model settings calibrate well to the 1Y smile across the strike range under study. In particular, a model setting with a local volatility that is piece-wise linear in the benchmark forward rate together with an uncorrelated CIR stochastic variance and one with a local volatility that is linear in the benchmark rate together with a correlated lognormal stochastic volatility with quadratic drift (QDLNSV) as in [2] calibrate well. We discuss why the later might be a preferable model. [1] Arun Kumar Polala and Bernhard Hientzsch. Parametric differential machine learning for pricing and calibration. arXiv preprint arXiv:2302.06682 , 2023. [2] Artur Sepp and Parviz Rakhmonov. A Robust Stochastic Volatility Model for Interest Rate Dynamics. Risk Magazine, 2023

In this paper we propose a semi-Lagrangian discontinuous Galerkin solver for the simulation of the scrape off layer for an electron-ion plasma. We use a time adaptive velocity space to deal with fast particles leaving the computational domain, a block structured mesh to resolve the sharp gradient in the plasma sheath, and limiters to avoid oscillations in the density function. In particular, we propose a limiter that can be computed directly from the information used in the semi-Lagrangian discontinuous Galerkin advection step. This limiter is particularly efficient on graphic processing units (GPUs) and compares favorable with limiters from the literature. We provide numerical results for a set of benchmark problems and compare different limiting strategies.

Computing the proximal operator of the $\ell_\infty$ norm, $\textbf{prox}_{\alpha ||\cdot||_\infty}(\mathbf{x})$, generally requires a sort of the input data, or at least a partial sort similar to quicksort. In order to avoid using a sort, we present an $O(m)$ approximation of $\textbf{prox}_{\alpha ||\cdot||_\infty}(\mathbf{x})$ using a neural network. A novel aspect of the network is that it is able to accept vectors of varying lengths due to a feature selection process that uses moments of the input data. We present results on the accuracy of the approximation, feature importance, and computational efficiency of the approach. We show that the network outperforms a "vanilla neural network" that does not use feature selection. We also present an algorithm with corresponding theory to calculate $\textbf{prox}_{\alpha ||\cdot||_\infty}(\mathbf{x})$ exactly, relate it to the Moreau decomposition, and compare its computational efficiency to that of the approximation., Comment: 30 pages, 5 figures, 6 tables

Zeitlin's model is a spatial discretization for the 2-D Euler equations on the flat 2-torus or the 2-sphere. Contrary to other discretizations, it preserves the underlying geometric structure, namely that the Euler equations describe Riemannian geodesics on a Lie group. Here we show how to extend Zeitlin's approach to the axisymmetric Euler equations on the 3-sphere. It is the first discretization of the 3-D Euler equations that fully preserves the geometric structure. Thus, this finite-dimensional model admits Riemannian curvature and Jacobi equations, which are discussed., Comment: 24 pages, 4 figures

We introduce and compare two domain decomposition based numerical methods, namely the Dirichlet-Neumann and Neumann-Neumann Waveform Relaxation methods (DNWR and NNWR respectively), tailored for solving partial differential equations (PDEs) incorporating time delay. Time delay phenomena frequently arise in various real-world systems, making their accurate modeling and simulation crucial for understanding and prediction. We consider a series of model problems, ranging from Parabolic, Hyperbolic to Neutral PDEs with time delay and apply the iterative techniques DNWR and NNWR for solving in parallel. We present the theoretical foundations, numerical implementation, and comparative performance analysis of these two methods. Through numerical experiments and simulations, we explore their convergence properties, computational efficiency, and applicability to various types of PDEs with time delay.

In the high-dimensional setting, Gaussian mixture kernel density estimates become increasingly suboptimal. In this work we aim to show that it is practical to instead use the optimal multivariate Epanechnikov kernel. We make use of this optimal Epanechnikov mixture kernel density estimate for the sequential filtering scenario through what we term the ensemble Epanechnikov mixture filter (EnEMF). We provide a practical implementation of the EnEMF that is as cost efficient as the comparable ensemble Gaussian mixture filter. We show on a static example that the EnEMF is robust to growth in dimension, and also that the EnEMF has a significant reduction in error per particle on the 40-variable Lorenz '96 system.

Whilst Paul de Casteljau is now famous for his fundamental algorithm of curve and surface approximation, little is known about his other findings. This article offers an insight into his results in geometry, algebra and number theory. Related to geometry, his classical algorithm is reviewed as an index reduction of a polar form. This idea is used to show de Casteljau's algebraic way of smoothing, which long went unnoticed. We will also see an analytic polar form and its use in finding the intersection of two curves. The article summarises unpublished material on metric geometry. It includes theoretical advances, e.g., the 14-point strophoid or a way to link Apollonian circles with confocal conics, and also practical applications such as a recurrence for conjugate mirrors in geometric optics. A view on regular polygons leads to an approximation of their diagonals by golden matrices, a generalisation of the golden ratio. Relevant algebraic findings include matrix quaternions (and anti-quaternions) and their link with Lorentz' equations. De Casteljau generalised the Euclidean algorithm and developed an automated method for approximating the roots of a class of polynomial equations. His contributions to number theory not only include aspects on the sum of four squares as in quaternions, but also a view on a particular sum of three cubes. After a review of a complete quadrilateral in a heptagon and its angles, the paper concludes with a summary of de Casteljau's key achievements. The article contains a comprehensive bibliography of de Casteljau's works, including previously unpublished material., Comment: 56 pages, 19 sections, 30 figures - part of a special issue on Paul de Casteljau

In this paper, we compute the expected logarithmic energy of solutions to the polynomial eigenvalue problem for random matrices. We generalize some known results for the Shub-Smale polynomials, and the spherical ensemble. These two processes are the two extremal particular cases of the polynomial eigenvalue problem, and we prove that the logarithmic energy lies between these two cases. In particular, the roots of the Shub-Smale polynomials are the ones with the lowest logarithmic energy of the family.

Let $\Omega\subset \mathbb{R}^d$ be a bounded domain. We consider the problem of how efficiently shallow neural networks with the ReLU$^k$ activation function can approximate functions from Sobolev spaces $W^s(L_p(\Omega))$ with error measured in the $L_q(\Omega)$-norm. Utilizing the Radon transform and recent results from discrepancy theory, we provide a simple proof of nearly optimal approximation rates in a variety of cases, including when $q\leq p$, $p\geq 2$, and $s \leq k + (d+1)/2$. The rates we derive are optimal up to logarithmic factors, and significantly generalize existing results. An interesting consequence is that the adaptivity of shallow ReLU$^k$ neural networks enables them to obtain optimal approximation rates for smoothness up to order $s = k + (d+1)/2$, even though they represent piecewise polynomials of fixed degree $k$.

A neural network architecture is presented that exploits the multilevel properties of high-dimensional parameter-dependent partial differential equations, enabling an efficient approximation of parameter-to-solution maps, rivaling best-in-class methods such as low-rank tensor regression in terms of accuracy and complexity. The neural network is trained with data on adaptively refined finite element meshes, thus reducing data complexity significantly. Error control is achieved by using a reliable finite element a posteriori error estimator, which is also provided as input to the neural network. The proposed U-Net architecture with CNN layers mimics a classical finite element multigrid algorithm. It can be shown that the CNN efficiently approximates all operations required by the solver, including the evaluation of the residual-based error estimator. In the CNN, a culling mask set-up according to the local corrections due to refinement on each mesh level reduces the overall complexity, allowing the network optimization with localized fine-scale finite element data. A complete convergence and complexity analysis is carried out for the adaptive multilevel scheme, which differs in several aspects from previous non-adaptive multilevel CNN. Moreover, numerical experiments with common benchmark problems from Uncertainty Quantification illustrate the practical performance of the architecture.

In this paper we propose a numerical procedure for the quantification of uncertainties in wave-structure interaction. We utilise the smoothed particle hydrodynamics (SPH) scheme for modelling the wave mechanics, coupled one-way with a finite element method (FEM) for the structural response. Physical wave flumes are extensively used in the study of hydrodynamics especially in wave-structure interaction (WSI) and prediction of forces to near-shore structures in disaster mitigation and offshore structures in the oil and gas, and more recently renewable energy sector. Over the years, numerical wave flumes have been developed extensively to enable the modelling of complex wave-structure interaction. However, most of these studies are deterministic and limited to using either simple flexible beams or rigid monolithic structures to model the structural part in the WSI. Additionally, uncertainties are commonly observed in both wave and structural parameters and need to be accounted for. This work presents a numerical framework to enable uncertainty quantification for wave-structure interaction problems in terms of the forces experienced by the structure. A one-way coupling between SPH with the FEM and uncertainty quantification (UQ) is proposed. We employ the so-called Tokyo wave flume geometry, which has potential for future surrogate modelling in WSI. The developed model is validated using numerical and experimental results from the literature and is used to demonstrate the prediction of probabilistic responses of structures under breaking and non-breaking wave scenarios.

We consider solving a probably ill-conditioned linear operator equation, where the operator is not modeled by physical laws but is specified via training pairs (consisting of images and data) of the input-output relation of the operator. We derive a stable method for computing the operator, which consists of first a Gram-Schmidt orthonormalization of images and a principal component analysis of the data. This two-step algorithm provides a spectral decomposition of the linear operator. Moreover, we show that both Gram-Schmidt and principal component analysis can be written as a deep neural network, which relates this procedure to de-and encoder networks. Therefore, we call the two-step algorithm a linear algebra network. Finally, we provide numerical simulations showing the strategy is feasible for reconstructing spectral functions and for solving operator equations without explicitly exploiting the physical model.

We present an approach to construct efficient sparse summation-by-parts (SBP) operators on triangles and tetrahedra with a tensor-product structure. The operators are constructed by splitting the simplices into quadrilateral or hexahedral subdomains, mapping tensor-product SBP operators onto the subdomains, and assembling back using a continuous-Galerkin-type procedure. These tensor-product split-simplex operators do not have repeated degrees of freedom at the interior interfaces between the split subdomains. Furthermore, they satisfy the SBP property by construction, leading to stable discretizations. The accuracy and sparsity of the operators substantially enhance the efficiency of SBP discretizations on simplicial meshes. The sparsity is particularly important for entropy-stable discretizations based on two-point flux functions, as it reduces the number of two-point flux computations. We demonstrate through numerical experiments that the operators exhibit efficiency surpassing that of the existing dense multidimensional SBP operators by more than an order of magnitude in many cases. This superiority is evident in both accuracy per degree of freedom and computational time required to achieve a specified error threshold., Comment: 28 pages, 15 figures

In this paper, we discuss the second-order finite element method (FEM) and finite difference method (FDM) for numerically solving elliptic cross-interface problems characterized by vertical and horizontal straight lines, piecewise constant coefficients, two homogeneous jump conditions, continuous source terms, and Dirichlet boundary conditions. For brevity, we consider a 2D simplified version where the intersection points of the interface lines coincide with grid points in uniform Cartesian grids. Our findings reveal interesting and important results: (1) When the coefficient functions exhibit either high jumps with low-frequency oscillations or low jumps with high-frequency oscillations, the finite element method and finite difference method yield similar numerical solutions. (2) However, when the interface problems involve high-contrast and high-frequency coefficient functions, the numerical solutions obtained from the finite element and finite difference methods differ significantly. Given that the widely studied SPE10 benchmark problem (see https://www.spe.org/web/csp/datasets/set02.htm) typically involves high-contrast and high-frequency permeability due to varying geological layers in porous media, this phenomenon warrants attention. Furthermore, this observation is particularly important for developing multiscale methods, as reference solutions for these methods are usually obtained using the standard second-order finite element method with a fine mesh, and analytical solutions are not available. We provide sufficient details to enable replication of our numerical results, and the implementation is straightforward. This simplicity ensures that readers can easily confirm the validity of our findings., Comment: 19 pages, 12 figures

Nonsmooth composite optimization problems under uncertainty are prevalent in various scientific and engineering applications. We consider risk-neutral composite optimal control problems, where the objective function is the sum of a potentially nonconvex expectation function and a nonsmooth convex function. To approximate the risk-neutral optimization problems, we use a Monte Carlo sample-based approach, study its asymptotic consistency, and derive nonasymptotic sample size estimates. Our analyses leverage problem structure commonly encountered in PDE-constrained optimization problems, including compact embeddings and growth conditions. We apply our findings to bang-bang-type optimal control problems and propose the use of a conditional gradient method to solve them effectively. We present numerical illustrations., Comment: 21 pages

This paper addresses the convergence analysis of a variant of the LevenbergMarquardt method (LMM) designed for nonlinear least-squares problems with non-zero residue. This variant, called LMM with Singular Scaling (LMMSS), allows the LMM scaling matrix to be singular, encompassing a broader class of regularizers, which has proven useful in certain applications. In order to establish local convergence results, a careful choice of the LMM parameter is made based on the gradient linearization error, dictated by the nonlinearity and size of the residual. Under completeness and local error bound assumptions we prove that the distance from an iterate to the set of stationary points goes to zero superlinearly and that the iterative sequence converges. Furthermore, we also study a globalized version of the method obtained using linesearch and prove that any limit point of the generated sequence is stationary. Some examples are provided to illustrate our theoretical results.

Interest in communication-avoiding orthogonalization schemes for high-performance computing has been growing recently. This manuscript addresses open questions about the numerical stability of various block classical Gram-Schmidt variants that have been proposed in the past few years. An abstract framework is employed, the flexibility of which allows for new rigorous bounds on the loss of orthogonality in these variants. We first analyze a generalization of (reorthogonalized) block classical Gram-Schmidt and show that a "strong" intrablock orthogonalization routine is only needed for the very first block in order to maintain orthogonality on the level of the unit roundoff. Then, using this variant, which has four synchronization points per block column, we remove the synchronization points one at a time and analyze how each alteration affects the stability of the resulting method. Our analysis shows that the variant requiring only one synchronization per block column cannot be guaranteed to be stable in practice, as stability begins to degrade with the first reduction of synchronization points. Our analysis of block methods also provides new theoretical results for the single-column case. In particular, it is proven that DCGS2 from [Bielich, D. et al. Par. Comput. 112 (2022)] and CGS-2 from [\'{S}wirydowicz, K. et al, Num. Lin. Alg. Appl. 28 (2021)] are as stable as Householder QR. Numerical examples from the BlockStab toolbox are included throughout, to help compare variants and illustrate the effects of different choices of intraorthogonalization subroutines.

We extend the divergence preserving cut finite element method presented in [T. Frachon, P. Hansbo, E. Nilsson, S. Zahedi, SIAM J. Sci. Comput., 46 (2024)] for the Darcy interface problem to unfitted outer boundaries. We impose essential boundary conditions on unfitted meshes with a stabilized Lagrange multiplier method. The stabilization term for the Lagrange multiplier is important for stability but it may perturb the approximate solution at the boundary. We study different stabilization terms from cut finite element discretizations of surface partial differential equations and trace finite element methods. To reduce the perturbation we use a Lagrange multiplier space of higher polynomial degree compared to previous work on unfitted discretizations. We propose a symmetric method that results in 1) optimal rates of convergence for the approximate velocity and pressure; 2) well-posed linear systems where the condition number of the system matrix scales as for fitted finite element discretizations; 3) optimal approximation of the divergence with pointwise divergence-free approximations of solenoidal velocity fields. The three properties are proven to hold for the lowest order discretization and numerical experiments indicate that these properties continue to hold also when higher order elements are used.

In a recently published article by G. Ambrus et al. a new upper bound for the density of an unit avoiding, periodic set is given as $0.2470$, the first upper bound $< 1/4$. A construction of Croft 1967 gave a lower bound $\delta_C = 0.22936$ for the density. To this date, no better construction with a higher lower bound has been given. In this article I give a construction planar sets with a higher density than Croft's tortoises. No explicit value for this density is given, it's just shown that Croft's density is a local minima of the density of a here constructed 1-parameter family of planar sets. So the densities are $> \delta_C$., Comment: 9 pages, comments welcome

Algorithms for numerical tasks in finite precision simultaneously seek to minimize the number of floating point operations performed, and also the number of bits of precision required by each floating point operation. This paper presents an algorithm for Hermitian diagonalization requiring only $\lg(1/\varepsilon)+O(\log(n)+\log\log(1/\varepsilon))$ bits of precision where $n$ is the size of the input matrix and $\varepsilon$ is the target error. Furthermore, it runs in near matrix multiplication time. In the general setting, the first complete analysis of the stability of a near matrix multiplication time algorithm for diagonalization is that of Banks et al. [BGVKS20]. They exhibit an algorithm for diagonalizing an arbitrary matrix up to $\varepsilon$ backward error using only $O(\log^4(n/\varepsilon)\log(n))$ bits of precision. This work focuses on the Hermitian setting, where we determine a dramatically improved bound on the number of bits needed. In particular, the result is close to providing a practical bound. The exact bit count depends on the specific implementation of matrix multiplication and QR decomposition one wishes to use, but if one uses suitable $O(n^3)$-time implementations, then for $\varepsilon=10^{-15},n=4000$, we show 92 bits of precision suffice (and 59 are necessary). By comparison, the same parameters in [BGVKS20] does not even show that 682,916,525,000 bits suffice.

We consider delay differential equations with a polynomially distributed delay. We derive an equivalent system of delay differential equations, which includes just two discrete delays. The stability of the equivalent system and its stationary solutions are investigated. Alternatively, a Gaussian quadrature generates a discretisation of the integral, which describes the distributed delay in the original delay differential equation. This technique yields an approximate differential equation with multiple discrete delays. We present results of numerical computations, where initial value problems of the differential equations are solved. Therein, the two approaches are compared., Comment: 15 pages, 8 figures

This paper describes a simplified model of an electric circuit with a DC-DC converter and a PID-regulator as a system of integral differential equations with an identically singular matrix multiplying the higher derivative of the desired vector-function. We use theoretical results on integral and differential equations and their systems to prove solvability of such a model and analyze its stability., Comment: submitted to COIA-2024. arXiv admin note: text overlap with arXiv:2408.06045

This survey examines the broad suite of methods and models for combining machine learning with physics knowledge for prediction and forecast, with a focus on partial differential equations. These methods have attracted significant interest due to their potential impact on advancing scientific research and industrial practices by improving predictive models with small- or large-scale datasets and expressive predictive models with useful inductive biases. The survey has two parts. The first considers incorporating physics knowledge on an architectural level through objective functions, structured predictive models, and data augmentation. The second considers data as physics knowledge, which motivates looking at multi-task, meta, and contextual learning as an alternative approach to incorporating physics knowledge in a data-driven fashion. Finally, we also provide an industrial perspective on the application of these methods and a survey of the open-source ecosystem for physics-informed machine learning., Comment: 56 pages, 8 figures, 2 tables

We present and analyze a framework for designing symplectic neural networks (SympNets) based on geometric integrators for Hamiltonian differential equations. The SympNets are universal approximators in the space of Hamiltonian diffeomorphisms, interpretable and have a non-vanishing gradient property. We also give a representation theory for linear systems, meaning the proposed P-SympNets can exactly parameterize any symplectic map corresponding to quadratic Hamiltonians. Extensive numerical tests demonstrate increased expressiveness and accuracy -- often several orders of magnitude better -- for lower training cost over existing architectures. Lastly, we show how to perform symbolic Hamiltonian regression with SympNets for polynomial systems using backward error analysis., Comment: 33 pages including appendices but not references, 7 figures

Using generating functions, we are proposing a unified approach to produce explicit formulas, which count the number of nodes in Smolyak grids based on various univariate quadrature or interpolation rules. Our approach yields, for instance, a new formula for the cardinality of a Smolyak grid, which is based on Chebyshev nodes of the first kind and it allows to recover certain counting-formulas previously found by Bungartz-Griebel, Kaarnioja, M\"uller-Gronbach, Novak-Ritter and Ullrich., Comment: 10 pages, 1 figure

This paper considers one of the fundamental parallel-in-time methods for the solution of ordinary differential equations, Parareal, and extends it by adopting a neural network as a coarse propagator. We provide a theoretical analysis of the convergence properties of the proposed algorithm and show its effectiveness for several examples, including Lorenz and Burgers' equations. In our numerical simulations, we further specialize the underpinning neural architecture to Random Projection Neural Networks (RPNNs), a 2-layer neural network where the first layer weights are drawn at random rather than optimized. This restriction substantially increases the efficiency of fitting RPNN's weights in comparison to a standard feedforward network without negatively impacting the accuracy, as demonstrated in the SIR system example.

In this paper, we derive a thermodynamically consistent non-isothermal diffuse interface model for phase transition and interface evolution involving heat transfer. This model is constructed by integrating concepts from classical irreversible thermodynamics with the energetic variational approach. By making specific assumptions about the kinematics of the temperature, we derive a non-isothermal Allen-Cahn equation. Through both asymptotic analysis and numerical simulations, we demonstrate that in the sharp interface limit, the non-isothermal Allen-Cahn equation converges to a two-phase Stefan-like problem, under a certain scaling of the melting/freezing energy. In this regime, the motion of the liquid-solid interface and the temperature interface coincide and are governed by the mean curvature, at least for a short time. The result provides a justification for the classical Stefan problem within a certain physical regime., Comment: 23 pages, 5 figures

Solving real-world nonlinear semiconductor device problems modeled by the drift-diffusion equations coupled with the Poisson equation (also known as the Poisson-Nernst-Planck equations) necessitates an accurate and efficient numerical scheme which can avoid non-physical oscillations even for problems with extremely sharp doping profiles. In this paper, we propose a flexible and high-order hybridizable discontinuous Galerkin (HDG) scheme with harmonic averaging (HA) technique to tackle these challenges. The proposed HDG-HA scheme combines the robustness of finite volume Scharfetter-Gummel (FVSG) method with the high-order accuracy and $hp$-flexibility offered by the locally conservative HDG scheme. The coupled Poisson equation and two drift-diffusion equations are simultaneously solved by the Newton method. Indicators based on the gradient of net doping $N$ and solution variables are proposed to switch between cells with HA technique and high-order conventional HDG cells, utilizing the flexibility of HDG scheme. Numerical results suggest that the proposed scheme does not exhibit oscillations or convergence issues, even when applied to heavily doped and sharp PN-junctions. Devices with circular junctions and realistic doping profiles are simulated in two dimensions, qualifying this scheme for practical simulation of real-world semiconductor devices.

We consider trajectory optimal control problems in which parameter uncertainty limits the applicability of control trajectories computed prior to travel. Hence, efficient trajectory adjustment is needed to ensure successful travel. However, it is often prohibitive or impossible to recalculate the optimal control in-transit due to strict time constraints or limited onboard computing resources. Thus, we propose a framework for quick and accurate trajectory approximations by using post-optimality sensitivity information. This allows the reduction of uncertain parameter space and an instantaneous approximation of the new optimal controller while using sensitivity data computed and stored pretransit.