Discrete approximation, which has been the prevailing scheme in stochastic programming in the past decade, has been extended to distributionally robust optimization (DRO) recently. In this paper, we conduct rigorous quantitative stability analysis of discrete approximation schemes for DRO, which measures the approximation error in terms of discretization sample size. For the ambiguity set defined through equality and inequality moment conditions, we quantify the discrepancy between the discretized ambiguity sets and the original set with respect to the Wasserstein metric. To establish the quantitative convergence, we develop a Hoffman error bound theory with Hoffman constant calculation criteria in a infinite dimensional space, which can be regarded as a byproduct of independent interest. For the ambiguity set defined by Wasserstein ball and moment conditions combined with Wasserstein ball, we present similar quantitative stability analysis by taking full advantage of the convex property inherently admitted by Wasserstein metric. Efficient numerical methods for specifically solving discrete approximation DRO problems with thousands of samples are also designed. In particular, we reformulate different types of discrete approximation problems into a class of saddle point problems with completely separable structures. The stochastic primal-dual hybrid gradient (PDHG) algorithm where in each iteration we update a random subset of the sampled variables is then amenable as a solution method for the reformulated saddle point problems. Some preliminary numerical tests are reported.

With the path integral approach, the thermal average in multi-electronic-state quantum systems can be approximated by the ring polymer representation on an extended configuration space, where the additional degrees of freedom are associated with the surface index of each bead. The primary goal of this work is to propose a more efficient sampling algorithm for the calculation of such thermal averages. We reformulate the extended ring polymer approximation according to the configurations of the surface indexes, and by introducing a proper reference measure, the reformulation is recast as a ratio of two expectations of function expansions. By quantitatively estimating the sub-estimators, and minimizing the total variance of the sampled average, we propose a multi-level Monte Carlo path integral molecular dynamics method (MLMC-PIMD) to achieve an optimal balance of computational cost and accuracy.

This paper is concerned with numerical solutions of time-fractional parabolic equations. Due to the Caputo time derivative being involved, the solutions of equations are usually singular near the initial time $t = 0$ even for a smooth setting. Based on a simple change of variable $s = t^β$, an equivalent $s$-fractional differential equation is derived and analyzed. Two types of finite difference methods based on linear and quadratic approximations in the $s$-direction are presented, respectively, for solving the $s$-fractional differential equation. We show that the method based on the linear approximation provides the optimal accuracy $\mathcal{O}(N ^{−(2−α)})$ where $N$ is the number of grid points in temporal direction. Numerical examples for both linear and nonlinear fractional equations are presented in comparison with $L1$ methods on uniform meshes and graded meshes, respectively. Our numerical results show clearly the accuracy and efficiency of the proposed methods.

Solving partial differential equations in high dimensions by deep neural networks has brought significant attentions in recent years. In many scenarios, the loss function is defined as an integral over a high-dimensional domain. Monte-Carlo method, together with a deep neural network, is used to overcome the curse of dimensionality, while classical methods fail. Often, a neural network outperforms classical numerical methods in terms of both accuracy and efficiency. In this paper, we propose to use quasi-Monte Carlo sampling, instead of Monte-Carlo method to approximate the loss function. To demonstrate the idea, we conduct numerical experiments in the framework of deep Ritz method. For the same accuracy requirement, it is observed that quasi-Monte Carlo sampling reduces the size of training data set by more than two orders of magnitude compared to that of Monte-Carlo method. Under some assumptions, we can prove that quasi-Monte Carlo sampling together with the deep neural network generates a convergent series with rate proportional to the approximation accuracy of quasi-Monte Carlo method for numerical integration. Numerically the fitted convergence rate is a bit smaller, but the proposed approach always outperforms Monte Carlo method.

In this paper, we suggest a new vertex interpolation algorithm to improve an existing cell-centered finite volume scheme for nonlinear diffusion problems on general meshes. The new vertex interpolation algorithm is derived by applying a special limit procedure to the well-known MPFA-O method. Since the MPFA-O method for 3D cases has been addressed in some studies, the new vertex interpolation algorithm can be extended to 3D cases naturally. More interesting is that the solvability of the corresponding local system is proved under some assumptions. Additionally, we modify the edge flux approximation by an edge-based discretization of diffusion coefficient, and thus the improved scheme is free of the so-called numerical heat-barrier issue suffered by many existing cell-centered or hybrid schemes. The final scheme allows arbitrary continuous or discontinuous diffusion coefficients and can be applicable to arbitrary star-shaped polygonal meshes. A second-order convergence rate for the approximate solution and a first-order accuracy for the flux are observed in numerical experiments. In the comparative experiments with some existing vertex interpolation algorithms, the new algorithm shows obvious improvement on highly distorted meshes.

In this paper, we analyze the convergence of the Peaceman-Rachford splitting method (PRSM) for a type of nonconvex and nonsmooth optimization with linear constraints, whose objective function is the sum of a proper lower semicontinuous function and a strongly convex differential function. When a suitable penalty parameter is chosen and the iterative point sequence is bounded, we show the global convergence of the PRSM. Furthermore, under the assumption that the associated function satisfies the Kurdyka-Łojasiewicz property, we prove the strong convergence of the PRSM. We also provide sufficient conditions guaranteeing the boundedness of the generated sequence. Finally, we present some preliminary numerical results to show the effectiveness of the PRSM and also give a comparison with the Douglas-Rachford splitting method.

Efficient and accurate Legendre spectral element methods for solving one-dimensional higher order differential equations with high oscillatory or steep gradient solutions are proposed. Some Sobolev orthogonal/biorthogonal basis functions corresponding to each subinterval are constructed, which reduce the non-zero entries of linear systems and computational cost. Numerical experiments exhibit the effectiveness and accuracy of the suggested approaches. 

In this article, a new compact difference scheme is proposed in exponential form to solve two-dimensional unsteady nonlinear Burgers' and Navier-Stokes equations of motion in polar cylindrical coordinates by using half-step discretization. At each time level by using only nine grid points in space, the proposed scheme gives accuracy of order four in space and two in time. The method is directly applicable to the equations having singularities at boundary points. Stability analysis is explained in detail and many benchmark problems like Burgers', Navier-Stokes and Taylor-vortex problems in polar cylindrical coordinates are solved to verify the accuracy and efficiency of the scheme.

In this paper, we propose and numerically investigate a superconvergent cluster recovery (SCR) method for the Crouzeix-Raviart (CR) element. The proposed recovery method reconstructs a $C^0$ linear gradient. A linear polynomial approximation is obtained by a least square fitting to the CR element approximation at certain sample points, and then taken derivatives to obtain the recovered gradient. The SCR recovery operator is superconvergent on uniform mesh of four patterns. Numerical examples show that SCR can produce a superconvergent gradient approximation for the CR element, and provide an asymptotically exact error estimator in the adaptive CR finite element method.

Conjugate gradient methods are interesting iterative methods that solve large scale unconstrained optimization problems. A lot of recent research has thus focussed on developing a number of conjugate gradient methods that are more effective. In this paper, we propose another hybrid conjugate gradient method as a linear combination of Dai-Yuan (DY) method and the Hestenes-Stiefel (HS) method. The sufficient descent condition and the global convergence of this method are established using the generalized Wolfe line search conditions. Compared to the other conjugate gradient methods, the proposed method gives good numerical results and is effective.

In this paper, we propose a mathematical model and present numerical simulations for ice melting phenomena. The model is based on the phase-field modeling for the crystal growth. To model ice melting, we ignore anisotropy in the crystal growth model and introduce a new melting term. The numerical solution algorithm is a hybrid method which uses both the analytic and numerical solutions. We perform various computational experiments. The computational results confirm the accuracy and efficiency of the proposed method for ice melting.