Extended RKN-type methods for numerical integration of perturbed oscillators

In this paper, extended Runge-Kutta-Nyström-type methods for the numerical integration of perturbed oscillators with low frequencies are presented, which inherit the framework of RKN methods and make full use of the special feature of the true flows for both the internal stages and the updates. Following the approach of J. Butcher, E. Hairer and G. Wanner, we develop a new kind of tree set to derive order conditions for the extended Runge-Kutta-Nyström-type methods. The numerical stability and phase properties of the new methods are analyzed. Numerical experiments are accompanied to show the applicability and efficiency of our new methods in comparison with some well-known high quality methods proposed in the scientific literature.     

A Runge-Kutta-Nyström pair for the numerical integration of perturbed oscillators

New Runge-Kutta-Nyström methods especially designed for the numerical integration of perturbed oscillators are presented in this paper. They are capable of exactly integrating the harmonic or unperturbed oscillator. We construct an embedded 4(3) RKN pair that is based on the FSAL property. The new method is much more efficient than previously derived RKN methods for some subclasses of problems.     

Computation of the eigenvalues of the Schrödinger equation by exponentially-fitted Runge-Kutta-Nyström methods

In this work we consider exponentially fitted and trigonometrically fitted Runge-Kutta-Nyström methods. These methods integrate exactly differential systems whose solutions can be expressed as linear combinations of the set of functions exp(wx), exp(−wx), or sin(wx), cos(wx), w∈ℜ. We modify existing RKN methods of fifth and sixth order. We apply these methods to the computation of the eigenvalues of the Schrödinger equation with different potentials as the harmonic oscillator, the doubly anharmonic oscillator and the exponential potential.     

Trigonometrically fitted explicit Numerov-type method for periodic IVPs with two frequencies

In this paper, new trigonometrically fitted Numerov type methods for the periodic initial problems are proposed. These methods are based on the original Numerov-type sixth order method with fifth internal stages motivated by Tsitouras (see [Ch. Tsitouras, Explicit Numerov type methods with reduced number of stages, Comput. Math. Appl. 45 (2003) 37-42]). Some parameters are added to these methods so that they can integrate exactly the combination of trigonometrically functions with two frequencies. Numerical stability and phase properties of the new methods are analyzed. Numerical experiments are carried out to show the efficiency and robustness of our new methods in comparison with the well known codes proposed in the scientific literature.     

A fourth-order symplectic exponentially fitted integrator

A numerical method for ordinary differential equations is called symplectic if, when applied to Hamiltonian problems, it preserves the symplectic structure in phase space, thus reproducing the main qualitative property of solutions of Hamiltonian systems. In a previous paper [G. Vanden Berghe, M. Van Daele, H. Van de Vyver, Exponential fitted Runge-Kutta methods of collocation type: fixed or variable knot points?, J. Comput. Appl. Math. 159 (2003) 217-239] some exponentially fitted RK methods of collocation type are proposed. In particular, three different versions of fourth-order exponentially fitted Gauss methods are described. It is well known that classical Gauss methods are symplectic. In contrast, the exponentially fitted versions given in [G. Vanden Berghe, M. Van Daele, H. Van de Vyver, Exponential fitted Runge-Kutta methods of collocation type: fixed or variable knot points?, J. Comput. Appl. Math. 159 (2003) 217-239] do not share this property. This paper deals with the construction of a fourth-order symplectic exponentially fitted modified Gauss method. The RK method is modified in the sense that two free parameters are added to the Buthcher tableau in order to retain symplecticity.     

A numerical differentiation library exploiting parallel architectures

We present a software library for numerically estimating first and second order partial derivatives of a function by finite differencing. Various truncation schemes are offered resulting in corresponding formulas that are accurate to order O(h), O(h²), and O(h⁴), h being the differencing step. The derivatives are calculated via forward, backward and central differences. Care has been taken that only feasible points are used in the case where bound constraints are imposed on the variables. The Hessian may be approximated either from function or from gradient values. There are three versions of the software: a sequential version, an OpenMP version for shared memory architectures and an MPI version for distributed systems (clusters). The parallel versions exploit the multiprocessing capability offered by computer clusters, as well as modern multi-core systems and due to the independent character of the derivative computation, the speedup scales almost linearly with the number of available processors/cores.

Program summary

Program title: NDL (Numerical Differentiation Library)Catalogue identifier: AEDG_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEDG_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlNo. of lines in distributed program, including test data, etc.: 73 030No. of bytes in distributed program, including test data, etc.: 630 876Distribution format: tar.gzProgramming language: ANSI FORTRAN-77, ANSI C, MPI, OPENMPComputer: Distributed systems (clusters), shared memory systemsOperating system: Linux, SolarisHas the code been vectorised or parallelized?: YesRAM: The library uses O(N) internal storage, N being the dimension of the problemClassification: 4.9, 4.14, 6.5Nature of problem: The numerical estimation of derivatives at several accuracy levels is a common requirement in many computational tasks, such as optimization, solution of nonlinear systems, etc. The parallel implementation that exploits systems with multiple CPUs is very important for large scale and computationally expensive problems.Solution method: Finite differencing is used with carefully chosen step that minimizes the sum of the truncation and round-off errors. The parallel versions employ both OpenMP and MPI libraries.Restrictions: The library uses only double precision arithmetic.Unusual features: The software takes into account bound constraints, in the sense that only feasible points are used to evaluate the derivatives, and given the level of the desired accuracy, the proper formula is automatically employed.Running time: Running time depends on the function's complexity. The test run took 15 ms for the serial distribution, 0.6 s for the OpenMP and 4.2 s for the MPI parallel distribution on 2 processors.     

A Mathematica program for the two-step twelfth-order method with multi-derivative for the numerical solution of a one-dimensional Schrödinger equation

In this paper, we present the detailed Mathematica symbolic derivation and the program which is used to integrate a one-dimensional Schrödinger equation by a new two-step numerical method. We add the fourth- and sixth-order derivatives to raise the precision of the traditional Numerov's method from fourth order to twelfth order, and to expand the interval of periodicity from (0,6) to the one of (0,9.7954) and (9.94792,55.6062). In the program we use an efficient algorithm to calculate the first-order derivative and avoid unnecessarily repeated calculation resulting from the multi-derivatives. We use the well-known Woods-Saxon's potential to test our method. The numerical test shows that the new method is not only superior to the previous lower order ones in accuracy, but also in the efficiency. This program is specially applied to the problem where a high accuracy or a larger step size is required.

Program summary

Title of program: ShdEq.nbCatalogue number: ADTTProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADTTProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions: noneComputer for which the program is designed and others on which it has been tested: The program has been designed for the microcomputer and been tested on the microcomputer.Computers: IBM PCOperating systems under which the program has been tested: Windows XPProgramming language used: Mathematica 4.2Memory required to execute with typical data: 51 712 bytesNo. of bytes in distributed program, including test data, etc.: 45 381No. of lines in distributed program, including test data, etc.: 7311Distribution format: tar gzip fileCPC Program Library subprograms used: noNature of physical problem: Numerical integration of one-dimensional or radial Schrödinger equation to find the eigenvalues for a bound states and phase shift for a continuum state.Method of solution: Using a two-step method twelfth-order method to integrate a Schrödinger equation numerically from both two ends and the connecting conditions at the matching point, an eigenvalue for a bound state or a resonant state with a given phase shift can be found.Restrictions on the complexity of the problem: The analytic form of the potential function and its high-order derivatives must be known.Typical running time: Less than one second.Unusual features of the program: Take advantage of the high-order derivatives of the potential function and efficient algorithm, the program can provide all the numerical solution of a given Schrödinger equation, either a bound or a resonant state, with a very high precision and within a very short CPU time. The program can apply to a very broad range of problems because the method has a very large interval of periodicity.References: [1] T.E. Simos, Proc. Roy. Soc. London A 441 (1993) 283.[2] Z. Wang, Y. Dai, An eighth-order two-step formula for the numerical integration of the one-dimensional Schrödinger equation, Numer. Math. J. Chinese Univ. 12 (2003) 146.[3] Z. Wang, Y. Dai, An twelfth-order four-step formula for the numerical integration of the one-dimensional Schrödinger equation, Internat. J. Modern Phys. C 14 (2003) 1087.     

P-stable linear symmetric multistep methods for periodic initial-value problems

In this paper we present a new kind of P-stable multistep methods for periodic initial-value problems. From the numerical results obtained by the new method to well-known periodic problems, show the superior efficiency, accuracy, stability of the method presented in this paper.     

Linearization methods in classical and quantum mechanics

The applicability and accuracy of linearization methods for initial-value problems in ordinary differential equations are verified on examples that include the nonlinear Duffing equation, the Lane-Emden equation, and scattering length calculations. Linearization methods provide piecewise linear ordinary differential equations which can be easily integrated, and provide accurate answers even for hypersingular potentials, for which perturbation methods diverge. It is shown that the accuracy of linearization methods can be substantially improved by employing variable steps which adjust themselves to the solution.     

Trigonometrically-fitted method for a periodic initial value problem with two frequencies

A second-order differential equation whose solution is periodic with two frequencies has important applications in many scientific fields. Nevertheless, it may exhibit ‘periodic stiffness’ for most of the available linear multi-step methods. The phenomena are similar to the popular Stömer-Cowell class of linear multi-step methods for one-frequency problems. According to the stability theory laid down by Lambert, ‘periodic stiffness’ appears in a two-frequency problem because the production of the step-length and the bigger angular frequency lies outside the interval of periodicity. On the other hand, for a two-frequency problem, even with a small step-length, the error in the numerical solution afforded by a P-stable trigonometrically-fitted method with one frequency would be too large for practical applications. In this paper we demonstrate that the interval of periodicity and the local truncation error of a linear multi-step method for a two-frequency problem can be greatly improved by a new trigonometric-fitting technique. A trigonometrically-fitted Numerov method with two frequencies is proposed and has been verified to be P-stable with vanishing local truncation error for a two-frequency test problem. Numerical results demonstrated that the proposed trigonometrically-fitted Numerov method with two frequencies has significant advantages over other types of Numerov methods for solving the ‘periodic stiffness’ problem.     

Trigonometrically-fitted method with the Fourier frequency spectrum for undamped Duffing equation

With non-linearities, the frequency spectrum of an undamped Duffing oscillator should be composed of odd multiples of the driving frequency which can be interpreted as resonance driving terms. It is expected that the frequency spectrum of the corresponding numerical solution with high accurateness should contain nearly the same components. Hence, to contain these Fourier components and to calculate the amplitudes of these components in a more accurate and efficient way is the key to develop a new numerical method with high stability, accuracy and efficiency for the Duffing equation. To explore the possibility of using trigonometrically-fitting technique to build a numerical method with resonance spectrum, we design four types of Numerov methods, in which the first one is the traditional Numerov method, which contains no Fourier component, the second one contains only the first resonance term, the third one contains the first two resonance terms, and the last one contains the first three resonance terms, and apply them to the well-known undamped Duffing equation with Dooren's parameters. The numerical results demonstrate that the Numerov method fitted with the Fourier components is much more stable, accurate and efficient than the one with no Fourier component. The accuracy of the fitted method with the first three Fourier components can attain 10⁻⁹ for a remarkable range of step sizes, including nearly infinite, except individual small range of instability, which is much higher than the one of the traditional Numerov method, with eight orders for step size of π/2.011.     

A new kind of high-efficient and high-accurate P-stable Obrechkoff three-step method for periodic initial-value problems

In order to improve the efficiency and accuracy of the previous Obrechkoff method, in this paper we put forward a new kind of P-stable three-step Obrechkoff method of O(h¹⁰) for periodic initial-value problems. By using a new structure and an embedded high accurate first-order derivative formula, we can avoid time-consuming iterative calculation to obtain the high-order derivatives. By taking advantage of new trigonometrically-fitting scheme we can make both the main structure and the first-order derivative formula to be P-stable. We apply our new method to three periodic problems and compare it with the previous three Obrechkoff methods. Numerical results demonstrate that our new method is superior over the previous ones in accuracy, efficiency and stability.     

Arbitrary high order finite-volume methods for electromagnetic wave propagation

Problems in electromagnetic wave propagation often require high accuracy approximations with low resolution computational grids. For non-stationary problems such schemes should possess the same approximation order in space and time. In the present article we propose for electromagnetic applications an explicit class of robust finite-volume (FV) schemes for the Maxwell equations. To achieve high accuracy we combine the FV method with the so-called ADER approach resulting in schemes which are arbitrary high order accurate in space and time. Numerical results and convergence investigations are shown for two and three-dimensional test cases on Cartesian grids, where the used FV-ADER schemes are up to 8th order accurate in both space and time.     

Multi-dimensional semi-Lagrangian scheme that guarantees exact conservation

A new numerical method that guarantees exact mass conservation is proposed to solve multi-dimensional hyperbolic equations in semi-Lagrangian form without directional splitting. The method is based on a concept of CIP scheme and keep the many good characteristics of the original CIP scheme. The CIP strategy is applied to the integral form of variable. Although the advection and non-advection terms are separately treated, the mass conservation is kept in a form of spatial profile inside a grid cell. Therefore, it retains various advantages of the semi-Lagrangian schemes with exact conservation that has been beyond the capability of conventional semi-Lagrangian schemes.     

Piecewise quasilinearization techniques for singular boundary-value problems

Piecewise quasilinearization methods for singular boundary-value problems in second-order ordinary differential equations are presented. These methods result in linear constant-coefficients ordinary differential equations which can be integrated analytically, thus yielding piecewise analytical solutions. The accuracy of the globally smooth piecewise quasilinear method is assessed by comparisons with exact solutions of several Lane-Emden equations, a singular problem of non-Newtonian fluid dynamics and the Thomas-Fermi equation. It is shown that the smooth piecewise quasilinearization method provides accurate solutions even near the singularity and is more precise than (iterative) second-order accurate finite difference discretizations. It is also shown that the accuracy of the smooth piecewise quasilinear method depends on the kind of singularity, nonlinearity and inhomogeneities of singular ordinary differential equations. For the Thomas-Fermi equation, it is shown that the piecewise quasilinearization method that provides globally smooth solutions is more accurate than that which only insures global continuity, and more accurate than global quasilinearization techniques which do not employ local linearization.     

Towards a more consistent discretization scheme for adaptive implicit RHD computations

We present a new implicit numerical discretization for the equations of radiation hydrodynamics (RHD) which is based on a more geometrical representation of a finite volume scheme suitable for spherical systems. In particular, the motion of the grid points is directly included by appropriate volume changes. Several examples illustrate the accuracy gained by this improved difference scheme.     

Arbitrarily precise numerical solutions of the one-dimensional Schrödinger equation

In this paper, how to overcome the barrier for a finite difference method to obtain the numerical solutions of a one-dimensional Schrödinger equation defined on the infinite integration interval accurate than the computer precision is discussed. Five numerical examples of solutions with the error less than 10⁻⁵⁰ and 10⁻³⁰ for the bound and resonant state, respectively, obtained by the Obrechkoff one-step method implemented in the multi precision mode, which include the harmonic oscillator, the Pöschl-Teller potential, the Morse potential and the Woods-Saxon potential, demonstrate that the finite difference method can yield the eigenvalues of a complex potential with an arbitrarily desired precision within a reasonable efficiency.     

FiEstAS sampling—a Monte Carlo algorithm for multidimensional numerical integration

This paper describes a new algorithm for Monte Carlo integration, based on the Field Estimator for Arbitrary Spaces (FiEstAS). The algorithm is discussed in detail, and its performance is evaluated in the context of Bayesian analysis, with emphasis on multimodal distributions with strong parameter degeneracies. Source code is available upon request.     

Numerical ‘health check’ for scientific codes: the CADNA approach

Scientific computation has unavoidable approximations built into its very fabric. One important source of error that is difficult to detect and control is round-off error propagation which originates from the use of finite precision arithmetic. We propose that there is a need to perform regular numerical ‘health checks’ on scientific codes in order to detect the cancerous effect of round-off error propagation. This is particularly important in scientific codes that are built on legacy software. We advocate the use of the CADNA library as a suitable numerical screening tool. We present a case study to illustrate the practical use of CADNA in scientific codes that are of interest to the Computer Physics Communications readership. In doing so we hope to stimulate a greater awareness of round-off error propagation and present a practical means by which it can be analyzed and managed.     

Tension spline approach for the numerical solution of nonlinear Klein-Gordon equation

The nonlinear Klein-Gordon equation describes a variety of physical phenomena such as dislocations, ferroelectric and ferromagnetic domain walls, DNA dynamics, and Josephson junctions. We derive approximate expressions for the dispersion relation of the nonlinear Klein-Gordon equation in the case of strong nonlinearities using a method based on the tension spline function and finite difference approximations. The resulting spline difference schemes are analyzed for local truncation error, stability and convergence. It has been shown that by suitably choosing the parameters, we can obtain two schemes of O(k²+k²h²+h²) and O(k²+k²h²+h⁴). In the end, some numerical examples are provided to demonstrate the effectiveness of the proposed schemes.