High order implicit method for odes stiff systems

This paper presents a new difference scheme for numerical solution of stiff systems of ODE’s. The present study is mainly motivated to develop an absolutely stable numerical method with a high order of approximation. In this work a double implicit A-stable difference scheme with the sixth order of approximation is suggested. Another purpose of this study is to introduce automatic choice of the integration step size of the difference scheme which is derived from the proposed scheme and the one step scheme of the fourth order of approximation. The algorithm was tested by means of solving the Kreiss problem and a chemical kinetics problem. The behavior of the gas explosive mixture (H ₂ + O₂) in a closed space with a mobile piston is considered in test problem 2. It is our conclusion that a hydrogen-operated engine will permit to decrease the emitted levels of hazardous atmospheric pollutants.     

Greedy sparse linear approximations of functionals from nodal data

In many numerical algorithms, integrals or derivatives of functions have to be approximated by linear combinations of function values at nodes. This ranges from numerical integration to meshless methods for solving partial differential equations. The approximations should use as few nodal values as possible and at the same time have a smallest possible error. For each fixed set of nodes and each fixed Hilbert space of functions with continuous point evaluation, e.g. a fixed Sobolev space, there is an error–optimal method available using the reproducing kernel of the space. But the choice of the nodes is usually left open. This paper shows how to select good nodes adaptively by a computationally cheap greedy method, keeping the error optimal in the above sense for each incremental step of the node selection. This is applied to interpolation, numerical integration, and numerical differentiation. The latter case is particularly important for the design of meshless methods with sparse generalized stiffness matrices. The greedy algorithm is described in detail, and numerical examples are provided. In contrast to the usual practice, the greedy method does not always use nearest neighbors for local approximations of function values and derivatives. Furthermore, it avoids multiple points from clusters and it is better conditioned than choosing nearest neighbors.     

Step size control in the numerical solution of stochastic differential equations

We introduce a variable step size algorithm for the pathwise numerical approximation of solutions to stochastic ordinary differential equations. The algorithm is based on a new pair of embedded explicit Runge-Kutta methods of strong order 1.5(1.0), where the method of strong order 1.5 advances the numerical computation and the difference between approximations defined by the two methods is used for control of the local error. We show that convergence of our method is preserved though the discretization times are not stopping times any more, and further, we present numerical results which demonstrate the effectiveness of the variable step size implementation compared to a fixed step size implementation.     

Adaptive partitioning techniques for ordinary differential equations

Many stiff systems of ordinary differential equations (ODEs) modeling practical problems can be partitioned into loosely coupled subsystems. In this paper the objective of the partitioning is to permit the numerical integration of one time step to be performed as the solution of a sequence of small subproblems. This reduces the computational complexity compared to solving one large system and permits efficient parallel execution under appropriate conditions. The subsystems are integrated using methods based on low order backward differentiation formulas.This paper presents an adaptive partitioning algorithm based on a classical graph algorithm and techniques for the efficient evaluation of the error introduced by the partitioning.The power of the adaptive partitioning algorithm is demonstrated by a real world example, a variable step-size integration algorithm which solves a system of ODEs originating from chemical reaction kinetics. The computational savings are substantial. In memory of Germund Dahlquist (1925–2005).AMS subject classification (2000) 65L06, 65Y05     

On finite element methods for convection dominated phenomena

The aim of this work is to give theoretical justification of several types of finite element approximations to the initial-boundary value problems of first order linear hyperbolic equations. Our approximate scheme is obtained by the piecewise linear continuous finite element method for space variable, x, and the Euler type step by step integration method for time variable, t. An artificial viscosity technique, up-stream type methods are considered within the frame work of L²-theory. The convergence and the error estimate of the approximate solutions to the true one are discussed.     

A new sequential approach for solving the integro-differential equation via Haar wavelet bases

In this work, we present a method for numerical approximation of fixed point operator, particularly for the mixed Volterra–Fredholm integro-differential equations. The main tool for error analysis is the Banach fixed point theorem. The advantage of this method is that it does not use numerical integration, we use the properties of rationalized Haar wavelets for approximate of integral. The cost of our algorithm increases accuracy and reduces the calculation, considerably. Some examples are provided toillustrate its high accuracy and numerical results are compared with other methods in the other papers.     

On the convergence of the Bermúdez-Moreno algorithm with constant parameters

Summary. Bermúdez-Moreno [5] presents a duality numerical algorithm for solving variational inequalities of the second kind. The performance of this algorithm strongly depends on the choice of two constant parameters. Assuming a further hypothesis of the inf-sup type, we present here a convergence theorem that improves on the one presented in [5]: we prove that the convergence is linear, and we give the expression of the asymptotic error constant and the explicit form of the optimal parameters, as a function of some constants related to the variational inequality. Finally, we present some numerical examples that confirm the theoretical results. Received June 28, 1999 / Revised version received February 19, 2001 / Published online October 17, 2001     

On the convergence of the sequences of Gerschgorin-like disks

Carstensen’s results from 1991, connected with Gerschgorin’s disks, are used to establish a theorem concerning the localization of polynomial zeros and to derive an a posteriori error bound method. The presented quasi-interval method possesses useful property of inclusion methods to produce disks containing all simple zeros of a polynomial. The centers of these disks behave as approximations generated by a cubic derivative free method where the use of quantities already calculated in the previous iterative step decreases the computational cost. We state initial convergence conditions that guarantee the convergence of error bound method and prove that the method has the order of convergence three. Initial conditions are computationally verifiable since they depend only on the polynomial coefficients, its degree and initial approximations. Some computational aspects and the possibility of implementation on parallel computers are considered, including two numerical examples.In honor of Professor Richard S. Varga.     

Computational aspects of continuous–discrete extended Kalman-filtering

This paper elaborates how the time update of the continuous–discrete extended Kalman-filter (EKF) can be computed in the most efficient way. The specific structure of the EKF-moment differential equations leads to a hybrid integration algorithm, featuring a new Taylor–Heun-approximation of the nonlinear vector field and a modified Gauss–Legendre-scheme, generating positive semidefinite solutions for the state error covariance. Furthermore, the order of consistency and stability behavior of the outlined procedure is investigated. The results are incorporated into an algorithm with adaptive controlled step size, assuring a fixed numerical precision with minimal computational effort.     

Computational inelasticity for loading conditions on multiple time scales by adaptive step size control

This contribution proposes an algorithm based on adaptive step size control for the simulation of inelastic solids and structures undergoing loading conditions at multiple time scales. Adaptivity in time integration of viscoelastic constitutive laws is directed by an refinement indicator which is constructed from integrators of different order, here a fourth-order Runge-Kutta (RK) method and linear Backward-Euler. The key novel aspect is that by virtue of an recently established consistency condition the higher order methods, p ≥ 2, can achieve their full nominal order without fulfilling the weak form of balance of linear momentum in the RK stages, but only at the end of the time interval. A representative numerical example illustrates the performance of the present adaptive method and underpins the computational savings compared with uniform time step sizes. (© 2017 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Approximate preservation of quadratic first integrals by explicit Runge–Kutta methods

The approximate preservation of quadratic first integrals (QFIs) of differential systems in the numerical integration with Runge–Kutta (RK) methods is studied. Conditions on the coefficients of the RK method to preserve all QFIs up to a given order are obtained, showing that the pseudo-symplectic methods studied by Aubry and Chartier (BIT 98(3):439–461, 1998) of algebraic order p preserve QFIs with order q = 2p. An expression of the error of conservation of QFIs by a RK method is given, and a new explicit six-stage formula with classical order four and seventh order of QFI-conservation is obtained by choosing their coefficients so that they minimize both local truncation and conservation errors. Several formulas with algebraic orders 3 and 4 and different orders of conservation have been tested with some problems with quadratic and general first integrals. It is shown that the new fourth-order explicit method preserves much better the qualitative properties of the flow than the standard fourth-order RK method at the price of two extra function evaluations per step and it is a practical and efficient alternative to the fully implicit methods required for a complete preservation of QFIs.     

Two‐level discretization of the Navier‐Stokes equations with r‐Laplacian subgridscale viscosity

The r‐Laplacian has played an important role in the development of computationally efficient models for applications, such as numerical simulation of turbulent flows. In this article, we examine two‐level finite element approximation schemes applied to the Navier‐Stokes equations with r‐Laplacian subgridscale viscosity, where r is the order of the power‐law artificial viscosity term. In the two‐level algorithm, the solution to the fully nonlinear coarse mesh problem is utilized in a single‐step linear fine mesh problem. When modeling parameters are chosen appropriately, the error in the two‐level algorithm is comparable to the error in solving the fully nonlinear problem on the fine mesh. We provide rigorous numerical analysis of the two‐level approximation scheme and derive scalings which vary based on the coefficient r, coarse mesh size H, fine mesh size h, and filter radius δ. We also investigate the two‐level algorithm in several computational settings, including the 3D numerical simulation of flow past a backward‐facing step at Reynolds number Re = 5100. In all numerical tests, the two‐level algorithm was proven to achieve the same order of accuracy as the standard one‐level algorithm, at a fraction of the computational cost. © 2011 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2011     

Solving monotone inclusions with linear multi-step methods

 In this paper a new class of proximal-like algorithms for solving monotone inclusions of the form T(x)∋0 is derived. It is obtained by applying linear multi-step methods (LMM) of numerical integration in order to solve the differential inclusion , which can be viewed as a generalization of the steepest decent method for a convex function. It is proved that under suitable conditions on the parameters of the LMM, the generated sequence converges weakly to a point in the solution set T ⁻¹ (0). The LMM is very similar to the classical proximal point algorithm in that both are based on approximately evaluating the resolvants of T. Consequently, LMM can be used to derive multi-step versions of many of the optimization methods based on the classical proximal point algorithm. The convergence analysis allows errors in the computation of the iterates, and two different error criteria are analyzed, namely, the classical scheme with summable errors, and a recently proposed more constructive criterion. Received: April 2001 / Accepted: November 2002 Published online: February 14, 2003 Key Words. proximal point algorithm – monotone operator – numerical integration – strong stability – relative error criterion Mathematics Subject Classification (1991): 20E28, 20G40, 20C20     

Adaptive stochastic numerical scheme in parallel random walk models for transport problems in shallow water

This paper deals with the simulation of transport of pollutants in shallow water using random walk models and develops several computation techniques to speed up the numerical integration of the stochastic differential equations (SDEs). This is achieved by using both random time stepping and parallel processing.We start by considering a basic stochastic Euler scheme for integration of the diffusion and drift terms of the SDEs, with a strong order 1 in the strong sense. The errors due to this scheme depend on the location of the pollutant; it is dominated by the diffusion term near boundaries, and by the deterministic drift further away from the boundaries. Using a pair of integration schemes, one of strong order 1.5 near the boundary and one of strong order 2.0 elsewhere, we can estimate the error and approximate an optimal step size for a given error tolerance. The resulting algorithm is developed such that it allows for complete flexibility of the step size, while guaranteeing the correct Brownian behaviour.Modelling pollutants by non-interacting particles enables the use of parallel processing in the simulation. We take advantage of this by implementing the algorithm using the MPI library. The inherent asynchronic nature of the particle simulation, in addition to the parallel processing, makes it difficult to get a coherent picture of the results at any given points. However, by inserting internal synchronisation points in the temporal discretisation, the code allows pollution snapshots and particle counts to be made at times specified by the user.     

Stabilized mixed finite element methods for the Navier‐Stokes equations with damping

In this paper, the stabilized mixed finite element methods are presented for the Navier‐Stokes equations with damping. The existence and uniqueness of the weak solutions are proven by use of the Brouwer fixed‐point theorem. Then, optimal error estimates for the H¹‐norm and L²‐norm of the velocity and the L²‐norm of the pressure are derived. Moreover, on the basis of the optimal L²‐norm error estimate of the velocity, a stabilized two‐step method is proposed, which is more efficient than the usual stabilized methods. Finally, two numerical examples are implemented to confirm the theoretical analysis.     

High‐order numerical solution of second‐order one‐dimensional hyperbolic telegraph equation using a shifted Gegenbauer pseudospectral method

We present a high‐order shifted Gegenbauer pseudospectral method (SGPM) to solve numerically the second‐order one‐dimensional hyperbolic telegraph equation provided with some initial and Dirichlet boundary conditions. The framework of the numerical scheme involves the recast of the problem into its integral formulation followed by its discretization into a system of well‐conditioned linear algebraic equations. The integral operators are numerically approximated using some novel shifted Gegenbauer operational matrices of integration. We derive the error formula of the associated numerical quadratures. We also present a method to optimize the constructed operational matrix of integration by minimizing the associated quadrature error in some optimality sense. We study the error bounds and convergence of the optimal shifted Gegenbauer operational matrix of integration. Moreover, we construct the relation between the operational matrices of integration of the shifted Gegenbauer polynomials and standard Gegenbauer polynomials. We derive the global collocation matrix of the SGPM, and construct an efficient computational algorithm for the solution of the collocation equations. We present a study on the computational cost of the developed computational algorithm, and a rigorous convergence and error analysis of the introduced method. Four numerical test examples have been carried out to verify the effectiveness, the accuracy, and the exponential convergence of the method. The SGPM is a robust technique, which can be extended to solve a wide range of problems arising in numerous applications. © 2015 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq 32: 307–349, 2016     

On an algorithm of Milne and Reynolds

The Milne-Reynolds averaging technique is extended to all weakly stable methods of numerical integration of ordinary differential equations, and a numerical example is presented. Also, a Milne-Reynolds average is given which reduces by a factor ofO(h ^2l+1) the unstable component of the error arising with Milne's methods without changing the order of the truncation error. The average is given explicitly forl=1,2.     

One-leg variable-coefficient formulas for ordinary differential equations and local–global step size control

In this paper we discuss a class of numerical algorithms termed one-leg methods. This concept was introduced by Dahlquist in 1975 with the purpose of studying nonlinear stability properties of multistep methods for ordinary differential equations. Later, it was found out that these methods are themselves suitable for numerical integration because of good stability. Here, we investigate one-leg formulas on nonuniform grids. We prove that there exist zero-stable one-leg variable-coefficient methods at least up to order 11 and give examples of two-step methods of orders 2 and 3. In this paper we also develop local and global error estimation techniques for one-leg methods and implement them with the local–global step size selection suggested by Kulikov and Shindin in 1999. The goal of this error control is to obtain automatically numerical solutions for any reasonable accuracy set by the user. We show that the error control is more complicated in one-leg methods, especially when applied to stiff problems. Thus, we adapt our local–global step size selection strategy to one-leg methods.     

On one class of numerical methods for solving stiff systems of differential equations

The paper is concerned with a modification of a Cowell-type method: at each step the solution is evaluated at several points and only some first of these values are retained. Stability of these methods is examined. In particular, among the methods of this group we single out one method that has the fourth order of accuracy and is stable for stiff systems with any choice of integration step.     

Cubature formulas for function spaces with moderate smoothness

We construct simple algorithms for high-dimensional numerical integration of function classes with moderate smoothness. These classes consist of square-integrable functions over the d-dimensional unit cube whose coefficients with respect to certain multiwavelet expansions decay rapidly. Such a class contains discontinuous functions on the one hand and, for the right choice of parameters, the quite natural d-fold tensor product of a Sobolev space H^s[0,1] on the other hand.The algorithms are based on one-dimensional quadrature rules appropriate for the integration of the particular wavelets under consideration and on Smolyak's construction. We provide upper bounds for the worst-case error of our cubature rule in terms of the number of function calls. We additionally prove lower bounds showing that our method is optimal in dimension d=1 and almost optimal (up to logarithmic factors) in higher dimensions. We perform numerical tests which allow the comparison with other cubature methods.