In the paper, for solving large sparse vertical linear complementarity problems, a twostep modulus-based matrix synchronous multisplitting parallel iteration method is constructed by using two-step multisplitting technique. The new method can be viewed as a generalization of the modulus-based matrix synchronous multisplitting iterative method and the two-step modulus-based matrix splitting iterative method in the existing literatures. Furthermore, under the assumption that the system matrices are H₊-matrices, the convergence analysis of the proposed method is given, and the convergence domain of the parameter matrix is obtained, which extends the convergence results of the existing methods. Finally, numerical experiments are carried out under the OpenMP framework for two numerical examples in the existing literatures. Numerical results show that the two-step multisplitting technique can improve the computational efficiency of the existing methods.

Numerical simulation of atmospheric model has been widely applied to weather prediction. Supercomputers play an important role in improving the accuracy and resolution of weather forecasting. Implicit solver of atmospheric model counts on many computation kernels. Their parallelization and optimization are critical to solver performance. Based on uniform recurrence relation (URR), this paper introduces a general framework to model the behavior of sparse triangular solve (SpTRSV) and incomplete matrix factorization (ILU). For SpTRSV and ILU, the extended framework leads to an speed up of 26×and 33×. This paper also parallelizes kernel functions in numerical simulation of atmospheric model, with acceleration ratio of 10-152×over serial implementations. The simulation achieves a sustained aggregated performance of 20.5PFlop/s in double precision, with strong- and weak-scaling efficiency above 56.81% and 41.87%, respectively.

Smoothed Particle Hydrodynamics (SPH) method is a meshless Lagrangian particle method, which is widely used in the field of fluid mechanics and the simulation of large deformation and impact load. Many scholars have done a lot of research on SPH algorithm to improve the calculation speed and accuracy of SPH algorithm. To solve the problem that the accuracy of particle approximation near the boundary of the existing SPH method decreases, this paper proposes an improved kernel approximation based on CSPH method and MSPH method. In the process of solving the field function, the first derivative approximation and the second derivative approximation, the equation with the second derivative term is optimized, which reduces the number of solutions of the second derivative approximation, and reduces the amount of calculation compared with MSPH method. In addition, based on the improved SPH algorithm, this paper established a two-dimensional numerical wave flume to simulate push-plate wave-making, and compared the wave parameters generated by SPH algorithm with the theoretical values through numerical simulation wave-making, which verified that the improved SPH method had a good simulation effect on wave generation and propagation, and provided an algorithm research foundation for the follow-up study of internal wave, freak wave and nonlinear wave interaction.

Ensemble Kalman filter based inversion algorithms have become increasingly popular as a method for solving Bayesian inverse problems. The method can be viewed as a derivative-free optimization approach for finding a solution of the inverse problem. Since the methods are based on the ensemble Kalman filter, it may not perform well for problems with strongly non-Gaussian observation models.To address the issue, in this work we provide an affine-mapping based variational Ensemble Kalman inversion algorithm, which identifies an affine-mapping from the prior ensemble to the posterior one. With numerical examples, we show that the proposed method has competitive performance against the standard ensemble Kalman inversion for problems with non-Gaussian observation models.

Bayesian inference becomes an increasingly popular tool for solving inverse problems largely thanks to its ability to characterize the uncertainty in the solution obtained. In many practical Bayesian inverse problems, multiple competing models may be available for describing the data and/or the parameter of interest. In this case, the Bayesian model selection approach is often used to choose the "best" model and this method relies on the knowledge of the normalizing constant in the posterior distribution. As such, estimating the normalizing constant becomes an important task in Bayesian inference, and this problem is computationally challenging for standard sampling methods. In this work we present a method to estimate the normalizing constant based on the Multicanonical Monte Carlo technique, an adaptive importance sampling scheme. The method can estimate the normalizing constant in a black-box manner, making it particularly suitable for problems with complex underlying models. With numerical examples, we demonstrate that the proposed method can efficiently and accurately compute an estimate of the normalizing constant.

Compressed Sensing (CS) is an effective method for fast Magnetic Resonance Imaging (MRI), which aims to reconstruct MRI from a small amount of under-sampled data space and accelerate data acquisition in MRI. In order to improve the reconstruction accuracy and speed of current MRI, in this paper, we combine the traditional model-driven compressed sensing method and data-driven deep learning method to propose a new algorithm–ADMMCNN. This proposed algorithm solves a subproblem of the alternating direction multiplier method (ADMM) by convolutional neural network, and then optimizes the penalty parameter and the update rate of Lagrange multiplizer by constructing the loss function and using the gradient descent method to obtain the trained ADMM model, which can effectively recover the MRI signal from the compressed sensing measurement. Compared with the traditional reconstruction algorithms, the proposed algorithm has fewer regularization parameters and lower computational complexity, at the same time, we introduce a trained convolutional neural network directly to greatly simplify the training process and computation, also we use GPU to accelerate the training. Finally, we give some experimental results to show that ADMM-CNN has faster calculation speed and better reconstruction accuracy compared with traditional reconstruction algorithms and other deep learning methods.

The block Kaczmarz method is an iterative algorithm for solving large-scale linear systems. At each iteration, the current iteration points will be orthogonally projected to the solution space of the constraint subset. In the paper, according to the principle of the block Kaczmarz method, we propose the randomized block Kaczmarz method with the Gaussian mixture model (GMM-RBK(k)). The partition of blocks is carried out by the Gaussian mixture model. In order to avoid the computation of pseudo-inverse or the solution of the least-squares problem, we propose the randomized average block Kaczmarz method with the Gaussian mixture model (GMM-RABK(k)). We prove that the two algorithms converge when the linear system is consistent, and give the corresponding formula of convergence rate. The effectiveness of the GMM-RBK(k) method and GMM-RABK(k) method is also verified in the final numerical experiments, and avoiding the calculation of pseudo-inverse, the GMM-RABK(k) method is superior to the GMM-RBK(k) method.

Using the method of separation of variables is given including integrated solution of one-dimensional heat conduction equation of point source, on this basis, study the one-dimensional heat conduction equation of point source strength inversion problem, the source strength inversion problem is transformed into optimization problem, the variational adjoint method are used to get the gradient expression, with the aid of Broyden class algorithm for the inversion, and the numerical simulation is given, The results show that the Broyden algorithm is feasible and effective.

Two Gauss-Seidel type methods with oblique direction for solving large-scale overdetermined inconsistent linear systems, namely, the greedy random Gauss-Seidel method with oblique direction (GRGSO) and the fast max-distance Gauss-Seidel method with oblique direction (FMDGSO), are proposed with the use of the strategy of selecting the working sequence of the coefficient matrix based on the greedy criterion and the max-distance criterion respectively. When the coefficient matrix is of full column rank, the theoretical results show that the two methods converge to the unique least square solution of the linear system. Particularly, when the columns of matrix A are close to linearly correlated, the numerical results show that the two methods have more advantages than the existing Gauss-Seidel method in solving performance.

Factor analysis is an important statistics tool in many practical fields. The estimation of the number of factors is one of the most key problems in factor analysis. In practice, the determination of the number of factors usually depends on the covariance estimation from observed data. In this paper, we use Kullback-Leibler divergence to characterize the uncertainty of the covariance estimation, proposed a minimum trace model based on KL divergence to estimate the number of factors, and develop a symmetric Gauss-Seidel decomposition based alternating direction method of multiplier for solving this problem. We also show that the proposed solution method is globally convergent. The numerical results are presented to demonstrate the effectiveness of the proposed model and the efficiency of the numerical algorithm.

Advection-diffusion equation (ADE) with fractional Laplacian has been developed to characterize some anomalous diffusion phenomena in nature and social systems. In this study, a new lattice Boltzmann model with BGK operator (LBGK) for solving 2D ADE with fractional Laplacian is proposed. First, an approximated formula of the governing equation can be obtained based on the Fourier transform of fractional Laplacian and the Gauss-Type quadrature formular. Then, by means of velocity discretization, time and space discretization, constructing appropriate equilibrium distribution function and force distribution functionan, an efficient LBGK model is developed for the ADE with fractional Laplacian. Through the Chapman-Enskog analysis, the macroscopic equation can be recovered from the present LBGK model, which proves the effectiveness of the model. Finally, the model is applied to a numerical example with analysis solution and the Allen-Cahn equation. The results further verify the correctness and effectiveness of the present model.

The waveform relaxation (WR) method is an important numerical method to solve differential equations. So far, the investigations of the behaviour of WR methods focus on the convergence and few studies of stability are published. In this paper we present the definition of WR methods for ordinary differential equations (ODEs). By applying the standard techniques used to study stable properties of classical numerical methods of ODEs, we investigate the stability properties of WR methods and give the sufficient conditions under which the continuous WR methods can preserve the stable properties of three standard test equations. The contraction properties of WR methods are investigated by Lyapunov techniques and some sufficient conditions of preserving the contraction properties of the test equations are obtained for continuous and discrete time WR methods.

In recent years, a series of important achievements have been made in experimental research on ground state solutions of Bose-Einstein condensates. In this paper, on the basis of relevant research results, Firstly, the ground state solution problem of Bose-Einstein condensate is transformed into an energy functional extremum problem by means of reductive and dimensionless methods, In discretizing this equation, the finite difference method and Fourier spectral method are used to discretize the one and two-dimensional cases of the energy functional equation. The second, In this paper, SQP optimization method is an efficient numerical optimization algorithm used to solve the ground state solution of Bose-Einstein condensate and the algorithm is used to simulate one and two-dimensional energy functional minimum problems. Finally, by analyzing the experimental data and images, it is concluded that the algorithm can improve the accuracy of energy function value.

An extrapolation full multigrid (EXFMG) method is proposed to solve two-dimensional elliptic interface problems. In this approach, a quite good initial guess is constructed by using Richardson extrapolation and spline interpolation, which accelerates V-cycle MG method for calculating the discretize system and greatly reduces the number of V-cycles required. Numerical results are given to show that our proposed EXFMG algorithm can keep less cost.

Three explicit Runge-Kutta schemes are constructed from the viewpoint of Newton-Cotes integration, including a 3-stage $2^{nd}$ order Runge-Kutta scheme with parameters and two 4-stage $3^{rd}$ order Runge-Kutta schemes with parameters. The commonly used RK3 and RK4 can be obtained by taking special parameters from our schemes. The accuracy and the stability of these schemes are presented. Numerical examples are given to verify the stability, effectiveness and high accuracy of the constructed schemes. Results show that our explicit Runge-Kutta schemes have better stability than that of the classic explicit Runge-Kutta schemes.

Applying the conjugate gradient method and linear projection operator, an iterative algorithm is presented to solve the least squares solution of linear matrix equation $AXB=C$ under any linear subspace. It is proved that the least squares solution, the minimum-norm least squares solution and the optimal approximation of the matrix equation $AXB=C$ can be obtained in finite iteration steps by the method without considering rounding errors. The numerical examples verify the efficiency of the algorithm. The merit of our method is that it is easy to implement in any linear subspace.

The generalized Dantzig selector problem is an effective approach to solve the parameter estimation, where any norm can be used to estimate. This paper adopts a dual alternating direction method of multipliers (dADMM for short) to solve the $\ell_1, \ell_2$ and $\ell_ {\infty}$ generalized Dantzig selector. The global convergence and local linear convergence rate of dADMM are presented. Numerical experiments demonstrate the effectiveness of the dADMM.

In this paper, we introduce a new augmented transformation, develop an improved partial Newton-correction method, establish and prove the close relationship between the new solutions of the boundary value problem for the fourth order nonlinear partial differential equation and their zero kernel space, remove the standard convergence assumption and make the proof more concise. We prove that our problem satisfies the conditions of the global separation theorem on Nehari manifold in different cases. The global separation theorem provides theoretical guarantee for our algorithm to find new solutions successfully. The Legendre Galerkin spectral method with interpolation projection is proposed for boundary value problems of two-dimensional nonlinear fourth order partial differential equations. By constructing interpolation operators and projection operators, the linear operators and nonlinear terms are optimized, and the algebraic equation of the original problem is obtained. It has the same condition number as the classical spectral method and achieves the same spectral accuracy. Finally, some computing results show that our method has the same high order of convergence as the classical spectral method or pseudo-spectral method, and requires less CPU time, can compute as many solutions as possible.

One-dimensional unsteady seepage model is the basic model in well-test research, existing solution includes two categories such as analytic solution and numerical solution. Analytic solution is depend on Laplace transform and numerical Inversion, and Numerical solution is depend on mesh generation and equation discretization. Iterative solution of non-linear equation is present which is direct solution for One-dimensional unsteady seepage equation. Based on steady state sequential replacement method and mass conservation equation, the unsteady seepage equation is solved by transformed as steady seepage equation, the nonlinear equation of investigation radius under different outer boundary were build, well-bore pressure and Investigation radius at each moment were obtained by Newton method. The error between iterative solution and analytic solution is small, and bigger error appeared in early stage. Compared with analytic solution and numerical solution, the iterative solution can be spread used in other One-dimensional unsteady seepage problems because of its simple algorithm, acceptable accuracy and fewer time-consuming.

Gegenbauer’s addition theorem is an important theory of Bessel functions. This paper studies the truncation error of Gegenbauer’s addition theorem. Based on the limiting forms of Bessel and Neumann functions as the argument tend to zero, an explicit bound of the truncation error and its convergence order are derived. Numerical experiments show that the estimates are valid and sharp.

In this paper, based on the active set identification, we propose a Proximal Newton method for $\ell_1$ problems. One advantage of this method is that it uses the good support identification property of the ISTA algorithm to determine the free variables and active set variables, and another advantage is that it uses the information of the partial Hessian matrix to update the free variables. Under appropriate conditions, we demonstrate the global convergence of the algorithm with a nonmonotonic line search strategy. The results of numerical experiments show that our proposed algorithm is effective.

In 2016, Karney proposed an exact sampling algorithm for the standard normal distribution. In this article, we present an exact sampling algorithm for the normal distribution of standard variance $\sqrt{1/(2\ln2)}$ and mean $0$, which can be called the binary Gaussian distribution, as its relative probability density function is given by $2^{-x^2}$ for $x\in\mathbb{R}$. Our proposed algorithm requires no floating-point arithmetic in practice, and can be regarded as the promotion of Karney's exact sampling technique. We give an estimate of the expected number of uniform deviates over the range $(0,1)$ used by this algorithm until outputting a sample value. Numerical experiments also demonstrate the effectiveness of the sampling algorithm. For any rational number $c$ greater than $1$ but less than Euler's number $e$, the idea of sampling exactly the binary Gaussian is generalized to a class of normal distributions of standard variance $\sqrt{1/(2\ln{c})}$ and mean $0$, called “$c$-ary Gaussian distributions”, and a similar complexity analysis is presented.

Multi-dimensional hypersingular integrals are widely used in many engineering fields such as elasticity and scattering of electromagnetic fields. In order to improve the calculation accuracy, we construct the formula of two-dimensional and three-dimensional hypersingular integrals. In this paper, the composite rectangle quadrature formula is used to approximate the part without singularity in the divided $N$ subinterval and the remaining part is solved by the analytic expression of the hypersingular integral. Based on the extrapolation, the modified composite rectangle quadrature formula of one-dimensional hypersingular integral is constructed. Finally, the modified rectangle quadrature formula is extended to the numerical quadrature of two-dimensional and three-dimensional surface hypersingular integrals. Numerical examples at the end of the paper verify the feasibility of the proposed method.

The high-dimensional differential and algebraic equations are common models of describing the running law of chips and electrical power systems, where the dimension of differential equations is too large to be solved effectively by using the classical numerical methods such as linear multistep methods and Runge-Kutta (RK) methods, ect. In order to solve such differential equation, the waveform relaxation (WR) method is proposed. In most cases, the differential equations mentioned above are stiff. Hence the implicit methods with very good stability property are needed to solve effectively them, and especially an A-stable method is needed. In addition, the RK method is the numerical scheme of ordinary differential equations which is employed most extensively. However, so far the works on A-stability of RK type of WR methods are not found. In this paper, we investigate A-stability of RK type of WR methods, and obtain the sufficient conditions of A-stability. The Gauss-Legendre method, Radau IA method, Radau IIA method, Lobatto IIIA method, Lobatto IIIB method and Lobatto IIIC method are A-stable Runge-Kutta methods used extensively. Under a common assumption, the obtained results show that there exist the splitting way such that the RK type of WR methods are A-stable when underlying methods are chosen as Radau IA methods or Radau IIA methods or Lobatto IIIC methods with low orders.

The adaptive finite element method (AFEM) for the P-Laplace equation in a newly defined distance was studied by D.J.Liu and Z.R. Chen ^[13] for the special case 2 ≤ p< ∞, where a posteriori error estimators with upper bound were obtained for conforming and nonconforming FEM. This paper focuses on the case when 1 < p< ∞, and the a posteriori error estimators with both upper and lower bounds are provided. The numerical investigation for benchmark problem provides the accuracy and robustness of the proposed a posteriori error estimators.

In this paper, we use the concatenating method to construct a local momentum-preserving algorithm for the Zakharov system based on the idea that the numerical algorithm should keep the intrinsic properties of the original problem as much as possible. The advantage of the proposed algorithm is that it preserves the local momentum conservation law and the local mass conservation law in any time-space regions, i.e., independence on any boundary conditions. Certainly, if the boundary condition is suitable, such as homogeneous boundary or periodic boundary, the proposed algorithm also preserves the global momentum conservation law and the global mass conservation law. Numerical experiments do not only show the second-order convergence, but also verify the stability in long term calculations and the conservation of invariants for the proposed algorithm.

In this paper, we consider the convergence analysis of two accelerated iterative algorithms for solving a nonsymmetric algebraic Riccati equation arising in transport theory. This equation has two parameters $\alpha\in[0,1), c\in(0,1]$. We prove two accelerated iterative algorithms have the same convergence rates, and show that two algorithms converge linearly in $(\alpha,c)\neq(0,1)$ and sublinearly in $(\alpha,c)=(0,1)$.

The two-subspace projection algorithm and the multi-step inertial randomized Kaczmarz algorithm are effective methods for solving the coherent linear system. By adding the soft shrinkage in the two-subspace projection algorithm and the multi-step inertial randomized Kaczmarz algorithm, this paper introduces the sparse two-subspace projection algorithm and the sparse multi-step inertial randomized Kaczmarz algorithm. The linear convergence rate in expectation of the proposed algorithms in noise and noiseless cases is presented. Finally, we give some numerical experiments to illustrate the effectiveness and advantage of our algorithms.

In this paper, we consider a inverse problem of sub-diffusion equations, i.e. the fractional parameter in the equation is optimized according to the spatial integral of the solution at a certain time. This paper shows that the inverse problem may have multiple solutions, when the initial value is 0. Based on the space-time finite element method, we propose a parameter optimization algorithm, and proved the global convergence of the algorithm by using the discrete Laplace transform technique. In addition, for the spherical region, we propose a fast algorithm by using the spectral decomposition of -Δ operator. Finally, two numerical examples verify the effectiveness of the algorithm.

The physical quantities in the electronic structure calculation of periodic systems involve the band integral over the Brillouin zone. For metals, the integrand is discontinuous when the Fermi surface crosses through the energy band, which makes it difficult to improve the calculation accuracy of general numerical integration methods. Based on smearing methods, we propose two numerical formats with higher order precision about the broadening parameters by first and second extrapolation respectively. The numerical method based on this kind of extrapolation scheme gives more precise integral approximation over the Brillouin zone with the premise of the same discrete $\mathbf{k}$-point in Brillouin region, and significantly improves the computational efficiency. The efficiency of the extrapolation methods is verified by numerical tests on two typical systems.