Central Finite Difference Research Articles

A new eighth-order implicit finite difference method to solve the three-dimensional Helmholtz equation

This paper introduces a novel method for deducing high-order finite difference schemes for the three-dimensional Helmholtz equation, Δp±λ2p=f. The new method can reach an eighth-order of accuracy. It is based on an implicit formulation by computing simultaneously the unknown variable and its corresponding derivatives. The scheme is obtained from Taylor series expansion and plane wave theory analysis, and it is constructed from a few modifications to the standard centered finite differences. A classical dispersion analysis and the error between the numerical wavenumber and the exact wavenumber is given. Numerical experiments are presented to verify the feasibility and accuracy of the proposed methods for a broad range of λ-values. Moreover, the numerical method is capable to solve with great precision the monotone and oscillatory Helmholtz problem with variable λ-values. Thus, the proposed formulation results in an attractive method, easy to implement, with high order of accuracy but nearly the same computational cost as those of an explicit formulation.

A two-dimensional finite element recursion relation for the transport equation using nine-diagonal solvers

A Galerkin-based finite element recursion relation is used to solve the heat transport equation in two-dimensions. The finite element method (FEM) is a powerful technique that is commonly used for solving complex engineering problems. However, the implementation of the FEM in multi-dimensional problems can be computationally expensive. A finite element recursion algorithm based on bilinear triangular, bilinear quadrilateral and quadratic Lagrangian approximations are employed to discretize the 2-D advection-diffusion equation. This algorithm is an extension of the 1-D Chapeau (linear element) technique, which employed a tridiagonal recursion expression common to the classical central finite-difference approach. The global matrix is nine-diagonal (for 2-D) and is solved using a modified strongly implicit procedure and a left-to-right sweep method.

Constructing a subgradient from directional derivatives for functions of two variables

For any scalar-valued bivariate function that is locally Lipschitz continuous and directionally differentiable, it is shown that a subgradient may always be constructed from the function's directional derivatives in the four compass directions, arranged in a so-called "compass difference". When the original function is nonconvex, the obtained subgradient is an element of Clarke's generalized gradient, but the result appears to be novel even for convex functions. The function is not required to be represented in any particular form, and no further assumptions are required, though the result is strengthened when the function is additionally L-smooth in the sense of Nesterov. For certain optimal-value functions and certain parametric solutions of differential equation systems, these new results appear to provide the only known way to compute a subgradient. These results also imply that centered finite differences will converge to a subgradient for bivariate nonsmooth functions. As a dual result, we find that any compact convex set in two dimensions contains the midpoint of its interval hull. Examples are included for illustration, and it is demonstrated that these results do not extend directly to functions of more than two variables or sets in higher dimensions.

Exponential Runge–Kutta Method for Two-Dimensional Nonlinear Fractional Complex Ginzburg–Landau Equations

In this work, we study numerically two-dimensional nonlinear spatial fractional complex Ginzburg–Landau equations. A centered finite difference method is exploited to discretize the spatial variables and leads to a system of the ordinary differential equation, in which the resulting coefficient matrix is complex symmetric and possesses the block Toeplitz structure. An exponential Runge–Kutta method is employed to solve such a system of the ordinary differential equation. Theoretically, the proposed method is second-order accuracy in space and fourth-order accuracy in time, respectively. In the practical implementation, the product of a block Toeplitz matrix exponential and a vector is calculated by the shift-invert Lanczos method. Meanwhile, the sectorial operator (the coefficient matrix) guarantees the fast approximation by the shift-invert Lanczos method. Numerical experiments are carried out to testify the theoretical results and demonstrate that the proposed method enjoys the excellent computational advantage.

Acoustic impedance sensitivity approach for defining absorbing material layout using the wave based method

This study used the wave based method (WBM) to investigate the optimum layout of an absorbing material on a vibrating surface. The structural dynamic behavior was examined using the finite element model, and the acoustic response was computed using the WBM. The thickness of the absorbing material on each structural element was considered as the design variable. The sensitivity of the acoustic response to the acoustic impedance was derived using wave based system matrices, and a polynomial function was used to describe the impedance variation with respect to the thickness of the absorbing material. Thereafter, the final design sensitivity of a vibro-acoustic model was derived with respect to the thickness of the damping material using the chain rule. The design sensitivity results were found to be in good agreement with the results of the central finite difference methods. The design optimization of the absorbing material layout was applied to a vibrating surface using the design sensitivity results under various excitations, and the results of the optimization showed a layout similar to the active intensity plot of the vibrating surface. The proposed method can easily be applied to reduce the noise produced by a vibrating structure.

CONUNDrum: A program for orbital-free density functional theory calculations

We present a new code for energy minimization, structure relaxation and evaluation of bulk parameters in the framework of orbital-free density functional theory (OF-DFT). The implementation is based on solving the Euler–Lagrange equation on an equidistant real space grid on which density dependent variables and derivatives are computed. Some potential components are computed in Fourier space. The code is able to use semilocal and non-local kinetic energy functionals (KEF) as well as neural network based KEFs thus facilitating testing and development of emerging machine-learned KEFs. For semi-local and machine-learned KEFs the kinetic energy potentials are evaluated with real-space differentiation of the components, which are partial derivatives of the KE with respect to the electron density, its gradient and Laplacian. Program summaryProgram title: CONUNDrum.CPC Library link to program files:http://dx.doi.org/10.17632/phnz2gg8mz.1Licensing provision: GNU GPL v3Programming language: C++External routines: Fastest Fourier Transform in the West (FFTW) library (http://www.fftw.org/)Nature of problem: Calculation of the electronic and structural properties of molecules and extended systems in the framework of the orbital-free density functional theory. Evaluation of the bulk parameters of solid compounds.Solution method: High-order central finite-difference method and fast Fourier transform are used for calculation of different total energy components. Density optimization is performed with the steepest descent or the Polak–Ribière variant of the non-linear conjugate-gradient method with a line search procedure based on the Armijo condition. A numerical approach is used for structural optimization — the total energies with respect to small variations in lattice geometries are computed directly, with subsequent evaluation of the force components via a high-order central-finite difference method. The same numerical procedure is used for evaluation of bulk properties.Restrictions: Local pseudopotentials.

Diagonal Approximation of the Hessian by Finite Differences for Unconstrained Optimization

A new quasi-Newton method with a diagonal updating matrix is suggested, where the diagonal elements are determined by forward or by central finite differences. The search direction is a direction of sufficient descent. The algorithm is equipped with an acceleration scheme. The convergence of the algorithm is linear. The preliminary computational experiments use a set of 75 unconstrained optimization test problems classified in five groups according to the structure of their Hessian: diagonal, block-diagonal, band (tri- or penta-diagonal), sparse and dense. Subject to the CPU time metric, intensive numerical experiments show that, for problems with Hessian in a diagonal, block-diagonal or band structure, the algorithm with diagonal approximation of the Hessian by finite differences is top performer versus the well-established algorithms: the steepest descent and the Broyden–Fletcher–Goldfarb–Shanno. On the other hand, as a by-product of this numerical study, we show that the Broyden–Fletcher–Goldfarb–Shanno algorithm is faster for problems with sparse Hessian, followed by problems with dense Hessian.

Balancing aspects of numerical dissipation, dispersion, and aliasing in time‐accurate simulations

SummaryThe current study looks at the selection of scheme elements that are well‐suited for long‐time integration of unsteady flows in the absence or under‐resolution of physical diffusion. A concerted assembly of numerical components are chosen relative to a target aliasing limit, which is taken as a best‐case scenario for overall spectral resolvability. High‐order and optimized difference stencils are employed in order to achieve accuracy; meanwhile, quasi skew‐symmetric splitting techniques for nonlinear transport terms are used in order to greatly improve robustness. Finally, tunable and scale‐discriminant artificial‐dissipation methods are incorporated for de‐aliasing purposes and as a means of further enhancing both accuracy and stability. Central finite difference methods are considered, and spectral characterizations of the scheme components are presented. Canonical test cases (the isentropic vortex [IV] and Taylor‐Green vortex problems) are chosen in order to highlight the benefits associated with the proposed approach for enhancing overall algorithm robustness and accuracy.

A numerical approach for assessing flow compressibility and transonic effect on airfoil aerodynamics

For free Mach numbers M∞ < 0.65, the standard methods for subsonic potential flows break down when the flow becomes supersonic locally. The transition of the potential equation from elliptic to hyperbolic indicates that the flow changes from a diffusive character to a propagation behavior. The elliptic numerical operators for the subsonic flow will not be able to correctly simulate the propagation properties of supersonic flow regions, which are governed by the equation of hyperbolic type. The presence of local supersonic regions significantly influences the aerodynamic characteristics of airfoil when compared to the predicted results using the standard methods for subsonic flows. To evaluate this transonic effect, the paper presents a method of solving the full potential equation using central and backward finite difference schemes corresponding to the subsonic and supersonic regions. The rotated difference scheme consists in differencing all the derivatives in the normal direction centrally, while the derivatives in the local flow direction are upwind differenced at supersonic points. The transonic effects on airfoil aerodynamics depending on the free Mach number, incidence angle and airfoil geometry were evaluated.

Mode shape curvature of multiple cracked beam and its use for crack identification in beam-like structures

The problem of using the modal curvature for crack detection is discussed in this paper based on an exact expression of mode shape and its curvature. Using the obtained herein exact expression for the mode shape and its curvature, it is demonstrated that the mode shape curvature is really more sensitive to crack than mode shape itself. Nevertheless, crack-induced change in the approximate curvature calculated from the exact mode shape by the central finite difference technique (Laplacian) is much greater in comparison with both the mode shape and curvature. It is produced by the fact, shown in this study, that miscalculation of the approximate curvature is straightforwardly dependent upon crack magnitude and resolution step of the finite difference approximation. Therefore, it can be confidently recommended to use the approximate curvature for multiple crack detection in beam by properly choosing the approximation mesh. The theoretical development has been illustrated by numerical results.

Simulating coupled surface–subsurface flows with ParFlow v3.5.0: capabilities, applications, and ongoing development of an open-source, massively parallel, integrated hydrologic model

Abstract. Surface flow and subsurface flow constitute a naturally linked hydrologic continuum that has not traditionally been simulated in an integrated fashion. Recognizing the interactions between these systems has encouraged the development of integrated hydrologic models (IHMs) capable of treating surface and subsurface systems as a single integrated resource. IHMs are dynamically evolving with improvements in technology, and the extent of their current capabilities are often only known to the developers and not general users. This article provides an overview of the core functionality, capability, applications, and ongoing development of one open-source IHM, ParFlow. ParFlow is a parallel, integrated, hydrologic model that simulates surface and subsurface flows. ParFlow solves the Richards equation for three-dimensional variably saturated groundwater flow and the two-dimensional kinematic wave approximation of the shallow water equations for overland flow. The model employs a conservative centered finite-difference scheme and a conservative finite-volume method for subsurface flow and transport, respectively. ParFlow uses multigrid-preconditioned Krylov and Newton–Krylov methods to solve the linear and nonlinear systems within each time step of the flow simulations. The code has demonstrated very efficient parallel solution capabilities. ParFlow has been coupled to geochemical reaction, land surface (e.g., the Common Land Model), and atmospheric models to study the interactions among the subsurface, land surface, and atmosphere systems across different spatial scales. This overview focuses on the current capabilities of the code, the core simulation engine, and the primary couplings of the subsurface model to other codes, taking a high-level perspective.

Seventh and ninth orders characteristic-wise alternative WENO finite difference schemes for hyperbolic conservation laws

In this work, the characteristic-wise alternative formulation of the seventh and ninth orders conservative weighted essentially non-oscillatory (AWENO) finite difference schemes is derived. The polynomial reconstruction procedure is applied to the conservative variables rather than the flux function of the classical WENO scheme. The numerical flux contains a low order term and high order derivative terms. The low order term can use arbitrary monotone fluxes that can enhance the resolution and reduce numerical dissipation of the fine scale structures while capturing shocks essentially non-oscillatory. The high order derivative terms are approximated by the central finite difference schemes. The improved performance in terms of accuracy, essentially non-oscillatory shock capturing and resolution for the complex shocked flow with fine scale structures in the classical one- and two-dimensional problems is demonstrated. However, the inclusion of the high order derivative terms is prone to generate Gibbs oscillations around a strong discontinuity and might result in a negative density and/or pressure. Therefore, a positivity-preserving limiter [Hu et al. J. Comput. Phys. 242 (2013)] is adopted to ensure the positive density and pressure in the shocked flows with extreme conditions, such as Mach 2000 jet flow problem.

Energy-stable predictor–corrector schemes for the Cahn–Hilliard equation

In this paper, we construct a new class of predictor–corrector time-stepping schemes for the Cahn–Hilliard equation, which are linear, second-order accurate in time, unconditionally energy stable, and uniquely solvable. Then, we present the stability and error estimates of the semi-discrete numerical schemes for solving the Cahn–Hilliard equation with general nonlinear bulk potentials. The semi-discrete scheme is further discretized using the compact central finite difference method. Several numerical examples are shown to verify the theoretical results. In particular, the numerical simulations show that the predictor–corrector schemes reach the second-order convergence rate at relatively larger time-step sizes than the classical linear schemes. The numerical strategies and theoretical tools developed in this article could be readily applied to study other phase-field models or models that can be cast as gradient flow problems.

Penyelesaian Numerik Masalah Syarat Batas Robin pada Persamaan Diferensial Cauchy-Euler

This research studied how the numerical solution of the Cauchy-Euler differential equation with Robin boundary conditions. There were several numerical methods that can be used to get the numerical solution of a boundary value problem, namely the finite-difference method, the shooting method, the collocation method, and others. In this study, the numerical solution of Robin's boundary condition problem was obtained by the center finite-difference and the shooting methods. From the two methods, the numerical error was compared to the exact solution. The simulation results shown that the shooting method produces a better numerical solution for approximating the completion of the Cauchy-Euler differential equation than the finite-difference method since it produced smaller numerical errors.

Computational Analysis of Natural Convection Heat Transfer on a Vertical Plate With Discrete Heat Sources

In the present study, buoyancy driven free convection flow along the vertical plate with discrete heat sources is analyzed. To illustrate free convection heat transfer along a vertical plate with finite discrete heat sources a 2-D steady-state model is considered. The thermos-physical properties of fluid are assumed constant except for the buoyancy terms, which are computed using the Boussinesq approximation of Navier-stokes equation. The two-dimensional Navier-stokes equations are solved using SIMPLER algorithm. The dimensionless equations are descretised and solved by using central difference finite difference approach. The development of the air flow caused by buoyancy induced free convective heat transfer has been studied through the progression of velocity and temperature fields. The results are obtained for various Grashof numbers Gr=103 , 104 and 5x 104 , and the influence of Grashof number on flow field has been studied. Average Nusselt number at the plate is also obtained. The effect of variation of Prandtl number at a given Grashof number is also studied.

Nonoverlapping Localized Exponential Time Differencing Methods for Diffusion Problems

In this paper, we propose nonoverlapping localized exponential time differencing (ETD) methods for diffusion problems. The model time-dependent diffusion equation is first reformulated on subdomains based on the nonoverlapping domain decomposition, in which Neumann boundary conditions are imposed on the interfaces for the subdomain problems and Dirichlet type conditions are enforced to form a space-time interface problem. After spatial discretization by standard central finite differences and temporal integration with the first or second order ETD methods, the fully discrete interface problem is obtained. Such an interface problem is then solved iteratively either at each time step or over the whole time interval: the former involves the solution of stationary problems in each subdomain at each iteration while the latter involves the solution of time-dependent subdomain problems at each iteration. For both approaches, we prove that localized ETD solutions conserve mass exactly and converge in time to the exact space semidiscrete solution. Numerical experiments in two dimensions are also presented to illustrate the performance of the proposed methods.

Operator-compensation methods with mass and energy conservation for solving the Gross-Pitaevskii equation

In this work, operator-compensation based high order methods are presented to solve the Gross-Pitaevskii equation modeling Bose-Einstein condensation (BEC). We begin with the high-order approximation to the Laplacian by the central finite difference method combined with operator-compensation technique, and the Crank-Nicolson and time-splitting methods are then proposed for the time discretization. We show that the Crank-Nicolson operator-compensation method keeps the mass and energy conservation, while the time-splitting operator-compensation method only conserves the mass. Then by extending the high-order operator-compensation methods for the first order derivative, the time-splitting finite difference method is developed to solve the Gross-Pitaevskii equation with an angular momentum rotation term. Numerical results are given to test accuracy and verify the conservation properties.

Central finite-difference of numerical solution for three-dimensional atmospheric transport equation

The energy transport equation is fundamental for the meteorological analysis; in this work, we analyze this equation in three dimensions using the methods of central finite differences; the analysis of convergence, consistency, and stability of the scheme shows a strong dependence of space and temporal variables. In conclusion, with the central finite differences was possible to predict the three-dimensional dynamics of the temperature.

Effects of compressibility and Atwood number on the single-mode Rayleigh-Taylor instability

In order to study the effect of compressibility on Rayleigh-Taylor (RT) instability, we numerically simulated the late-time evolution of two-dimensional single-mode RT instability for isothermal background stratification with different isothermal Mach numbers and Atwood numbers (At) using a high-order central compact finite difference scheme. It is found that the initial density stratification caused by compressibility plays a stabilizing role, while the expansion-compression effect of flow plays a destabilizing role. For the case of small Atwood number, the density difference between the two sides of the interface is small, and the density distribution of the upper and lower layers is nearly symmetrical. The initial density stratification plays a dominant role, and the expansion-compression effect has little influence. With the increase in the Atwood number, the stabilization effect of initial density stratification decreases, and the instability caused by the expansion-compression effect becomes more significant. The flow structures of bubbles and spikes are quite different at medium Atwood number. The effect of compressibility on the bubble velocity is strong at large At. The bubble height is approximately a quadratic function of time at potential flow growth stage. The average bubble acceleration is nearly proportional to the square of Mach number at At = 0.9.

New exact and numerical solutions with their stability for Ito integro-differential equation via Riccati–Bernoulli sub-ODE method

This paper points out several exact travelling wave solutions for -dimensional Ito integro-differential equation via the Riccati–Bernoulli sub-ODE approach. We also aim to develop a numerical solution of the respective equation using central finite difference formulas. The stability of the presented exact and numerical solutions is also deduced using the Hamiltonian system and Von Neumann's concept, respectively. Moreover, the numerical schemes are studied in terms of their accuracy. The relative error arising from executing the numerical method is exhibited. We compare our results with others published in some articles. The accomplished numerical solutions are successfully compared with the analytical ones. The used processes can be extended to solve more integrable problems as well as non-integrable ones.