Mathematics

We propose a unified meshless method to solve classical and fractional PDE problems with (?��?) α 2 for α??0,2] . The classical ( α=2 ) and fractional ( α<2 ) Laplacians, one local and the other nonlocal, have distinct properties. Therefore, their numerical methods and computer implementations are usually incompatible. We notice that for any α?? , the Laplacian (?��?) α 2 of generalized inverse multiquadric (GIMQ) functions can be analytically written by the Gauss hypergeometric function, and thus propose a GIMQ-based method. Our method unifies the discretization of classical and fractional Laplacians and also bypasses numerical approximation to the hypersingular integral of fractional Laplacian. These two merits distinguish our method from other existing methods for the fractional Laplacian. Extensive numerical experiments are carried out to test the performance of our method. Compared to other methods, our method can achieve high accuracy with fewer number of unknowns, which effectively reduces the storage and computational requirements in simulations of fractional PDEs. Moreover, the meshfree nature makes it free of geometric constraints and enables simple implementation for any dimension d?? . Additionally, two approaches of selecting shape parameters, including condition number-indicated method and random-perturbed method, are studied to avoid the ill-conditioning issues when large number of points.

By using the Onsager principle as an approximation tool, we give a novel derivation for the moving finite element method for gradient flow equations. We show that the discretized problem has the same energy dissipation structure as the continuous one. This enables us to do numerical analysis for the stationary solution of a nonlinear reaction diffusion equation using the approximation theory of free-knot piecewise polynomials. We show that under certain conditions the solution obtained by the moving finite element method converges to a local minimizer of the total energy when time goes to infinity. The global minimizer, once it is detected by the discrete scheme, approximates the continuous stationary solution in optimal order. Numerical examples for a linear diffusion equation and a nonlinear Allen-Cahn equation are given to verify the analytical results.

The variational discrete element method developed in [28] for dynamic elasto-plastic computations is adapted to compute the deformation of elastic Cosserat materials. In addition to cellwise displacement degrees of freedom (dofs), cellwise rotational dofs are added. A reconstruction is devised to obtain P 1 non-conforming polynomials in each cell and thus constant strains and stresses in each cell. The method requires only the usual macroscopic parameters of a Cosserat material and no microscopic parameter. Numerical examples show the robustness of the method for both static and dynamic computations in two and three dimensions.

In this work, an r -linearly converging adaptive solver is constructed for parabolic evolution equations in a simultaneous space-time variational formulation. Exploiting the product structure of the space-time cylinder, the family of trial spaces that we consider are given as the spans of wavelets-in-time and (locally refined) finite element spaces-in-space. Numerical results illustrate our theoretical findings.

In this paper, we report the results obtained from the acceleration of multi-binary64-type multiple precision matrix multiplication with AVX2. We target double-double (DD), triple-double (TD), and quad-double (QD) precision arithmetic designed by certain types of error-free transformation (EFT) arithmetic. Furthermore, we implement SIMDized EFT functions, which simultaneously compute with four binary64 numbers on x86_64 computing environment, and by using help of them, we also develop SIMDized DD, TD, and QD additions and multiplications. In addition, AVX2 load/store functions were adopted to efficiently speed up reading and storing matrix elements from/to memory. Owing to these combined techniques, our implemented multiple precision matrix multiplications have been accelerated more than three times compared with non-accelerated ones. Our accelerated matrix multiplication modifies the performance of parallelization with OpenMP.

In this paper we consider utilizing a residual neural network (ResNet) to solve ordinary differential equations. Stochastic gradient descent method is applied to obtain the optimal parameter set of weights and biases of the network. We apply forward Euler, Runge-Kutta2 and Runge-Kutta4 finite difference methods to generate three sets of targets training the ResNet and carry out the target study. The well trained ResNet behaves just as its counterpart of the corresponding one-step finite difference method. In particular, we carry out (1) the architecture study in terms of number of hidden layers and neurons per layer to find the optimal ResNet structure; (2) the target study to verify the ResNet solver behaves as accurate as its finite difference method counterpart; (3) solution trajectory simulation. Even the ResNet solver looks like and is implemented in a way similar to forward Euler scheme, its accuracy can be as high as any one step method. A sequence of numerical examples are presented to demonstrate the performance of the ResNet solver.

This article considers Hamiltonian mechanical systems with potential functions admitting jump discontinuities. The focus is on accurate and efficient numerical approximations of their solutions, which will be defined via the laws of reflection and refraction. Despite of the success of symplectic integrators for smooth mechanical systems, their construction for the discontinuous ones is nontrivial, and numerical convergence order can be impaired too. Several rather-usable numerical methods are proposed, including: a first-order symplectic integrator for general problems, a third-order symplectic integrator for problems with only one linear interface, arbitrarily high-order reversible integrators for general problems (no longer symplectic), and an adaptive time-stepping version of the previous high-order method. Interestingly, whether symplecticity leads to favorable long time performance is no longer clear due to discontinuity, as traditional Hamiltonian backward error analysis does not apply any more. Therefore, at this stage, our recommended default method is the last one. Various numerical evidence, on the order of convergence, long time performance, momentum map conservation, and consistency with the computationally-expensive penalty method, are supplied. A complex problem, namely the Sauteed Mushroom, is also proposed and numerically investigated, for which multiple bifurcations between trapped and ergodic dynamics are observed.

In practical computations, the (preconditioned) conjugate gradient (P)CG method is the iterative method of choice for solving systems of linear algebraic equations Ax=b with a real symmetric positive definite matrix A . During the iterations it is important to monitor the quality of the approximate solution x k so that the process could be stopped whenever x k is accurate enough. One of the most relevant quantities for monitoring the quality of x k is the squared A -norm of the error vector x??x k . This quantity cannot be easily evaluated, however, it can be estimated. Many of the existing estimation techniques are inspired by the view of CG as a procedure for approximating a certain Riemann--Stieltjes integral. The most natural technique is based on the Gauss quadrature approximation and provides a lower bound on the quantity of interest. The bound can be cheaply evaluated using terms that have to be computed anyway in the forthcoming CG iterations. If the squared A -norm of the error vector decreases rapidly, then the lower bound represents a tight estimate. In this paper we suggest a heuristic strategy aiming to answer the question of how many forthcoming CG iterations are needed to get an estimate with the prescribed accuracy. Numerical experiments demonstrate that the suggested strategy is efficient and robust.

Mathematical modelling of ionic electrodiffusion and water movement is emerging as a powerful avenue of investigation to provide new physiological insight into brain homeostasis. However, in order to provide solid answers and resolve controversies, the accuracy of the predictions is essential. Ionic electrodiffusion models typically comprise non-trivial systems of non-linear and highly coupled partial and ordinary differential equations that govern phenomena on disparate time scales. Here, we study numerical challenges related to approximating these systems. We consider a homogenized model for electrodiffusion and osmosis in brain tissue and present and evaluate different associated finite element-based splitting schemes in terms of their numerical properties, including accuracy, convergence, and computational efficiency for both idealized scenarios and for the physiologically relevant setting of cortical spreading depression (CSD). We find that the schemes display optimal convergence rates in space for problems with smooth manufactured solutions. However, the physiological CSD setting is challenging: we find that the accurate computation of CSD wave characteristics (wave speed and wave width) requires a very fine spatial and fine temporal resolution.

We study time-harmonic scattering in R n ( n=2,3 ) by a planar screen (a "crack" in the context of linear elasticity), assumed to be a non-empty bounded relatively open subset Γ of the hyperplane R n−1 ×{0} , on which impedance (Robin) boundary conditions are imposed. In contrast to previous studies, Γ can have arbitrarily rough (possibly fractal) boundary. To obtain well-posedness for such Γ we show how the standard impedance boundary value problem and its associated system of boundary integral equations must be supplemented with additional solution regularity conditions, which hold automatically when ∂Γ is smooth. We show that the associated system of boundary integral operators is compactly perturbed coercive in an appropriate function space setting, strengthening previous results. This permits the use of Mosco convergence to prove convergence of boundary element approximations on smoother "prefractal" screens to the limiting solution on a fractal screen. We present accompanying numerical results, validating our theoretical convergence results, for three-dimensional scattering by a Koch snowflake and a square snowflake.