Hybrid BDF methods for the numerical solutions of ordinary differential equations

In this article, we have presented the details of hybrid methods which are based on backward differentiation formula (BDF) for the numerical solutions of ordinary differential equations (ODEs). In these hybrid BDF, one additional stage point (or off-step point) has been used in the first derivative of the solution to improve the absolute stability regions. Stability domains of our presented methods have been obtained showing that all these new methods, we say HBDF, of order p, p = 2,4,..., 12, are A(α)-stable whereas they have wide stability regions comparing with those of some known methods like BDF, extended BDF (EBDF), modified EBDF (MEBDF), adaptive EBDF (A-EBDF), and second derivtive Enright methods. Numerical results are also given for five test problems.     

Exact optimal values of step-size coefficients for boundedness of linear multistep methods

Linear multistep methods (LMMs) applied to approximate the solution of initial value problems—typically arising from method-of-lines semidiscretizations of partial differential equations—are often required to have certain monotonicity or boundedness properties (e.g., strong-stability-preserving, total-variation-diminishing or total-variation-boundedness properties). These properties can be guaranteed by imposing step-size restrictions on the methods. To qualitatively describe the step-size restrictions, one introduces the concept of step-size coefficient for monotonicity (SCM, also referred to as the strong-stability-preserving (SSP) coefficient) or its generalization, the step-size coefficient for boundedness (SCB). An LMM with larger SCM or SCB is more efficient, and the computation of the maximum SCM for a particular LMM is now straightforward. However, it is more challenging to decide whether a positive SCB exists, or determine if a given positive number is a SCB. Theorems involving sign conditions on certain linear recursions associated to the LMM have been proposed in the literature that allow us to answer the above questions: the difficulty with these theorems is that there are in general infinitely many sign conditions to be verified. In this work, we present methods to rigorously check the sign conditions. As an illustration, we confirm some recent numerical investigations concerning the existence of positive SCBs in the BDF and in the extrapolated BDF (EBDF) families. As a stronger result, we determine the optimal values of the SCBs as exact algebraic numbers in the BDF family (with 1 ≤ k ≤ 6 steps) and in the Adams–Bashforth family (with 1 ≤ k ≤ 3 steps).     

Stability of fully discrete schemes with interpolation-type fractional formulas for distributed-order subdiffusion equations

Two fully discrete methods are investigated for simulating the distributed-order sub-diffusion equation in Caputo’s form. The fractional Caputo derivative is approximated by the Caputo’s BDF1 (called L1 early) and BDF2 (or L1-2 when it was first introduced) approximations, which are constructed by piecewise linear and quadratic interpolating polynomials, respectively. It is shown that the first scheme, using the BDF1 formula, possesses the discrete minimum-maximum principle and nonnegativity preservation property such that it is stable and convergent in the maximum norm. The method using the BDF2 formula is shown to be stable and convergent in the discrete H ¹ norm by using the discrete energy method. For problems of distributed order within a certain region, the method is also proven to preserve the discrete maximum principle and nonnegativity property. Extensive numerical experiments are provided to show the effectiveness of numerical schemes, and to examine the initial singularity of the solution. The applicability of our numerical algorithms to a problem with solution which lacks the smoothness near the initial time is examined by employing a class of power-type nonuniform meshes.     

Diagonalizable Extended Backward Differentiation Formulas

We generalize the extended backward differentiation formulas (EBDFs) introduced by Cash and by Psihoyios and Cash so that the system matrix in the modified Newton process can be block-diagonalized, enabling an efficient parallel implementation. The purpose of this paper is to justify the use of diagonalizable EBDFs on parallel computers and to offer a starting point for the development of a variable stepsize-variable order method. We construct methods which are L-stable up to order p = 6 and which have the same computational complexity per processor as the conventional BDF methods. Numerical experiments with the order 6 method show that a speedup factor of between 2 and 4 on four processors can be expected.     

Convergence of parallel multistep hybrid methods for singular perturbation problems

Parallel multistep hybrid methods (PHMs) can be implemented in parallel with two processors, accordingly have almost the same computational speed per integration step as BDF methods of the same order with the same stepsize. But PHMs have better stability properties than BDF methods of the same order for stiff differential equations. In the present paper, we give some results on error analysis of A(α)-stable PHMs for the initial value problems of ordinary differential equations in singular perturbation form. Our convergence results are similar to those of linear multistep methods (such as BDF methods), i.e. the convergence orders are equal to their classical convergence orders, and no order reduction occurs. Some numerical examples also confirm our results.     

Numerical solution of a thermomechanical coupled model governing glacier evolution

In this paper the numerical solution of a highly nonlinear model for the thermomechanical behavior of polythermal glaciers is presented. The modeling follows the shallow ice approximation (SIA) for glaciers introduced in Fowler (1997) [13]. The model has been extended to incorporate additional moving boundaries and other nonlinear features. Moreover, a fixed domain formulation is proposed to avoid the computational drawbacks of a time-dependent domain in the numerical simulation with front tracking methods. In this setting, the coupled problem is decomposed into different nonlinear problems which allow one to obtain sequentially the profile evolution, the velocity field, the glacier surface and atmospheric temperatures, basal magnitudes and the temperature distribution inside the ice mass. A fixed point iteration algorithm converges to the solution of the nonlinear coupled problem. Among different numerical methods involved in the solution of the subproblems, characteristic schemes for time discretization, finite elements for spatial discretization, duality methods for the nonlinearities associated to maximal monotone operators and a Newton scheme for the nonlinear viscous term are proposed. Several numerical simulation examples illustrate the performance of the numerical methods and the behavior of the involved physical magnitudes.     

New class of hybrid BDF methods for the computation of numerical solutions of IVPs

A new class of hybrid BDF-like methods is presented for solving systems of ordinary differential equations (ODEs) by using the second derivative of the solution in the stage equation of class 2 + 1hybrid BDF-like methods to improve the order and stability regions of these methods. An off-step point, together with two step points, has been used in the first derivative of the solution, and the stability domains of the new methods have been obtained by showing that these methods are A-stable for order p, p =?3,4,5,6,7and A(α)-stable for order p, 8 ≤ p ≤?14. The numerical results are also given for four test problems by using variable and fixed step-size implementations.     

Parallel iterative linear solvers for multistep Runge-Kutta methods

This paper deals with solving stiff systems of differential equations by implicit Multistep Runge-Kutta (MRK) methods. For this type of methods, nonlinear systems of dimension sd arise, where s is the number of Runge-Kutta stages and d the dimension of the problem. Applying a Newton process leads to linear systems of the same dimension, which can be very expensive to solve in practice. With a parallel iterative linear system solver, especially designed for MRK methods, we approximate these linear systems by s systems of dimension d, which can be solved in parallel on a computer with s processors. In terms of Jacobian evaluations and LU-decompositions, the k-steps-stage MRK applied with this technique is on s processors equally expensive as the widely used k-step Backward Differentiation Formula on 1 processor, whereas the stability properties are better than that of BDF. A simple implementation of both methods shows that, for the same number of Newton iterations, the accuracy delivered by the new method is higher than that of BDF.     

Parallel iteration of the extended backward differentiation formulas

The extended backward differentiation formulas (EBDFs) and theirmodified form (MEBDF) were proposed by Cash in the 1980s forsolving initial value problems (IVPs) for stiff systems of ordinarydifferential equations (ODEs). In a recent performance evaluationof various IVP solvers, including a variable-step-variable-orderimplementation of the MEBDF method by Cash, it turned out thatthe MEBDF code often performs more efficiently than codes likeRADAU5, DASSL and VODE. This motivated us to look at possibleparallel implementations of the MEBDF method. Each MEBDF stepessentially consists of successively solving three non-linearsystems by means of modified Newton iteration using the sameJacobian matrix. In a direct implementation of the MEBDF methodon a parallel computer system, the only scope for (coarse grain)parallelism consists of a number of parallel vector updates.However, all forward–backward substitutions and all right-hand-sideevaluations have to be done in sequence. In this paper, ourstarting point is the original (unmodified) EBDF method. Asa consequence, two different Jacobian matrices are involvedin the modified Newton method, but on a parallel computer system,the effective Jacobian-evaluation and the LU decomposition costsare not increased. Furthermore, we consider the simultaneoussolution, rather than the successive solution, of the threenon-linear systems, so that in each iteration the forward–backwardsubstitutions and the right-hand-side evaluations can be doneconcurrently. A mutual comparison of the performance of theparallel EBDF approach and the MEBDF approach shows that wecan expect a speed-up factor of about 2 on three processors.     

Fast equal and biased distance fields for medial axis transform with meshing in mind

A method for robust and efficient medial axis transform (MAT) of arbitrary domains using distance solutions (or level sets) is presented. The distance field, d, is calculated by solving the hyperbolic-natured eikonal equation. The solution is obtained on Cartesian grids. Both the fast-marching method (FMM) and fast-sweeping method (FSM) are used to calculate d. Medial axis point clouds are then extracted based on the distance solution via a simple criterion: the Laplacian or the Hessian determinant of d(x). These point clouds in the pixel/voxel space are further thinned to single pixel wide so that medial axis curves or surfaces can be connected and splined. As an alternative to other methods, the current d-MAT procedure bypasses difficulties that are usually encountered by pure geometric methods (e.g. the Voronoi approach), especially in three dimensions, and provides better accuracy than pure thinning methods. It is also shown that the d-MAT approach provides the potential to sculpt/control the MAT form for specialized solution purposes. Various examples are given to demonstrate the current approach.     

On explicit two-derivative Runge-Kutta methods

The theory of Runge-Kutta methods for problems of the form y′?=?f(y) is extended to include the second derivative y′′?=?g(y):?=?f′(y)f(y). We present an approach to the order conditions based on Butcher’s algebraic theory of trees (Butcher, Math Comp 26:79–106, 1972), and derive methods that take advantage of cheap computations of the second derivatives. Only explicit methods are considered here where attention is given to the construction of methods that involve one evaluation of f and many evaluations of g per step. Methods with stages up to five and of order up to seven including some embedded pairs are presented. The first part of the paper discusses a theoretical formulation used for the derivation of these methods which are also of wider applicability. The second part presents experimental results for non-stiff and mildly stiff problems. The methods include those with the computation of one second derivative (plus many first derivatives) per step, and embedded methods for changing stepsize as well as those involving one first derivative (plus many second derivatives) per step. The experiments have been performed on standard problems and comparisons made with some standard explicit Runge-Kutta methods.     

A class of A(α)-stable numerical methods for stiff problems in ordinary differential equations

The A(α)-stable numerical methods (ANMs) for the number of steps k ≤ 7 for stiff initial value problems (IVPs) in ordinary differential equations (ODEs) are proposed. The discrete schemes proposed from their equivalent continuous schemes are obtained. The scaled time variable t in a continuous method, which determines the discrete coefficients of the discrete method, is chosen in such a way as to ensure that the discrete scheme attains a high order and A(α)-stability. We select the value of α for which the schemes proposed are absolutely stable. The new algorithms are found to have a comparable accuracy with that of the backward differentiation formula (BDF) discussed in [12] which implements the Ode15s in the Matlab suite.     

Controllability and the linear equation Ax=b

The linear equation Ax=b, with An×n, is considered, and it is shown that the controllable subspace of (A, b) can intersect the solution set of this equation in at most one point. Some additional properties of the controllable subspace related to the solution set are presented, and the role of controllability in establishing nonsingularity of A is exhibited.     

Efficient parallel multivalue hybrid methods for stiff differential equations

A class of efficient parallel multivalue hybrid methods for stiff differential equations are presented, which are all extremely stable at infinity,A-stable for orders 1–3 and A(α)-stable for orders 4–8. Each method of the class can be performed parallelly using two processors with each processor having almost the same computational amount per integration step as a backward differentiation formula (BDF) of the same order with the same stepsize performed in serial, whereas the former has not only much better stability properties but also a computational accuracy higher than the corresponding BDF. Theoretical analysis and numerical experiments show that the methods constructed in this paper are superior in many respects not only to BDFs but also to some other commonly used methods.     

Entropy-stable and entropy-dissipative approximations of a fourth-order quantum diffusion equation

Structure-preserving numerical schemes for a nonlinear parabolic fourth-order equation, modeling the electron transport in quantum semiconductors, with periodic boundary conditions are analyzed. First, a two-step backward differentiation formula (BDF) semi-discretization in time is investigated. The scheme preserves the nonnegativity of the solution, is entropy stable and dissipates a modified entropy functional. The existence of a weak semi-discrete solution and, in a particular case, its temporal second-order convergence to the continuous solution is proved. The proofs employ an algebraic relation which implies the G-stability of the two-step BDF. Second, an implicit Euler and $q$ -step BDF discrete variational derivative method are considered. This scheme, which exploits the variational structure of the equation, dissipates the discrete Fisher information (or energy). Numerical experiments show that the discrete (relative) entropies and Fisher information decay even exponentially fast to zero.     

On the stability at all times of linearly extrapolated BDF2 timestepping for multiphysics incompressible flow problems

We prove long‐time stability of linearly extrapolated BDF2 (BDF2LE) timestepping methods, together with finite element spatial discretizations, for incompressible Navier‐Stokes equations (NSE) and related multiphysics problems. For the NSE, Boussinesq, and magnetohydrodynamics schemes, we prove unconditional long time L² stability, provided external forces (and sources) are uniformly bounded in time. We also provide numerical experiments to compare stability of BDF2LE to linearly extrapolated Crank‐Nicolson scheme for NSE, and find that BDF2LE has better stability properties, particularly for smaller viscosity values. © 2016 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 33: 999–1017, 2017     

Algorithms for solutions of extended general mixed variational inequalities and fixed points

In this paper, we introduce and study a new class of extended general nonlinear mixed variational inequalities and a new class of extended general resolvent equations and establish the equivalence between the extended general nonlinear mixed variational inequalities and implicit fixed point problems as well as the extended general resolvent equations. Then by using this equivalent formulation, we discuss the existence and uniqueness of solution of the problem of extended general nonlinear mixed variational inequalities. Applying the aforesaid equivalent alternative formulation and a nearly uniformly Lipschitzian mapping S, we construct some new resolvent iterative algorithms for finding an element of set of the fixed points of nearly uniformly Lipschitzian mapping S which is the unique solution of the problem of extended general nonlinear mixed variational inequalities. We study convergence analysis of the suggested iterative schemes under some suitable conditions. We also suggest and analyze a class of extended general resolvent dynamical systems associated with the extended general nonlinear mixed variational inequalities and show that the trajectory of the solution of the extended general resolvent dynamical system converges globally exponentially to the unique solution of the extended general nonlinear mixed variational inequalities. The results presented in this paper extend and improve some known results in the literature.     

Using constraint preconditioners with regularized saddle-point problems

The problem of finding good preconditioners for the numerical solution of a certain important class of indefinite linear systems is considered. These systems are of a 2 by 2 block (KKT) structure in which the (2,2) block (denoted by -C) is assumed to be nonzero. In Constraint preconditioning for indefinite linear systems, SIAM J. Matrix Anal. Appl. 21 (2000), Keller, Gould and Wathen introduced the idea of using constraint preconditioners that have a specific 2 by 2 block structure for the case of C being zero. We shall give results concerning the spectrum and form of the eigenvectors when a preconditioner of the form considered by Keller, Gould and Wathen is used but the system we wish to solve may have C ≠0. In particular, the results presented here indicate clustering of eigenvalues and, hence, faster convergence of Krylov subspace iterative methods when the entries of C are small; such a situations arise naturally in interior point methods for optimization and we present results for such problems which validate our conclusions.     

Variational approach to some damped Dirichlet problems with impulses

In this paper, we study the existence of solutions for damped nonlinear impulsive differential equations with Dirichlet boundary conditions. By using critical point theory and variational methods, we give some new criteria to guarantee that the impulsive problems have at least one solution. Some recent results are extended and significantly improved. Finally, some examples are presented to illustrate our main results. Copyright © 2013 John Wiley & Sons, Ltd.     

Approximation of weak adjoints by reverse automatic differentiation of BDF methods

We shed light on the relation between the discrete adjoints of multistep backward differentiation formula (BDF) methods and the solution of the adjoint differential equation. To this end, we develop a functional-analytic framework based on a constrained variational problem and introduce the notion of weak adjoint solutions of ordinary differential equations. We devise a Petrov-Galerkin finite element (FE) interpretation of the BDF method and its discrete adjoint scheme obtained by reverse internal numerical differentiation. We show how the FE approximation of the weak adjoint is computed by the discrete adjoint scheme and prove its convergence in the space of normalized functions of bounded variation. We also show convergence of the discrete adjoints to the classical adjoints on the inner time interval. Finally, we give numerical results for non-adaptive and fully adaptive BDF schemes. The presented framework opens the way to carry over techniques on global error estimation from FE methods to BDF methods.