Accelerating numerical solution of stochastic differential equations with CUDA

Numerical integration of stochastic differential equations is commonly used in many branches of science. In this paper we present how to accelerate this kind of numerical calculations with popular NVIDIA Graphics Processing Units using the CUDA programming environment. We address general aspects of numerical programming on stream processors and illustrate them by two examples: the noisy phase dynamics in a Josephson junction and the noisy Kuramoto model. In presented cases the measured speedup can be as high as 675× compared to a typical CPU, which corresponds to several billion integration steps per second. This means that calculations which took weeks can now be completed in less than one hour. This brings stochastic simulation to a completely new level, opening for research a whole new range of problems which can now be solved interactively.

Program summary

Program title: SDECatalogue identifier: AEFG_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEFG_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: Gnu GPL v3No. of lines in distributed program, including test data, etc.: 978No. of bytes in distributed program, including test data, etc.: 5905Distribution format: tar.gzProgramming language: CUDA CComputer: any system with a CUDA-compatible GPUOperating system: LinuxRAM: 64 MB of GPU memoryClassification: 4.3External routines: The program requires the NVIDIA CUDA Toolkit Version 2.0 or newer and the GNU Scientific Library v1.0 or newer. Optionally gnuplot is recommended for quick visualization of the results.Nature of problem: Direct numerical integration of stochastic differential equations is a computationally intensive problem, due to the necessity of calculating multiple independent realizations of the system. We exploit the inherent parallelism of this problem and perform the calculations on GPUs using the CUDA programming environment. The GPU's ability to execute hundreds of threads simultaneously makes it possible to speed up the computation by over two orders of magnitude, compared to a typical modern CPU.Solution method: The stochastic Runge-Kutta method of the second order is applied to integrate the equation of motion. Ensemble-averaged quantities of interest are obtained through averaging over multiple independent realizations of the system.Unusual features: The numerical solution of the stochastic differential equations in question is performed on a GPU using the CUDA environment.Running time: < 1 minute     

A stable algorithm for Hankel transforms using hybrid of Block-pulse and Legendre polynomials

A new numerical method, based on hybrid of Block-pulse and Legendre polynomials for numerical evaluation of Hankel transform is proposed in this paper. Hybrid of Block-pulse and Legendre polynomials are used as a basis to expand a part of the integrand, rf(r), appearing in the Hankel transform integral. Thus transforming the integral into a Fourier-Bessel series. Truncating the series, an efficient algorithm is obtained for the numerical evaluations of the Hankel transforms of order ν>−1. The method is quite accurate and stable, as illustrated by given numerical examples with varying degree of random noise terms εθi added to the data function f(r), where θi is a uniform random variable with values in [−1,1]. Finally, an application of the proposed method is given in solving the heat equation in an infinite cylinder with a radiation condition.     

Fast spherical Bessel transform via fast Fourier transform and recurrence formula

We propose a new method for the numerical evaluation of the spherical Bessel transform. A formula is derived for the transform by using an integral representation of the spherical Bessel function and by changing the integration variable. The resultant algorithm consists of a set of the Fourier transforms and numerical integrations over a linearly spaced grid of variable k in Fourier space. Because the k-dependence appears in the upper limit of the integration range, the integrations can be performed effectively in a recurrence formula. Several types of atomic orbital functions are transformed with the proposed method to illustrate its accuracy and efficiency, demonstrating its applicability for transforms of general order with high accuracy.     

Application of the Feynman-Kac path integral method in finding excited states of quantum systems

Group theory considerations and properties of a continuous path are used to define a failure tree procedure for finding eigenvalues of the Schrödinger equation using stochastic methods. The procedure is used to calculate the lowest excited state eigenvalues of eigenfunctions possessing anti-symmetric nodal regions in configuration space using the Feynman-Kac path integral method. Within this method the solution of the imaginary time Schrödinger equation is approximated by random walk simulations on a discrete grid constrained only by symmetry considerations of the Hamiltonian. The required symmetry constraints on random walk simulations are associated with a given irreducible representation and are found by identifying the eigenvalues for the irreducible representation corresponding to symmetric or antisymmetric eigenfunctions for each group operator. The method provides exact eigenvalues of excited states in the limit of infinitesimal step size and infinite time. The numerical method is applied to compute the eigenvalues of the lowest excited states of the hydrogenic atom that transform as Γ² and Γ⁴ irreducible representations. Numerical results are compared with exact analytical results.     

Feynman Integral Evaluation by a Sector decomposiTion Approach (FIESTA)

We present a new program performing the sector decomposition and integrating the expression afterwards. The program takes a set of propagators and a set of indices as input and returns the epsilon-expansion of the corresponding integral.

Program summary

Program title: FIESTACatalogue identifier: AECP_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AECP_v1_0.htmlProgram obtainable from: CPC Program Library, Queen's University, Belfast, N. IrelandLicensing provisions: GPL v2No. of lines in distributed program, including test data, etc.: 88 281No. of bytes in distributed program, including test data, etc.: 6 153 480Distribution format: tar.gzProgramming language: Wolfram Mathematica 6.0 [3] and CComputer: from a desktop PC to supercomputerOperating system: Unix, Linux, WindowsRAM: depends on the complexity of the problemClassification: 4.4, 4.12, 5, 6.5External routines: QLink [1], Vegas [2]Nature of problem: The sector decomposition approach to evaluating Feynman integrals falls apart into the sector decomposition itself, where one has to minimize the number of sectors; the pole resolution and epsilon expansion; and the numerical integration of the resulting expression.Solution method: The sector decomposition is based on a new strategy. The sector decomposition, pole resolution and epsilon-expansion are performed in Wolfram Mathematica 6.0 [3]. The data is stored on hard disk via a special program, QLink [1]. The expression for integration is passed to the C-part of the code, that parses the string and performs the integration by the Vegas algorithm [2]. This part of the evaluation is perfectly parallelized on multi-kernel computers.Restrictions: The complexity of the problem is mostly restricted by the CPU time required to perform the evaluation of the integral, however there is currently a limit of maximum 11 positive indices in the integral; this restriction is to be removed in future versions of the code.Running time: Depends on the complexity of the problem.References:[1] http://qlink08.sourceforge.net, open source.[2] G.P. Lepage, The Cornell preprint CLNS-80/447, 1980.[3] http://www.wolfram.com/products/mathematica/index.html2.     

Numerical simulation of two-dimensional sine-Gordon solitons via a split cosine scheme

This paper proposes a split cosine scheme for simulating solitary solutions of the sine-Gordon equation in two dimensions, as it arises, for instance, in rectangular large-area Josephson junctions. The dispersive nonlinear partial differential equation allows for soliton-type solutions, a ubiquitous phenomenon in a large variety of physical problems. The semidiscretization approach first leads to a system of second-order nonlinear ordinary differential equations. The system is then approximated by a nonlinear recurrence relation which involves a cosine function. The numerical solution of the system is obtained via a further application of a sequential splitting in a linearly implicit manner that avoids solving the nonlinear scheme at each time step and allows an efficient implementation of the simulation in a locally one-dimensional fashion. The new method has potential applications in further multi-dimensional nonlinear wave simulations. Rigorous analysis is given for the numerical stability. Numerical demonstrations for colliding circular solitons are given.     

On the hilbert transform of bounded bandlimited signals

In this paper we analyze the Hilbert transform and existence of the analytical signal for the space B _?? ^?? of bandlimited signals that are bounded on the real axis. Originally, the theory was developed for signals in L ²(?) and then extended to larger signal spaces. While it is well known that the common integral representation of the Hilbert transform may diverge for some signals in B _?? ^?? and that the Hilbert transform is not a bounded operator on B _?? ^?? , it is nevertheless possible to define the Hilbert transform for the space B _?? ^?? . We use a definition that is based on the H ¹-BMO(?) duality. This abstract definition, which can be used for general bounded signals, gives no constructive procedure to compute the Hilbert transform. However, for the practically important special case of bounded bandlimited signals, we can provide such an explicit procedure by giving a closed-form expression for the Hilbert transform. Further, it is shown that the Hilbert transform of a signal in B _?? ^?? is still bandlimited but not necessarily bounded. With these results we continue the work of [1, 2].     

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.     

希尔伯特-黄变换在表音故障诊断的应用

针对钟表表音信号提取与分析存在较大困难的问题,提出了双麦克降噪采集装置和基于经验模态分解的手表表音希尔伯特-黄变换(Hilbert-Huang)方法。该方法提取机械手表的故障信号进行经验模态分解,进而对内禀函数进行希尔伯特变换得到希尔伯特谱和希尔伯特边际谱。仿真实验结果表明,边际谱能识别出故障信息,该方法能够定位误差并实现故障诊断。     

Characterization of the peak value behavior of the Hilbert transform of bounded bandlimited signals

The peak value of a signal is a characteristic that has to be controlled in many applications. In this paper we analyze the peak value of the Hilbert transform for the space $\mathcal{B}_\pi ^\infty $ of bounded bandlimited signals. It is known that for this space the Hilbert transform cannot be calculated by the common principal value integral, because there are signals for which it diverges everywhere. Although the classical definition fails for $\mathcal{B}_\pi ^\infty $ , there is a more general definition of the Hilbert transform, which is based on the abstract H ¹-BMO(?) duality. It was recently shown in [1] that, in addition to this abstract definition, there exists an explicit formula for calculating the Hilbert transform. Based on this formula we study properties of the Hilbert transform for the space $\mathcal{B}_\pi ^\infty $ of bounded bandlimited signals. We analyze its asymptotic growth behavior, and thereby solve the peak value problem of the Hilbert transform for this space. Further, we obtain results for the growth behavior of the Hilbert transform for the subspace $\mathcal{B}_{\pi ,0}^\infty $ of bounded bandlimited signals that vanish at infinity. By studying the properties of the Hilbert transform, we continue the work [2].     

基于Hilbert-Huang Transform的心音信号谱分析

心音信号是一种典型的非平稳信号,传统信号处理方法的应用受到很大限制.针对此本文提出了基于Hilbert-Huang Transform(HHT) 的心音信号的分析方法,对冠心病患者的心音信号进行了分析.通过把心音信号分解为内蕴模式函数,利用Hilbert变换建立了心音信号的时间-频率-能量三维Hilbert谱分布以及边界谱分布;Hilbert谱及其边界谱在时域以及频域以较高的分辨率表征了心音信号的时频变化特性,揭示了冠心病患者心音信号的病理特征;为冠心病的早期无损诊断奠定了坚实基础,临床实践中有较大的指导价值.     

Numerical solution of the Hilbert transform for phase calculation from an amplitude spectrum

A method for the numerical solution of the Hilbert transform integral to obtain the phase corresponding to a given amplitude spectrum is presented. The method, which is based on the assumption that the amplitude spectrum at high and low frequencies can be approximated by constant slopes, can be used to calculate the phase over the entire frequency range for both lowpass, bandpass, and highpass characteristics. The numerical solution is carried out on a minicomputer by a Fortran IV program, and the calculation error can be brought down to the level of the truncation error. The usefulness of the method combined with FFT has been shown by calculating the step response of an amplifier from its measured amplitude spectrum.     

Minimax state estimation for linear discrete-time differential-algebraic equations

This paper presents a state estimation approach for an uncertain linear equation with a non-invertible operator in Hilbert space. The approach addresses linear equations with uncertain deterministic input and noise in the measurements, which belong to a given convex closed bounded set. A new notion of a minimax observable subspace is introduced. By means of the presented approach, new equations describing the dynamics of a minimax recursive estimator for discrete-time non-causal differential-algebraic equations (DAEs) are presented. For the case of regular DAEs it is proved that the estimator’s equation coincides with the equation describing the seminal Kalman filter. The properties of the estimator are illustrated by a numerical example.     

Time resolution of Josephson balanced comparators with shunt resistance

In this paper, the time resolution of Josephson balanced comparators with shunt resistance is theoretically estimated. A general expression is obtained for the time resolution of the linearly growing strobe pulse. It is shown that the time resolution deteriorates by a factor of (1+2R _N /R _S ^8/15 (R _S is the shunt resistance, and R _N is the normal resistance of the Josephson junction), while shunting improves the sensitivity of balanced comparators, as is proven in [20]. Estimates show that it is possible to reach subpicosec time resolution for balanced Josephson comparators with shunt resistance.     

An improved deconvolution method for bremsstrahlung spectra from hot plasmas

The influence of tip region cross sections and the energy resolution of X-ray detectors on the deconvolution of X-ray spectra from hot plasmas into the electronic distribution function is discussed. From the discussion an improved deconvolution method is derived starting from an integral equation of the second kind. Here numerical discontinuities can be handled by a differentiation process, the reliability of which increases with increasing refinement of the measurements.     

Vscape V1.1.0. An interactive tool for metastable vacua

Vscape is an interactive tool for studying the one-loop effective potential of an ungauged supersymmetric model of chiral multiplets. The program allows the user to define a supersymmetric model by specifying the superpotential. The F-terms and the scalar and fermionic mass matrices are calculated symbolically. The program then allows you to search numerically for (meta)stable minima of the one-loop effective potential. Additional commands enable you to further study specific minima, by, e.g., computing the mass spectrum for those vacua. Vscape combines the flexibility of symbolic software, with the speed of a numerical package.

Program summary

Program title:Vscape 1.1.1Catalogue identifier: ADZW_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADZW_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.: 80 507No. of bytes in distributed program, including test data, etc.: 6 708 938Distribution format: tar.gzProgramming language: C++Computer: Pentium 4 PC Computers: need (GNU) C++ compiler, Linux standard GNU installation (./configure; make; make install). A precompiled Windows XP version is included in the distribution packageOperating system: Linux, Windows XP using cygwinRAM: 10 MBWord size: 32 bitsClassification: 11.6External routines: GSL (http://www.gnu.org/software/gsl/), CLN (http://www.ginac.de/CLN/), GiNaC (http://directory.fsf.org/GiNaC.html)Nature of problem:Vscape is an interactive tool for studying the one-loop effective potential of an ungauged supersymmetric model of chiral multiplets. The program allows the user to define a supersymmetric model by specifying the superpotential. The F-terms and the scalar and fermionic mass matrices are calculated symbolically. The program then allows you to search numerically for (meta)stable minima of the one-loop effective potential. Additional commands enable you to further study specific minima, by, e.g., computing the mass spectrum for those vacua. Vscape combines the flexibility of symbolic software with the speed of a numerical package.Solution method: Coleman-Weinberg potential is computed using numerical matrix diagonalization. Minima of the one-loop effective potential are found using the Nelder and Mead simplex algorithm. The one-loop effective potential can be studied using numerical differentiation. Symbolic users interface implemented using flex and bison.Restrictions:N=1 supersymmetric chiral models onlyUnusual features: GiNaC (+CLN), GSL, ReadLib (not essential)Running time: Interactive users interface. Most commands execute in a few ms. Computationally intensive commands execute in order of minutes, depending on the complexity of the user defined model.     

Calcolo numerico delle trasformate di hilbert ed applicazioni

The well known sampling theorm (7 § 1) is obtained simply by application of a Parseval equality, under the assumption that the Fourier transformF(ω) off(t) is zero ontside a finite interval and is bounded and integrable. By operating on the samples off(t) it is also possible to compute its Hilbert transformg(t). IfF(ω) is not zero outside a finite interval, the series on the right side of 13 § 1 approximates tof(t) and the series 6 § 2 approximates tog(t). Bounds for the errors are given. A formula particularly useful for the numerical calculus of the Hilbert transforms is found (8 § 2) and applications are suggested for the computation of special functions (§ 3). Finally the Wiener-Hopf integral equation is numerically solved (§ 4) by making use of Hilbert transforms.     

EEMD分解在电力系统故障信号检测中的应用

针对经验模态分解(EMD)的希尔伯特-黄变换(HHT)在电力系统故障信号检测问题,应用存在的模态混叠会导致扰动信号检测失效,为此提出一种基于聚类经验模型分解(EEMD)的故障信号检测的方法。方法通过多次对目标数据加入随机白噪声序列以保证不同区域信号映射的完整性,并且克服了传统EMD分解造成的模态混叠问题,通过EEMD方法提取信号的固有模态函数(IMF),再进行Hilbert变换,利用Hilbert谱对故障暂态和扰动时刻进行检测,通过瞬时频率实现对故障暂态和扰动时刻的准确定位。通过数字仿真分析表明,方法是准确有效的。     

Convergence and accuracy of the path integral approach for elastostatics

This paper addresses convergence rate and accuracy of a numerical technique for linear elastostatics based on a path integral formulation [Int. J. Numer. Math. Eng. 47 (2000) 1463]. The computational implementation combines a simple polynomial approximation of the displacement field with an approximate statement of the exact evolution equations, which is designated as functional integral method. A convergence analysis is performed for some simple nodal arrays. This is followed by two different estimations of the optimum parameter ζ: one is based on statistical arguments and the other on inspection of third order residuals. When the eight closest neighbors to a node are used for polynomial approximation the optimum parameter is found to depend on Poisson's ratio and lie in the range 0.5<ζ<1.5. Two well established numerical methods are then recovered as specific instances of the FIM. The strong formulation––point collocation––corresponds to the limit ζ=0 while bilinear finite elements corresponds exactly to the choice ζ=0.5. The use of the optimum parameter provides better precision than the other two methods with similar computational cost. Other nodal arrays are also studied both in two and three dimensions and the performance of the FIM compared with the corresponding finite element and collocation schemes. Finally, the implementation of FIM on unstructured meshes is discussed, and a numerical example solving Laplace equation is analyzed. It is shown that FIM compares favorably with FEM and offers a number of advantages.