Date: Thu, 9 May 2024 | Total: 17
We propose a computationally efficient and systematically convergent approach for elastodynamics simulations. We recast the second-order dynamical equation of elastodynamics into an equivalent first-order system of coupled equations, so as to express the solution in the form of a Magnus expansion. With any spatial discretization, it entails computing the exponential of a matrix acting upon a vector. We employ an adaptive Krylov subspace approach to inexpensively and and accurately evaluate the action of the exponential matrix on a vector. In particular, we use an apriori error estimate to predict the optimal Kyrlov subspace size required for each time-step size. We show that the Magnus expansion truncated after its first term provides quadratic and superquadratic convergence in the time-step for nonlinear and linear elastodynamics, respectively. We demonstrate the accuracy and efficiency of the proposed method for one linear (linear cantilever beam) and three nonlinear (nonlinear cantilever beam, soft tissue elastomer, and hyperelastic rubber) benchmark systems. For a desired accuracy in energy, displacement, and velocity, our method allows for $10-100\times$ larger time-steps than conventional time-marching schemes such as Newmark-$\beta$ method. Computationally, it translates to a $\sim$$1000\times$ and $\sim$$10-100\times$ speed-up over conventional time-marching schemes for linear and nonlinear elastodynamics, respectively.
In this paper, we present a randomized extension of the deep splitting algorithm introduced in [Beck, Becker, Cheridito, Jentzen, and Neufeld (2021)] using random neural networks suitable to approximately solve both high-dimensional nonlinear parabolic PDEs and PIDEs with jumps having (possibly) infinite activity. We provide a full error analysis of our so-called random deep splitting method. In particular, we prove that our random deep splitting method converges to the (unique viscosity) solution of the nonlinear PDE or PIDE under consideration. Moreover, we empirically analyze our random deep splitting method by considering several numerical examples including both nonlinear PDEs and nonlinear PIDEs relevant in the context of pricing of financial derivatives under default risk. In particular, we empirically demonstrate in all examples that our random deep splitting method can approximately solve nonlinear PDEs and PIDEs in 10'000 dimensions within seconds.
We propose a novel score-based particle method for solving the Landau equation in plasmas, that seamlessly integrates learning with structure-preserving particle methods [arXiv:1910.03080]. Building upon the Lagrangian viewpoint of the Landau equation, a central challenge stems from the nonlinear dependence of the velocity field on the density. Our primary innovation lies in recognizing that this nonlinearity is in the form of the score function, which can be approximated dynamically via techniques from score-matching. The resulting method inherits the conservation properties of the deterministic particle method while sidestepping the necessity for kernel density estimation in [arXiv:1910.03080]. This streamlines computation and enhances scalability with dimensionality. Furthermore, we provide a theoretical estimate by demonstrating that the KL divergence between our approximation and the true solution can be effectively controlled by the score-matching loss. Additionally, by adopting the flow map viewpoint, we derive an update formula for exact density computation. Extensive examples have been provided to show the efficiency of the method, including a physically relevant case of Coulomb interaction.
In this work, we analyze the convergence rate of randomized quasi-Monte Carlo (RQMC) methods under Owen's boundary growth condition [Owen, 2006] via spectral analysis. Specifically, we examine the RQMC estimator variance for the two commonly studied sequences: the lattice rule and the Sobol' sequence, applying the Fourier transform and Walsh--Fourier transform, respectively, for this analysis. Assuming certain regularity conditions, our findings reveal that the asymptotic convergence rate of the RQMC estimator's variance closely aligns with the exponent specified in Owen's boundary growth condition for both sequence types. We also provide guidance on choosing the importance sampling density to minimize RQMC estimator variance.
This paper is concerned with inverse crack scattering problems for time-harmonic acoustic waves. We prove that a piecewise linear crack with the sound-soft boundary condition in two dimensions can be uniquely determined by the far-field data corresponding to a single incident plane wave or point source. We propose two non-iterative methods for imaging the location and shape of a crack. The first one is a contrast sampling method, while the second one is a variant of the classical factorization method but only with one incoming wave. Newton's iteration method is then employed for getting a more precise reconstruction result. Numerical examples are presented to show the effectiveness of the proposed hybrid method.
In this paper, the SQP method applied to a hyperbolic PDE-constrained optimization problem is considered. The model arises from the acoustic full waveform inversion in the time domain. The analysis is mainly challenging due to the involved hyperbolicity and second-order bilinear structure. This notorious character leads to an undesired effect of loss of regularity in the SQP method, calling for a substantial extension of developed parabolic techniques. We propose and analyze a novel strategy for the well-posedness and convergence analysis based on the use of a smooth-in-time initial condition, a tailored self-mapping operator, and a two-step estimation process along with Stampacchia's method for second-order wave equations. Our final theoretical result is the R-superlinear convergence of the SQP method.
We study topology optimization governed by the incompressible Navier-Stokes flows using a phase field model. Novel stabilized semi-implicit schemes for the gradient flows of Allen-Cahn and Cahn-Hilliard types are proposed for solving the resulting optimal control problem. Unconditional energy stability is shown for the gradient flow schemes in continuous and discrete spaces. Numerical experiments of computational fluid dynamics in 2d and 3d show the effectiveness and robustness of the optimization algorithms proposed.
This paper presents an analysis of properties of two hybrid discretization methods for Gaussian derivatives, based on convolutions with either the normalized sampled Gaussian kernel or the integrated Gaussian kernel followed by central differences. The motivation for studying these discretization methods is that in situations when multiple spatial derivatives of different order are needed at the same scale level, they can be computed significantly more efficiently compared to more direct derivative approximations based on explicit convolutions with either sampled Gaussian kernels or integrated Gaussian kernels. While these computational benefits do also hold for the genuinely discrete approach for computing discrete analogues of Gaussian derivatives, based on convolution with the discrete analogue of the Gaussian kernel followed by central differences, the underlying mathematical primitives for the discrete analogue of the Gaussian kernel, in terms of modified Bessel functions of integer order, may not be available in certain frameworks for image processing, such as when performing deep learning based on scale-parameterized filters in terms of Gaussian derivatives, with learning of the scale levels. In this paper, we present a characterization of the properties of these hybrid discretization methods, in terms of quantitative performance measures concerning the amount of spatial smoothing that they imply, as well as the relative consistency of scale estimates obtained from scale-invariant feature detectors with automatic scale selection, with an emphasis on the behaviour for very small values of the scale parameter, which may differ significantly from corresponding results obtained from the fully continuous scale-space theory, as well as between different types of discretization methods.
Adaptive finite elements combined with geometric multigrid solvers are one of the most efficient numerical methods for problems such as the instationary Navier-Stokes equations. Yet despite their efficiency, computations remain expensive and the simulation of, for example, complex flow problems can take many hours or days. GPUs provide an interesting avenue to speed up the calculations due to their very large theoretical peak performance. However, the large degree of parallelism and non-standard API make the use of GPUs in scientific computing challenging. In this work, we develop a GPU acceleration for the adaptive finite element library Gascoigne and study its effectiveness for different systems of partial differential equations. Through the systematic formulation of all computations as linear algebra operations, we can employ GPU-accelerated linear algebra libraries, which simplifies the implementation and ensures the maintainability of the code while achieving very efficient GPU utilizations. Our results for a transport-diffusion equation, linear elasticity, and the instationary Navier-Stokes equations show substantial speedups of up to 20X compared to multi-core CPU implementations.
This paper deals with a numerical analysis of plastic deformation under various conditions, utilizing Radial Basis Function (RBF) approximation. The focus is on the elasto-plastic von Mises problem under plane-strain assumption. Elastic deformation is modelled using the Navier-Cauchy equation. In regions where the von Mises stress surpasses the yield stress, corrections are applied locally through a return mapping algorithm. The non-linear deformation problem in the plastic domain is solved using the Picard iteration. The solutions for the Navier-Cauchy equation are computed using the Radial Basis Function-Generated Finite Differences (RBF-FD) meshless method using only scattered nodes in a strong form. Verification of the method is performed through the analysis of an internally pressurized thick-walled cylinder subjected to varying loading conditions. These conditions induce states of elastic expansion, perfectly-plastic yielding, and plastic yielding with linear hardening. The results are benchmarked against analytical solutions and traditional Finite Element Method (FEM) solutions. The paper also showcases the robustness of this approach by solving case of thick-walled cylinder with cut-outs. The results affirm that the RBF-FD method produces results comparable to those obtained through FEM, while offering substantial benefits in managing complex geometries without the necessity for conventional meshing, along with other benefits of meshless methods.
In this article, a posteriori error analysis of the elliptic obstacle problem is addressed using hybrid high-order methods. The method involve cell unknowns represented by degree-$r$ polynomials and face unknowns represented by degree-$s$ polynomials, where $r=0$ and $s$ is either $0$ or $1$. The discrete obstacle constraints are specifically applied to the cell unknowns. The analysis hinges on the construction of a suitable Lagrange multiplier, a residual functional and a linear averaging map. The reliability and the efficiency of the proposed a posteriori error estimator is discussed, and the study is concluded by numerical experiments supporting the theoretical results.
It is often the case that, while the numerical solution of the non-linear dispersive equation $\mathrm{i}\partial_t u(t)=\mathcal{H}(u(t),t)u(t)$ represents a formidable challenge, it is fairly easy and cheap to solve closely related linear equations of the form $\mathrm{i}\partial_t u(t)=\mathcal{H}_1(t)u(t)+\widetilde{\mathcal H}_2(t)u(t)$, where $\mathcal{H}_1(t)+\mathcal{H}_2(v,t)=\mathcal{H}(v,t)$. In that case we advocate an iterative linearisation procedure that involves fixed-point iteration of the latter equation to solve the former. A typical case is when the original problem is a nonlinear Schr\"odinger or Gross--Pitaevskii equation, while the `easy' equation is linear Schr\"odinger with time-dependent potential. We analyse in detail the iterative scheme and its practical implementation, prove that each iteration increases the order, derive upper bounds on the speed of convergence and discuss in the case of nonlinear Schr\"odinger equation with cubic potential the preservation of structural features of the underlying equation: the $\mathrm{L}_2$ norm, momentum and Hamiltonian energy. A key ingredient in our approach is the use of the Magnus expansion in conjunction with Hermite quadratures, which allows effective solutions of the linearised but non-autonomous equations in an iterative fashion. The resulting Magnus--Hermite methods can be combined with a wide range of numerical approximations to the matrix exponential. The paper concludes with a number of numerical experiments, demonstrating the power of the proposed approach.
Randomized algorithms have proven to perform well on a large class of numerical linear algebra problems. Their theoretical analysis is critical to provide guarantees on their behaviour, and in this sense, the stochastic analysis of the randomized low-rank approximation error plays a central role. Indeed, several randomized methods for the approximation of dominant eigen- or singular modes can be rewritten as low-rank approximation methods. However, despite the large variety of algorithms, the existing theoretical frameworks for their analysis rely on a specific structure for the covariance matrix that is not adapted to all the algorithms. We propose a general framework for the stochastic analysis of the low-rank approximation error in Frobenius norm for centered and non-standard Gaussian matrices. Under minimal assumptions on the covariance matrix, we derive accurate bounds both in expectation and probability. Our bounds have clear interpretations that enable us to derive properties and motivate practical choices for the covariance matrix resulting in efficient low-rank approximation algorithms. The most commonly used bounds in the literature have been demonstrated as a specific instance of the bounds proposed here, with the additional contribution of being tighter. Numerical experiments related to data assimilation further illustrate that exploiting the problem structure to select the covariance matrix improves the performance as suggested by our bounds.
We show that even though the Discontinuous Galerkin Spectral Element Method is stable for hyperbolic boundary-value problems, and the overset domain problem is well-posed in an appropriate norm, the energy of the approximation is bounded by data only for fixed polynomial order and time. In the absence of dissipation, coupling of the overlapping domains is destabilizing by allowing positive eigenvalues in the system to be integrated in time. This coupling can be stabilized in one space dimension by using the upwind numerical flux. To help provide additional dissipation, we introduce a novel penalty method that applies dissipation at arbitrary points within the overlap region and depends only on the difference between the solutions. We present numerical experiments in one space dimension to illustrate the implementation of the well-posed penalty formulation, and show spectral convergence of the approximations when dissipation is applied.
We consider a stable unique continuation problem for the wave equation where the initial data is lacking and the solution is reconstructed using measurements in some subset of the bulk domain. Typically fairly sophisticated space-time methods have been used in previous work to obtain stable and accurate solutions to this reconstruction problem. Here we propose to solve the problem using a standard discontinuous Galerkin method for the temporal discretization and continuous finite elements for the space discretization. Error estimates are established under a geometric control condition. We also investigate two preconditioning strategies which can be used to solve the arising globally coupled space-time system by means of simple time-stepping procedures. Our numerical experiments test the performance of these strategies and highlight the importance of the geometric control condition for reconstructing the solution beyond the data domain.
We reexamine Smale's alpha theory as a way to certify a numerical solution to an analytic system. For a given point and a system, Smale's alpha theory determines whether Newton's method applied to this point shows the quadratic convergence to an exact solution. We introduce the alpha theory computation using interval arithmetic to avoid costly exact arithmetic. As a straightforward variation of the alpha theory, our work improves computational efficiency compared to software employing the traditional alpha theory.
Algorithms operating on real numbers are implemented as floating-point computations in practice, but floating-point operations introduce roundoff errors that can degrade the accuracy of the result. We propose $\Lambda_{num}$, a functional programming language with a type system that can express quantitative bounds on roundoff error. Our type system combines a sensitivity analysis, enforced through a linear typing discipline, with a novel graded monad to track the accumulation of roundoff errors. We prove that our type system is sound by relating the denotational semantics of our language to the exact and floating-point operational semantics. To demonstrate our system, we instantiate $\Lambda_{num}$ with error metrics proposed in the numerical analysis literature and we show how to incorporate rounding operations that faithfully model aspects of the IEEE 754 floating-point standard. To show that $\Lambda_{num}$ can be a useful tool for automated error analysis, we develop a prototype implementation for $\Lambda_{num}$ that infers error bounds that are competitive with existing tools, while often running significantly faster. Finally, we consider semantic extensions of our graded monad to bound error under more complex rounding behaviors, such as non-deterministic and randomized rounding.