A Fourth Order Numerical Method for the Primitive Equations Formulated in Mean Vorticity

A fourth-order finite difference method is proposed and studied for the primitive equations (PEs) of large-scale atmospheric and oceanic flow based on mean vorticity formulation. Since the vertical average of the horizontal velocity field is divergence-free, we can introduce mean vorticity and mean stream function which are connected by a 2-D Poisson equation. As a result, the PEs can be reformulated such that the prognostic equation for the horizontal velocity is replaced by evolutionary equations for the mean vorticity field and the vertical derivative of the horizontal velocity. The mean vorticity equation is approximated by a compact difference scheme due to the difficulty of the mean vorticity boundary condition, while fourth-order long-stencil approximations are utilized to deal with transport type equations for computational convenience. The numerical values for the total velocity field (both horizontal and vertical) are statically determined by a discrete realization of a differential equation at each fixed horizontal point. The method is highly efficient and is capable of producing highly resolved solutions at a reasonable computational cost. The full fourth-order accuracy is checked by an example of the reformulated PEs with force terms. Additionally, numerical results of a large-scale oceanic circulation are presented.     

Learning to Discretize: Solving 1D Scalar Conservation Laws via Deep Reinforcement Learning

Conservation laws are considered to be fundamental laws of nature. It has broad applications in many fields, including physics, chemistry, biology, geology, and engineering. Solving the differential equations associated with conservation laws is a major branch in computational mathematics. The recent success of machine learning, especially deep learning in areas such as computer vision and natural language processing, has attracted a lot of attention from the community of computational mathematics and inspired many intriguing works in combining machine learning with traditional methods. In this paper, we are the first to view numerical PDE solvers as an MDP and to use (deep) RL to learn new solvers. As proof of concept, we focus on 1-dimensional scalar conservation laws. We deploy the machinery of deep reinforcement learning to train a policy network that can decide on how the numerical solutions should be approximated in a sequential and spatial-temporal adaptive manner. We will show that the problem of solving conservation laws can be naturally viewed as a sequential decision-making process, and the numerical schemes learned in such a way can easily enforce long-term accuracy. Furthermore, the learned policy network is carefully designed to determine a good local discrete approximation based on the current state of the solution, which essentially makes the proposed method a meta-learning approach. In other words, the proposed method is capable of learning how to discretize for a given situation mimicking human experts. Finally, we will provide details on how the policy network is trained, how well it performs compared with some state-of-the-art numerical solvers such as WENO schemes, and supervised learning based approach L3D and PINN, and how well it generalizes.     

Solution of a specific class of nonlinear fractional optimal control problems including multiple delays

This research provides a new framework based on a hybrid of block‐pulse functions and Legendre polynomials for the numerical examination of a special class of scalar nonlinear fractional optimal control problems involving delay. The concepts of the fractional derivative and the fractional integral are employed in the Caputo sense and the Riemann‐Liouville sense, respectively. In accordance with the notion of the Riemann‐Liouville integral, we derive a new integral operator related to the proposed basis called the operational matrix of fractional integration. By employing two significant operators, namely, the delay operator and the integral operator connected to the hybrid basis, the system dynamics of the primal optimal control problem converts to another system involving algebraic equations. Consequently, the optimal control problem under study is reduced to a static optimization one that is solved by existing well‐established optimization procedures. Some new theoretical results regarding the new basis are obtained. Various kinds of fractional optimal control problems containing delay are examined to measure the accuracy of the new method. The simulation results justify the merits and superiority of the devised procedure over the existing optimization methods in the literature.     

Relaxation Schemes for the $M_1$ Model with Space-Dependent Flux: Application to Radiotherapy Dose Calculation

Because of stability constraints, most numerical schemes applied to hyperbolic systems of equations turn out to be costly when the flux term is multiplied by some very large scalar. This problem emerges with the $M_1$ system of equations in the field of radiotherapy when considering heterogeneous media with very disparate densities. Additionally, the flux term of the $M_1$ system is non-linear, and in order for the model to be well-posed the numerical solution needs to fulfill conditions called realizability. In this paper, we propose a numerical method that overcomes the stability constraint and preserves the realizability property. For this purpose, we relax the $M_1$ system to obtain a linear flux term. Then we extend the stencil of the difference quotient to obtain stability. The scheme is applied to a radiotherapy dose calculation example.     

A Tailored Finite Point Method for Solving Steady MHD Duct Flow Problems with Boundary Layers

In this paper we propose a development of the finite difference method, called the tailored finite point method, for solving steady magnetohydrodynamic (MHD) duct flow problems with a high Hartmann number. When the Hartmann number is large, the MHD duct flow is convection-dominated and thus its solution may exhibit localized phenomena such as the boundary layer. Most conventional numerical methods can not efficiently solve the layer problem because they are lacking in either stability or accuracy. However, the proposed tailored finite point method is capable of resolving high gradients near the layer regions without refining the mesh. Firstly, we devise the tailored finite point method for the scalar inhomogeneous convection-diffusion problem, and then extend it to the MHD duct flow which consists of a coupled system of convection-diffusion equations. For each interior grid point of a given rectangular mesh, we construct a finite-point difference operator at that point with some nearby grid points, where the coefficients of the difference operator are tailored to some particular properties of the problem. Numerical examples are provided to show the high performance of the proposed method.     

A Distributed Optimal Control Problem with Averaged Stochastic Gradient Descent

In this work, we study a distributed optimal control problem, in which the governing system is given by second-order elliptic equations with log-normal coefficients. To lessen the curse of dimensionality that originates from the representation of stochastic coefficients, the Monte Carlo finite element method is adopted for numerical discretization where a large number of sampled constraints are involved. For the solution of such a large-scale optimization problem, stochastic gradient descent method is widely used but has slow convergence asymptotically due to its inherent variance. To remedy this problem, we adopt an averaged stochastic gradient descent method which performs stably even with the use of relatively large step sizes and small batch sizes. Numerical experiments are carried out to validate our theoretical findings.     

A Well-Posed and Discretely Stable Perfectly Matched Layer for Elastic Wave Equations in Second Order Formulation

We present a well-posed and discretely stable perfectly matched layer for the anisotropic (and isotropic) elastic wave equations without first re-writing the governing equations as a first order system. The new model is derived by the complex coordinate stretching technique. Using standard perturbation methods we show that complex frequency shift together with a chosen real scaling factor ensures the decay of eigen-modes for all relevant frequencies. To buttress the stability properties and the robustness of the proposed model, numerical experiments are presented for anisotropic elastic wave equations. The model is approximated with a stable node-centered finite difference scheme that is second order accurate both in time and space.     

An Element Decomposition Method for the Helmholtz Equation

It is well-known that the traditional full integral quadrilateral element fails to provide accurate results to the Helmholtz equation with large wave numbers due to the "pollution error" caused by the numerical dispersion. To overcome this deficiency, this paper proposed an element decomposition method (EDM) for analyzing 2D acoustic problems by using quadrilateral element. In the present EDM, the quadrilateral element is first subdivided into four sub-triangles, and the local acoustic gradient in each sub-triangle is obtained using linear interpolation function. The acoustic gradient field of the whole quadrilateral is then formulated through a weighted averaging operation, which means only one integration point is adopted to construct the system matrix. To cure the numerical instability of one-point integration, a variation gradient item is complemented by variance of the local gradients. The discretized system equations are derived using the generalized Galerkin weak form. Numerical examples demonstrate that the EDM can achieves better accuracy and higher computational efficiency. Besides, as no mapping or coordinate transformation is involved, restrictions on the shape elements can be easily removed, which makes the EDM works well even for severely distorted meshes.     

Energy‐optimal trajectory planning for the Pendubot and the Acrobot

In this paper, we study the trajectory‐planning problem for planar underactuated robot manipulators with two revolute joints in the presence of gravity. This problem is studied as an optimal control problem. We consider both possible models of planar underactuated robot manipulators, namely with the elbow joint not actuated, called Pendubot, and with the shoulder joint not actuated, called Acrobot. Our method consists of a numerical resolution of a reformulation of the optimal control problem, as an unconstrained calculus of variations problem, in which the dynamic equations of the mechanical system are regarded as constraints and treated by means of special derivative multipliers. We solve the resulting calculus of variations problem by using a numerical approach based on the Euler–Lagrange necessary condition in integral form. In which, time is discretized, and admissible variations for each variable are approximated using a linear combination of the piecewise continuous basis functions of time. In this way, a general method for the solution of optimal control problems with constraints is obtained, which, in particular, can deal with nonintegrable differential constraints arising from the dynamic models of underactuated planar robot manipulators with two revolute joints in the presence of gravity. Copyright © 2012 John Wiley & Sons, Ltd.     

Efficient Dual Algorithms for Image Segmentation Using TV-Allen-Cahn Type Models

Variational image segmentation based on the Mumford and Shah model [31], together with implementation by the piecewise constant level-set method (PCLSM) [26], leads to fully nonlinear Total Variation (TV)-Allen-Cahn equations. The commonly-used numerical approaches usually suffer from the difficulties not only with the non-differentiability of the TV-term, but also with directly evolving the discontinuous piecewise constant-structured solutions. In this paper, we propose efficient dual algorithms to overcome these drawbacks. The use of a splitting-penalty method results in TV-Allen-Cahn type models associated with different "double-well" potentials, which allow for the implementation of the dual algorithm of Chambolle [8]. Moreover, we present a new dual algorithm based on an edge-featured penalty of the dual variable, which only requires to solve a vectorial Allen-Cahn type equation with linear ∇(div)-diffusion rather than fully nonlinear diffusion in the Chambolle's approach. Consequently, more efficient numerical algorithms such as time-splitting method and Fast Fourier Transform (FFT) can be implemented. Various numerical tests show that two dual algorithms are much faster and more stable than the primal gradient descent approach, and the new dual algorithm is at least as efficient as the Chambolle's algorithm but is more accurate. We demonstrate that the new method also provides a viable alternative for image restoration.     

A Multigrid Method for a Model of the Implicit Immersed Boundary Equations

Explicit time stepping schemes for the immersed boundary method require very small time steps in order to maintain stability. Solving the equations that arise from an implicit discretization is difficult. Recently, several different approaches have been proposed, but a complete understanding of this problem is still emerging. A multigrid method is developed and explored for solving the equations in an implicit-time discretization of a model of the immersed boundary equations. The model problem consists of a scalar Poisson equation with conformation-dependent singular forces on an immersed boundary. This model does not include the inertial terms or the incompressibility constraint. The method is more efficient than an explicit method, but the efficiency gain is limited. The multigrid method alone may not be an effective solver, but when used as a preconditioner for Krylov methods, the speed-up over the explicit-time method is substantial. For example, depending on the constitutive law for the boundary force, with a time step 100 times larger than the explicit method, the implicit method is about 15-100 times more efficient than the explicit method. A very attractive feature of this method is that the efficiency of the multigrid preconditioned Krylov solver is shown to be independent of the number of immersed boundary points.     

A gradient flow approach to computing a non-linear discrete time quadratic optimal feedback gain matrix

In this paper we propose an approach to solving infinite planning horizon quadratic optimal regulator problems with linear static state feedback for discrete time systems. The approach is based on solving a sequence of approximate problems constructed by combining a finite horizon problem with an infinite horizon linear problem. A gradient-flow algorithm is derived to solve the approximate problems. As part of this, a new algorithm is derived for computing the gradient of the cost functional, based on a system of difference equations to be solved completely forward in time. Two numerical examples are presented.     

A Kernel Based Unconditionally Stable Scheme for Nonlinear Parabolic Partial Differential Equations

In this paper, a class of high order numerical schemes is proposed to solve the nonlinear parabolic equations with variable coefficients. This method is based on our previous work [11] for convection-diffusion equations, which relies on a special kernel-based formulation of the solutions and successive convolution. However, disadvantages appear when we extend the previous method to our equations, such as inefficient choice of parameters and unprovable stability for high-dimensional problems. To overcome these difficulties, a new kernel-based formulation is designed to approach the spatial derivatives. It maintains the good properties of the original one, including the high order accuracy and unconditionally stable for one-dimensional problems, hence allowing much larger time step evolution compared with other explicit schemes. In addition, without extra computational cost, the proposed scheme can enlarge the available interval of the special parameter in the formulation, leading to less errors and higher efficiency. Moreover, theoretical investigations indicate that it is unconditionally stable for multi-dimensional problems as well. We present numerical tests for one- and two-dimensional scalar and system, demonstrating the designed high order accuracy and unconditionally stable property of the scheme.     

A New Stable Version of the SPH Method in Lagrangian Coordinates

The purpose of this paper is to solve some of the trouble spots of the classical SPH method by proposing an alternative approach. First, we focus on the problem of the stability for two different SPH schemes, one is based on the approach of Vila [25] and the other is proposed in this article which mimics the classical 1D Lax Wendroff scheme. In both approaches the classical SPH artificial viscosity term is removed preserving nevertheless the linear stability of the methods, demonstrated via the von Neumann stability analysis. Moreover, the issue of the consistency for the equations of gas dynamics is analyzed. An alternative approach is proposed that consists of using Godunov-type SPH schemes in Lagrangian coordinates. This not only provides an improvement in accuracy of the numerical solutions, but also assures that the consistency condition on the gradient of the kernel function is satisfied using an equidistant distribution of particles in Lagrangian mass coordinates. Three different Riemann solvers are implemented for the first-order Godunov type SPH schemes in Lagrangian coordinates, namely the Godunov flux based on the exact Riemann solver, the Rusanov flux and a new modified Roe flux, following the work of Munz [17]. Some well-known numerical 1D shock tube test cases [22] are solved, comparing the numerical solutions of the Godunov-type SPH schemes in Lagrangian coordinates with the first-order Godunov finite volume method in Eulerian coordinates and the standard SPH scheme with Monaghan's viscosity term.     

Numerical Methods for Solving the Hartree-Fock Equations of Diatomic Molecules II

In order to solve the partial differential equations that arise in the Hartree-Fock theory for diatomic molecules and in molecular theories that include electron correlation, one needs efficient methods for solving partial differential equations. In this article, we present numerical results for a two-variable model problem of the kind that arises when one solves the Hartree-Fock equations for a diatomic molecule. We compare results obtained using the spline collocation and domain decomposition methods with third-order Hermite splines to results obtained using the more-established finite difference approximation and the successive over-relaxation method. The theory of domain decomposition presented earlier is extended to treat regions that are divided into an arbitrary number of subregions by families of lines parallel to the two coordinate axes. While the domain decomposition method and the finite difference approach both yield results at the micro-Hartree level, the finite difference approach with a 9-point difference formula produces the same level of accuracy with fewer points. The domain decomposition method has the strength that it can be applied to problems with a large number of grid points. The time required to solve a partial differential equation for a fine grid with a large number of points goes down as the number of partitions increases. The reason for this is that the length of time necessary for solving a set of linear equations in each subregion is very much dependent upon the number of equations. Even though a finer partition of the region has more subregions, the time for solving the set of linear equations in each subregion is very much smaller. This feature of the theory may well prove to be a decisive factor for solving the two-electron pair equation, which – for a diatomic molecule – involves solving partial differential equations with five independent variables. The domain decomposition theory also makes it possible to study complex molecules by dividing them into smaller fragments thatare calculated independently. Since the domain decomposition approach makes it possible to decompose the variable space into separate regions in which the equations are solved independently, this approach is well-suited to parallel computing.     

Computation of Shock Profiles in Radiative Hydrodynamics

This article is devoted to the construction of a numerical scheme to solve the equations of radiative hydrodynamics. We use this numerical procedure to compute shock profiles and illustrate some earlier theoretical results about their smoothness and monotonicity properties. We first consider a scalar toy model, then we extend our analysis to a more realistic system for the radiative hydrodynamics that couples the Euler equations and an elliptic equation.     

Convergent Adaptive Finite Element Method Based on Centroidal Voronoi Tessellations and Superconvergence

We present a novel adaptive finite element method (AFEM) for elliptic equations which is based upon the Centroidal Voronoi Tessellation (CVT) and superconvergent gradient recovery. The constructions of CVT and its dual Centroidal Voronoi Delaunay Triangulation (CVDT) are facilitated by a localized Lloyd iteration to produce almost equilateral two dimensional meshes. Working with finite element solutions on such high quality triangulations, superconvergent recovery methods become particularly effective so that asymptotically exact a posteriori error estimations can be obtained. Through a seamless integration of these techniques, a convergent adaptive procedure is developed. As demonstrated by the numerical examples, the new AFEM is capable of solving a variety of model problems and has great potential in practical applications.     

Analysis and Numerical Simulation of Hyperbolic Shallow Water Moment Equations

Shallow Water Moment Equations allow for vertical changes in the horizontal velocity, so that complex shallow flows can be described accurately. However, we show that these models lack global hyperbolicity and the loss of hyperbolicity already occurs for small deviations from equilibrium. This leads to instabilities in a numerical test case. We then derive new Hyperbolic Shallow Water Moment Equations based on a modification of the system matrix. The model can be written in analytical form and hyperbolicity can be proven for a large number of equations. A second variant of this model is obtained by generalizing the modification with the help of additional parameters. Numerical tests of a smooth periodic problem and a dam break problem using the new models yield accurate and fast solutions while guaranteeing hyperbolicity.     

Optimal control of radiative heat transfer in glass cooling with restrictions on the temperature gradient

This paper is motivated by an optimal boundary control problem for the cooling process of molten and already formed glass down to room temperature. The high temperatures at which glass is processed demand to include radiative heat transfer in the computational model. Since the complete radiative heat transfer equations are too complex for optimization purposes, we use simplified approximations of spherical harmonics coupled with a practically relevant frequency bands model. The optimal control problem is considered as a partial differential algebraic equation (PDAE)‐constrained optimization problem with box constraints on the control. In this paper, we augment the objective by a functional depending on the state gradient, which forces a minimization of thermal stress inside the glass. To guarantee consistent and grid‐independent values of the reduced objective gradient at the end of the cooling process, we pursue two approaches. The first includes the temperature gradient with a time‐dependent linearly decreasing weight. In the second approach, we augment the objective functional by the final state tracking and final state gradient optimization. To determine an optimal boundary control, we apply a projected gradient method with the Armijo step size rule. The reduced objective gradient is computed by the continuous adjoint approach. The arising time‐dependent PDAEs are numerically solved by variable step size one‐step methods of Rosenbrock type in time and adaptive multilevel finite elements in space. We present two‐dimensional numerical results for an infinitely long glass block and compare the two different approaches derived to ensure consistency at the end of the cooling process. Copyright © 2011 John Wiley & Sons, Ltd.     

One-Dimensional Blood Flow with Discontinuous Properties and Transport: Mathematical Analysis and Numerical Schemes

In this paper we consider the one-dimensional blood flow model with discontinuous mechanical and geometrical properties, as well as passive scalar transport, proposed in [E.F. Toro and A. Siviglia. Flow in collapsible tubes with discontinuous mechanical properties: mathematical model and exact solutions. Communications in Computational Physics. 13(2), 361-385, 2013], completing the mathematical analysis by providing new propositions and new proofs of relations valid across different waves. Next we consider a first order DOT Riemann solver, proposing an integration path that incorporates the passive scalar and proving the well-balanced properties of the resulting numerical scheme for stationary solutions. Finally we describe a novel and simple well-balanced, second order, non-linear numerical scheme to solve the equations under study; by using suitable test problems for which exact solutions are available, we assess the well-balanced properties of the scheme, its capacity to provide accurate solutions in challenging flow conditions and its accuracy.