An efficient swap algorithm for the lattice Boltzmann method

During the last decade, the lattice-Boltzmann method (LBM) as a valuable tool in computational fluid dynamics has been increasingly acknowledged. The widespread application of LBM is partly due to the simplicity of its coding. The most well-known algorithms for the implementation of the standard lattice-Boltzmann equation (LBE) are the two-lattice and two-step algorithms. However, implementations of the two-lattice or the two-step algorithm suffer from high memory consumption or poor computational performance, respectively. Ultimately, the computing resources available decide which of the two disadvantages is more critical. Here we introduce a new algorithm, called the swap algorithm, for the implementation of LBE. Simulation results demonstrate that implementations based on the swap algorithm can achieve high computational performance and have very low memory consumption. Furthermore, we show how the performance of its implementations can be further improved by code optimization.     

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.     

Generalized dynamic programming: A variables freezing method

A variables freezing method is designed for widening the possibilities of numerical application of generalized dynamic programming algorithms. Reducing the size of the memory under certain conditions, the method surmounts the “dimensional curse.” This paper is the continuation of [1–3].     

A volume-of-fluid method for simulation of compressible axisymmetric multi-material flow

A two-dimensional Eulerian hydrodynamic method for the numerical simulation of inviscid compressible axisymmetric multi-material flow in external force fields for the situation of pure fluids separated by macroscopic interfaces is presented. The method combines an implicit Lagrangian step with an explicit Eulerian advection step. Individual materials obey separate energy equations, fulfill general equations of state, and may possess different temperatures. Material volume is tracked using a piecewise linear volume-of-fluid method. An overshoot-free logically simple and economic material advection algorithm for cylinder coordinates is derived, in an algebraic formulation. New aspects arising in the case of more than two materials such as the material ordering strategy during transport are presented. One- and two-dimensional numerical examples are given.     

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.     

A Mathematica program for the approximate analytical solution to a nonlinear undamped Duffing equation by a new approximate approach

According to Mickens [R.E. Mickens, Comments on a Generalized Galerkin's method for non-linear oscillators, J. Sound Vib. 118 (1987) 563], the general HB (harmonic balance) method is an approximation to the convergent Fourier series representation of the periodic solution of a nonlinear oscillator and not an approximation to an expansion in terms of a small parameter. Consequently, for a nonlinear undamped Duffing equation with a driving force Bcos(ωx), to find a periodic solution when the fundamental frequency is identical to ω, the corresponding Fourier series can be written as

     

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.     

An efficient linearity and bound preserving conservative interpolation (remapping) on polyhedral meshes

An accurate conservative interpolation (remapping) algorithm is an essential component of most Arbitrary Lagrangian-Eulerian (ALE) methods. In this paper, we describe an efficient linearity and bound preserving method for polyhedral meshes. The algorithm is based on reconstruction, approximate integration and conservative redistribution. We validate our method with a suite of numerical examples, analyzing the results from the viewpoint of accuracy and order of convergence.     

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.     

A simplified gradient evaluation on non-orthogonal meshes; application to a plasma torch simulation method

A numerical method is proposed to evaluate gradients on non-uniform, non-orthogonal, three-dimensional structured meshes of hexahedra, as commonly used by finite volume methods. The method uses isoparametric transforms on tetrahedra to evaluate the gradient on a regular mesh and transform it back to the general mesh. It provides second-order accuracy, even on highly non-orthogonal meshes. Results of stationary three-dimensional numerical simulations of a direct current plasma torch, making use of the proposed method, are presented.     

How would you integrate the equations of motion in dissipative particle dynamics simulations?

In this work we assess the quality and performance of several novel dissipative particle dynamics integration schemes that have not previously been tested independently. Based on a thorough comparison we identify the respective methods of Lowe and Shardlow as particularly promising candidates for future studies of large-scale properties of soft matter systems.     

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.     

Oedometric test, Bauer's law and the micro-macro connection for a dry sand

What is the relationship between the macroscopic parameters of the constitutive equation for a granular soil and the microscopic forces between grains? In order to investigate this connection, we have simulated by molecular dynamics the oedometric compression of a granular soil (a dry and bad-graded sand) and computed the hypoplastic parameters hs (the granular skeleton hardness) and η (the exponent in the compression law) by following the same procedure than in experiments, that is by fitting the Bauer's law e/e₀=exp(−n(3p/hs)), where p is the pressure and e₀ and e are the initial and present void ratios. The micro-mechanical simulation includes elastic and dissipative normal forces plus slip, rolling and static friction between grains. By this way we have explored how the macroscopic parameters change by modifying the grains stiffness, V; the dissipation coefficient, γn; the static friction coefficient, μs; and the dynamic friction coefficient, μk. Cumulating all simulations, we obtained an unexpected result: the two macroscopic parameters seems to be related by a power law, hs=0.068(4)η9.88(3). Moreover, the experimental result for a Guamo sand with the same granulometry fits perfectly into this power law. Is this relation real? What is the final ground of the Bauer's Law? We conclude by exploring some hypothesis.     

Cluster integration method in Lagrangian particle dynamics

An efficient and robust approach is proposed in order to conduct numerical simulations of collisional particle dynamics in the Lagrangian framework. Clusters of particles are made of particles that interact or may interact during the next global time-step. Potential collision partners are found by performing a test move, that follows the patterns of a hard-sphere model. The clusters are integrated separately and the collisional forces between particles are given by a soft-sphere collision model. However, the present approach also allows longer range inter-particle forces. The integration of the clusters can be done by any one-step ordinary differential equation solver, but for dilute particle systems, the variable step-size Runge-Kutta solvers as the Dormand and Prince scheme [J. Comput. Appl. Math. 6 (1980) 19] are superior. The cluster integration method is applied on sedimentation of 5000 particles in a two-dimensional box. A significant speed-up is achieved. Compared to a traditional discrete element method with the forward Euler scheme, a speed-up factor of three orders of magnitude in the dilute regime and two orders of magnitude in the dense regime were observed. As long as the particles are dilute, the Dormand and Prince scheme is ten times faster than the classical fourth-order Runge-Kutta solver with fixed step size.     

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.     

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.     

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.     

Extended Jacobian elliptic function algorithm with symbolic computation to construct new doubly-periodic solutions of nonlinear differential equations

With the aid of computerized symbolic computation, the extended Jacobian elliptic function expansion method and its algorithm are presented by using some relations among ten Jacobian elliptic functions and are very powerful to construct more new exact doubly-periodic solutions of nonlinear differential equations in mathematical physics. The new (2+1)-dimensional complex nonlinear evolution equations is chosen to illustrate our algorithm such that sixteen families of new doubly-periodic solutions are obtained. When the modulus m→1 or 0, these doubly-periodic solutions degenerate as solitonic solutions including bright solitons, dark solitons, new solitons as well as trigonometric function solutions.     

Computing zeros of analytic functions in the complex plane without using derivatives

We present a package in Fortran 90 which solves f(z)=0, where z∈W⊂C without requiring the evaluation of derivatives, f^′(z). W is bounded by a simple closed curve and f(z) must be holomorphic within W.We have developed and tested the package to support our work in the modeling of high frequency and optical wave guiding and resonant structures. The respective eigenvalue problems are particularly challenging because they require the high precision computation of all multiple complex roots of f(z) confined to the specified finite domain. Generally f(z), despite being holomorphic, does not have explicit analytical form thereby inhibiting evaluation of its derivatives.

Program summary

Title of program:EZEROCatalogue identifier:ADXY_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXY_v1_0Program obtainable from:CPC Program Library, Queen's University of Belfast, N. IrelandComputer:IBM compatible desktop PCOperating system:Fedora Core 2 Linux (with 2.6.5 kernel)Programming languages used:Fortran 90No. of bits in a word:32No. of processors used:oneHas the code been vectorized:noNo. of lines in distributed program, including test data, etc.:21045Number of bytes in distributed program including test data, etc.:223 756Distribution format:tar.gzPeripherals used:noneMethod of solution:Our package uses the principle of the argument to count the number of zeros encompassed by a contour and then computes estimates for the zeros. Refined results for each zero are obtained by application of the derivative-free Halley method with or without Aitken acceleration, as the user wishes.     

On the numerical treatment of an ordinary differential equation arising in one-dimensional non-Fickian diffusion problems

In a recent study, Chen and Liu [Comput. Phys. Comm. 150 (2003) 31] considered a one-dimensional, linear non-Fickian diffusion problem with a potential field, which, upon application of the Laplace transform, resulted in a second-order linear ordinary differential equation which was solved by means of a control-volume finite difference method that employs exponential shape functions. It is first shown that this formulation does not properly account for the spatial dependence of the drift forces and results in oscillatory solutions near the left boundary when these forces are large. A piecewise linearized method that provides piecewise analytical solutions, is exact in exact arithmetic for constant coefficients, homogeneous, second-order linear ordinary differential equations and results in three-point finite difference equations is then proposed. Numerical simulations indicate that the piecewise linearized method is free from unphysical oscillations and more accurate than that of Chen and Liu, especially for large drift forces. The method is then applied to non-Fickian diffusion problems with non-constant drift forces in order to determine the effects of the potential field on the concentration distribution.