Zhengyang XueUniversity of Science and Technology of ChinaYinhua XiaUniversity of Science and Technology of ChinaChen LiState Key Laboratory of AerodynamicsXianxu YuanState Key Laboratory of Aerodynamics
Numerical Analysis and Scientific Computingmathscidoc:2208.25003
Guillaume BalDepartment of Applied Physics and Applied Mathematics, Columbia University, 10027 New York, USA.Wenjia JingD´epartement de Math´ematiques et Applications, Ecole Normale Sup´erieure, 45 Rue d’Ulm, 75230 Paris Cedex 05, France
Analysis of PDEsNumerical Analysis and Scientific ComputingProbabilitymathscidoc:2206.03008
This paper analyzes the random fluctuations obtained by a heterogeneous multi-scale first-order finite element method applied to solve elliptic equations with a random potential. Several multiscale numerical algorithms have been shown to correctly capture the homogenized limit of solutions of elliptic equations with coefficients modeled as stationary and ergodic random fields. Because theoretical results are available in the continuum setting for such equations, we consider here the case of a secondorder elliptic equations with random potential in two dimensions of space. We show that the random fluctuations of such solutions are correctly estimated by the heterogeneous multi-scale algorithm when appropriate fine-scale problems are solved on subsets that cover the whole computational domain. However, when the fine-scale problems are solved over patches that do not cover the entire domain, the random fluctuations may or may not be estimated accurately. In the case of random potentials with short-range interactions, the variance of the random fluctuations is amplified as the inverse of the fraction of the medium covered by the patches. In the case of random potentials with long-range interactions, however, such an amplification does not occur and random fluctuations are correctly captured independent of the (macroscopic) size of the patches. These results are consistent with those obtained in  for more general equations in the one-dimensional setting and provide indications on the loss in accuracy that results from using coarser, and hence computationally less intensive, algorithms.
Habib AmmariDepartment of Mathematics and Applications, Ecole Normale Supérieure, 45 Rue d’Ulm, 75005, Paris, FranceThomas BoulierDepartment of Mathematics and Applications, Ecole Normale Supérieure, 45 Rue d’Ulm, 75005, Paris, FranceJosselin GarnierLaboratoire de Probabilités et Modèles Aléatoires & Laboratoire Jacques-Louis Lions, Université Paris VII, 75205, Paris Cedex 13, FranceWenjia JingDepartment of Mathematics and Applications, Ecole Normale Supérieure, 45 Rue d’Ulm, 75005, Paris, FranceHyeonbae KangDepartment of Mathematics, Inha University, Incheon, 402-751, KoreaHan WangDepartment of Mathematics and Applications, Ecole Normale Supérieure, 45 Rue d’Ulm, 75005, Paris, France
Analysis of PDEsMathematical PhysicsNumerical Analysis and Scientific ComputingOptimization and Controlmathscidoc:2206.03006
Foundations of Computational Mathematics, 14, 27-62, 2013.9
The aim of this paper is to provide a fast and efficient procedure for (real-time) target identification in imaging based on matching on a dictionary of precomputed generalized polarization tensors (GPTs). The approach is based on some important properties of the GPTs and new invariants. A new shape representation is given and numerically tested in the presence of measurement noise. The stability and resolution of the proposed identification algorithm is numerically quantified. We compare the proposed GPT-based shape representation with a moment-based one.
Guillaume BalDepartment of Applied Physics and Applied Mathematics, Columbia University, New York 10027Wenjia JingDépartement de Mathématiques et Applications, Ecole Normale Supérieure, 45 Rue d’Ulm, 75230 ParisCedex 05, France
Numerical Analysis and Scientific ComputingProbabilitymathscidoc:2206.25009
We analyze the random fluctuations of several multiscale algorithms, such as the multiscale finite element method (MsFEM) and the finite element heterogeneous multiscale method (HMM), that have been developed to solve partial differential equations with highly heterogeneous coefficients. Such multiscale algorithms are often shown to correctly capture the homogenization limit when the highly oscillatory random medium is stationary and ergodic. This paper is concerned with the random fluctuations of the solution about the deterministic homogenization limit. We consider the simplified setting of the one-dimensional elliptic equation, where the theory of random fluctuations is well understood. We develop a fluctuation theory for the multiscale algorithms in the presence of random environments with short-range and long-range correlations. For a given mesh size h, we show that the fluctuations converge in distribution in the space of continuous paths to Gaussian processes as the correlation length ε→0. We next derive the limit of such Gaussian processes as h→0 and compare this limit with the distribution of the random fluctuations of the continuous model. When such limits agree, we conclude that the multiscale algorithm captures the random fluctuations accurately and passes the corrector test. This property serves as an interesting benchmark to assess the behavior of the multiscale algorithm in practical situations where the assumptions necessary for the theory of homogenization are not met. What we find is that the computationally more expensive methods MsFEM, and HMM with a choice of parameter δ=h, correctly capture the random fluctuations both for short-range and long-range oscillations in the medium. The less expensive method HMM with δ<h correctly captures the fluctuations for long-range oscillations and strongly amplifies their size in media with short-range oscillations. We present a modified scheme with an intermediate computational cost that captures the random fluctuations in all cases.
In this article, we propose a numerical method based on sparse Gaussian processes (SGPs) to solve nonlinear partial differential equations (PDEs). The SGP algorithm is based on a Gaussian process (GP) method, which approximates the solution of a PDE with the maximum a posteriori probability estimator of a GP conditioned on the PDE evaluated at a finite number of sample points. The main bottleneck of the GP method lies in the inversion of a covariance matrix, whose cost grows cubically with respect to the size of samples. To improve the scalability of the GP method while retaining desirable accuracy, we draw inspiration from SGP approximations, where inducing points are introduced to summarize the information of samples. More precisely, our SGP method uses a Gaussian prior associated with a low-rank kernel generated by inducing points randomly selected from samples. In the SGP method, the size of the matrix to be inverted is proportional to the number of inducing points, which is much less than the size of the samples. The numerical experiments show that the SGP method using less than half of the uniform samples as inducing points achieves comparable accuracy to the GP method using the same number of uniform samples, which significantly reduces the computational cost. We give the existence proof for the approximation to the solution of a PDE and provide rigorous error analysis.
Chenchen MouDepartment of Mathematics, City University of Hong Kong, Hong Kong, ChinaXianjin YangYau Mathematical Sciences Center, Tsinghua University, Haidian District, Beijing, 100084, China; Beijing Institute of Mathematical Sciences and Applications, Huairou District, Beijing, 101408, ChinaChao ZhouDepartment of Mathematics and Risk Management Institute, National University of Singapore, Singapore
Numerical Analysis and Scientific ComputingOptimization and Controlmathscidoc:2206.25007
Journal of Computational Physics, 460, (1), 111188, 2022.7
In this article, we propose two numerical methods, the Gaussian Process (GP) method and the Fourier Features (FF) algorithm, to solve mean field games (MFGs). The GP algorithm approximates the solution of a MFG with maximum a posteriori probability estimators of GPs conditioned on the partial differential equation (PDE) system of the MFG at a finite number of sample points. The main bottleneck of the GP method is to compute the inverse of a square gram matrix, whose size is proportional to the number of sample points. To improve the performance, we introduce the FF method, whose insight comes from the recent trend of approximating positive definite kernels with random Fourier features. The FF algorithm seeks approximated solutions in the space generated by sampled Fourier features. In the FF method, the size of the matrix to be inverted depends only on the number of Fourier features selected, which is much less than the size of sample points. Hence, the FF method reduces the precomputation time, saves the memory, and achieves comparable accuracy to the GP method. We give the existence and the convergence proofs for both algorithms. The convergence argument of the GP method does not depend on any monotonicity condition, which suggests the potential applications of the GP method to solve MFGs with non-monotone couplings in future work. We show the efficacy of our algorithms through experiments on a stationary MFG with a non-local coupling and on a time-dependent planning problem. We believe that the FF method can also serve as an alternative algorithm to solve general PDEs.
Kai JiangSchool of Mathematics and Computational Science, Xiangtan University, Xiangtan, Hunan, P.R. China, 411105Wei SiSchool of Mathematics and Computational Science, Xiangtan University, Xiangtan, Hunan, P.R. China, 411105Chenglong BaoYau Mathematical Sciences Center, Tsinghua University, Beijing, P. R. China, 100084
Numerical Analysis and Scientific Computingmathscidoc:2206.25006
Computing stationary states is an important topic for phase field crystal (PFC) models. Great efforts have been made for energy dissipation of the numerical schemes when using gradient flows. However, it is always time-consuming due to the requirement of small effective time steps. In this paper, we propose an adaptive accelerated proximal gradient method for finding the stationary states of PFC models. The energy dissipation is guaranteed and the convergence property is established for the discretized energy functional. Moreover, the connections between generalized proximal operator with classical (semi-)implicit and explicit schemes for gradient flow are given. Extensive numerical experiments, including two three dimensional periodic crystals in Landau-Brazovskii (LB) model and a two dimensional quasicrystal in Lifshitz-Petrich (LP) model, demonstrate that our approach has adaptive time steps which lead to significant acceleration over semi-implicit methods for computing complex structures. Furthermore, our result reveals a deep physical mechanism of the simple LB model via which the sigma phase is first discovered.
Chenglong BaoDepartment of Mathematics, National University of Singapore, Singapore 117543, SingaporeHui JiDepartment of Mathematics, National University of Singapore, Singapore 117543, SingaporeZuowei ShenDepartment of Mathematics, National University of Singapore, Singapore 117543, Singapore
Numerical Analysis and Scientific Computingmathscidoc:2206.25005
Applied and Computational Harmonic Analysis, 38, (3), 510-523, 2015.5
Sparse modeling/approximation of images plays an important role in image restoration. Instead of using a fixed system to sparsely model any input image, a more promising approach is using a system that is adaptive to the input image. A non-convex variational model is proposed in  for constructing a tight frame that is optimized for the input image, and an alternating scheme is used to solve the resulting non-convex optimization problem. Although it showed good empirical performance in image denoising, the proposed alternating iteration lacks convergence analysis. This paper aims at providing the convergence analysis of the method proposed in . We first established the sub-sequence convergence property of the iteration scheme proposed in , i.e., there exists at least one convergent sub-sequence and any convergent sub-sequence converges to a stationary point of the minimization problem. Moreover, we showed that the original method can be modified to have sequence convergence, i.e., the modified algorithm generates a sequence that converges to a stationary point of the minimization problem.
Jae Kyu ChoiInstitute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240, China Chenglong BaoYau Mathematical Sciences Center, Tsinghua University, Beijing 100084, ChinaXiaoqun ZhangInstitute of Natural Sciences, School of Mathematical Sciences, and MOE-LSC, Shanghai Jiao Tong University, Shanghai 200240, China
Numerical Analysis and Scientific Computingmathscidoc:2206.25004
SIAM Journal on Imaging Sciences, 11, (2), 1179-1204, 2018.5
Recent technical advances lead to the coupling of PET and MRI scanners, enabling one to acquire functional and anatomical data simultaneously. In this paper, we propose a tight frame based PET-MRI joint reconstruction model via the joint sparsity of tight frame coefficients. In addition, a nonconvex balanced approach is adopted to take the different regularities of PET and MRI images into account. To solve the nonconvex and nonsmooth model, a proximal alternating minimization algorithm is proposed, and the global convergence is present based on the Kurdyka--Łojasiewicz property. Finally, the numerical experiments show that our proposed models achieve better performance over the existing PET-MRI joint reconstruction models.
Chenglong BaoYau Mathematical Sciences Center, Tsinghua University, Beijing, 100084 China Jae Kyu ChoiSchool of Mathematical Sciences, Tongji University, Shanghai, 200092 ChinaBin DongBeijing International Center for Mathematical Research and Laboratory for Biomedical Image Analysis, Beijing Institute of Big Data Research, Peking University, Beijing, 100871 China
Numerical Analysis and Scientific Computingmathscidoc:2206.25003
SIAM Journal on Imaging Sciences, 12, (1), 492-520, 2019.2
Quantitative susceptibility mapping (QSM) uses the phase data in magnetic resonance signals to visualize a three-dimensional susceptibility distribution by solving the magnetic field to susceptibility inverse problem. Due to the presence of zeros of the integration kernel in the frequency domain, QSM is an ill-posed inverse problem. Although numerous regularization-based models have been proposed to overcome this problem, incompatibility in the field data, which leads to deterioration of the recovery, has not received enough attention. In this paper, we show that the data acquisition process of QSM inherently generates a harmonic incompatibility in the measured local field. Based on this discovery, we propose a novel regularization-based susceptibility reconstruction model with an additional sparsity-based regularization term on the harmonic incompatibility. Numerical experiments show that the proposed method achieves better performance than existing approaches.
Kai JiangSchool of Mathematics and Computational Science, Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, Xiangtan, Hunan, China, 411105Wei SiSchool of Mathematics and Computational Science, Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, Xiangtan, Hunan, China, 411105Chang ChenSchool of Mathematics and Computational Science, Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, Xiangtan, Hunan, China, 411105Chenglong BaoYau Mathematical Sciences Center, Tsinghua University, Beijing, China, 100084
Numerical Analysis and Scientific Computingmathscidoc:2206.25002
Finding the stationary states of a free energy functional is an important problem in phase field crystal (PFC) models. Many efforts have been devoted to designing numerical schemes with energy dissipation and mass conservation properties. However, most existing approaches are time-consuming due to the requirement of small effective step sizes. In this paper, we discretize the energy functional and propose efficient numerical algorithms for solving the constrained nonconvex minimization problem. A class of gradient-based approaches, which are the so-called adaptive accelerated Bregman proximal gradient (AA-BPG) methods, is proposed, and the convergence property is established without the global Lipschitz constant requirements. A practical Newton method is also designed to further accelerate the local convergence with convergence guarantee. One key feature of our algorithms is that the energy dissipation and mass conservation properties hold during the iteration process. Moreover, we develop a hybrid acceleration framework to accelerate the AA-BPG methods and most of the existing approaches through coupling with the practical Newton method. Extensive numerical experiments, including two three-dimensional periodic crystals in the Landau--Brazovskii (LB) model and a two-dimensional quasicrystal in the Lifshitz--Petrich (LP) model, demonstrate that our approaches have adaptive step sizes which lead to a significant acceleration over many existing methods when computing complex structures.
Chenglong BaoYau Mathematical Sciences Center, Tsinghua University, Beijing, 100084, China; Yanqi Lake Beijing Institute of Mathematical Sciences and Applications, ChinaChang ChenYau Mathematical Sciences Center, Tsinghua University, Beijing, 100084, ChinaKai JiangSchool of Mathematics and Computational Science, Hunan Key Laboratory for Computation and Simulation in Science and Engineering, Xiangtan University, Xiangtan, Hunan, 411105, China
Numerical Analysis and Scientific Computingmathscidoc:2206.25001
In this paper, we compute the stationary states of the multicomponent phase-field crystal model by formulating it as a block constrained minimization problem. The original infinite-dimensional non-convex minimization problem is approximated by a finite-dimensional constrained non-convex minimization problem after an appropriate spatial discretization. To efficiently solve the above optimization problem, we propose a so-called adaptive block Bregman proximal gradient (AB-BPG) algorithm that fully exploits the problem's block structure. The proposed method updates each order parameter alternatively, and the update order of blocks can be chosen in a deterministic or random manner. Besides, we choose the step size by developing a practical linear search approach such that the generated sequence either keeps energy dissipation or has a controllable subsequence with energy dissipation. The convergence property of the proposed method is established without the requirement of global Lipschitz continuity of the derivative of the bulk energy part by using the Bregman divergence. The numerical results on computing stationary ordered structures in binary, ternary, and quinary component coupled-mode Swift-Hohenberg models have shown a significant acceleration over many existing methods.
Jie DuYau Mathematical Sciences Center, Tsinghua University, Beijing, 100084, China; Yanqi Lake Beijing Institute of Mathematical Sciences and Applications, Beijing, 101408, ChinaYang YangDepartment of Mathematical Sciences, Michigan Technological University, Houghton, MI, 49931, USA
Numerical Analysis and Scientific Computingmathscidoc:2205.25022
Communications on Applied Mathematics and Computation, 2021.4
In this paper, we apply high-order finite difference (FD) schemes for multispecies and multireaction detonations (MMD). In MMD, the density and pressure are positive and the mass fraction of the ith species in the chemical reaction, say zi, is between 0 and 1, with ∑z_i=1. Due to the lack of maximum-principle, most of the previous bound-preserving technique cannot be applied directly. To preserve those bounds, we will use the positivity-preserving technique to all the z′is and enforce ∑z_i=1 by constructing conservative schemes, thanks to conservative time integrations and consistent numerical fluxes in the system. Moreover, detonation is an extreme singular mode of flame propagation in premixed gas, and the model contains a significant stiff source. It is well known that for hyperbolic equations with stiff source, the transition points in the numerical approximations near the shocks may trigger spurious shock speed, leading to wrong shock position. Intuitively, the high-order weighted essentially non-oscillatory (WENO) scheme, which can suppress oscillations near the discontinuities, would be a good choice for spatial discretization. However, with the nonlinear weights, the numerical fluxes are no longer “consistent”, leading to nonconservative numerical schemes and the bound-preserving technique does not work. Numerical experiments demonstrate that, without further numerical techniques such as subcell resolutions, the conservative FD method with linear weights can yield better numerical approximations than the nonconservative WENO scheme.
Jie DuYau Mathematical Sciences Center, Tsinghua University, Beijing, 100084, China; Yanqi Lake Beijing Institute of Mathematical Sciences and Applications, Beijing, 101408, ChinaEric ChungDepartment of Mathematics, The Chinese University of Hong Kong, Hong Kong SAR, ChinaYang YangDepartment of Mathematical Sciences, Michigan Technological University, Houghton, MI, 49931, USA
Numerical Analysis and Scientific Computingmathscidoc:2205.25021
Communications on Applied Mathematics and Computation, 4, 353-379, 2022.4
In this paper, we study the classical Allen-Cahn equations and investigate the maximum-principle-preserving (MPP) techniques. The Allen-Cahn equation has been widely used in mathematical models for problems in materials science and fluid dynamics. It enjoys the energy stability and the maximum-principle. Moreover, it is well known that the Allen-Cahn equation may yield thin interface layer, and nonuniform meshes might be useful in the numerical solutions. Therefore, we apply the local discontinuous Galerkin (LDG) method due to its flexibility on h-p adaptivity and complex geometry. However, the MPP LDG methods require slope limiters, then the energy stability may not be easy to obtain. In this paper, we only discuss the MPP technique and use numerical experiments to demonstrate the energy decay property. Moreover, due to the stiff source given in the equation, we use the conservative modified exponential Runge-Kutta methods and thus can use relatively large time step sizes. Thanks to the conservative time integration, the bounds of the unknown function will not decay. Numerical experiments will be given to demonstrate the good performance of the MPP LDG scheme.
Jie DuYau Mathematical Sciences Center, Tsinghua University, Beijing, ChinaEric T. ChungDepartment of Mathematics, The Chinese University of Hong Kong, Sha Tin, Hong Kong SARMing Fai LamDepartment of Mathematics, The Chinese University of Hong Kong, Sha Tin, Hong Kong SARXiao-Ping WangDepartment of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong SAR
Numerical Analysis and Scientific Computingmathscidoc:2205.25020
Journal of Scientific Computing, 76, 1547–1577, 2018.3
In this paper, we present a discontinuous Galerkin method with staggered hybridization to discretize a class of nonlinear Stokes equations in two dimensions. The utilization of staggered hybridization is new and this approach combines the features of traditional hybridization method and staggered discontinuous Galerkin method. The main idea of our method is to use hybrid variables to impose the staggered continuity conditions instead of enforcing them in the approximation space. Therefore, our method enjoys some distinctive advantages, including mass conservation, optimal convergence and preservation of symmetry of the stress tensor. We will also show that, one can obtain superconvergent and strongly divergence-free velocity by applying a local postprocessing technique on the approximate solution. We will analyze the stability and derive a priori error estimates of the proposed scheme. The resulting nonlinear system is solved by using the Newton’s method, and some numerical results will be demonstrated to confirm the theoretical rates of convergence and superconvergence.
Eric T. ChungDepartment of Mathematics, The Chinese University of Hong Kong, Hong Kong, SAR, People’s Republic of ChinaJie DuDepartment of Mathematics, The Chinese University of Hong Kong, Hong Kong, SAR, People’s Republic of ChinaMan Chun YuenDepartment of Mathematics, The Chinese University of Hong Kong, Hong Kong, SAR, People’s Republic of China
Numerical Analysis and Scientific Computingmathscidoc:2205.25019
Journal of Scientific Computing, 70, 766-792, 2016.8
Staggered grid techniques are attractive ideas for flow problems due to their more enhanced conservation properties. Recently, a staggered discontinuous Galerkin method is developed for the Stokes system. This method has several distinctive advantages, namely high order optimal convergence as well as local and global conservation properties. In addition, a local postprocessing technique is developed, and the postprocessed velocity is superconvergent and pointwisely divergence-free. Thus, the staggered discontinuous Galerkin method provides a convincing alternative to existing schemes. For problems with corner singularities and flows in porous media, adaptive mesh refinement is crucial in order to reduce the computational cost. In this paper, we will derive a computable error indicator for the staggered discontinuous Galerkin method and prove that this indicator is both efficient and reliable. Moreover, we will present some numerical results with corner singularities and flows in porous media to show that the proposed error indicator gives a good performance.