Discontinuities handled with events in Assimulo

Often integrating ordinary differential equations or differential algebraic equations (DAE) do not constitute the problem alone. A common complement is finding the root of an algebraic function (an event function) that depends on the states of the problem. This formulation of a model enables the possibility of including discontinuities, an important part of the Functional Mock-up Interface standard which allows hybrid models of differential algebraic equations. The problem of root-finding during integration is however difficult. Both in a theoretical aspect and as a software problem. An implementation of software for root-finding is done in Assimulo, a Python/Cython wrapper for integrators. The implementation takes the Functional Mock-up Interface standard into consideration. The implementation is made usable for a wide variety of integration algorithms and is also verified and benchmarked with advanced industrial models, showing good indications of being robust and scaling well for large systems

Connecting the Dots: Towards Continuous Time Hamiltonian Monte Carlo

Continuous time Hamiltonian Monte Carlo is introduced, as a powerful alternative to Markov chain Monte Carlo methods for continuous target distributions. The method is constructed in two steps: First Hamiltonian dynamics are chosen as the deterministic dynamics in a continuous time piecewise deterministic Markov process. Under very mild restrictions, such a process will have the desired target distribution as an invariant distribution. Secondly, the numerical implementation of such processes, based on adaptive numerical integration of second order ordinary differential equations is considered. The numerical implementation yields an approximate, yet highly robust algorithm that, unlike conventional Hamiltonian Monte Carlo, enables the exploitation of the complete Hamiltonian trajectories (hence the title). The proposed algorithm may yield large speedups and improvements in stability relative to relevant benchmarks, while incurring numerical errors that are negligible relative to the overall Monte Carlo errors

Engineering and programming manual: Two-dimensional kinetic reference computer program (TDK)

The Two Dimensional Kinetics (TDK) computer program is a primary tool in applying the JANNAF liquid rocket thrust chamber performance prediction methodology. The development of a methodology that includes all aspects of rocket engine performance from analytical calculation to test measurements, that is physically accurate and consistent, and that serves as an industry and government reference is presented. Recent interest in rocket engines that operate at high expansion ratio, such as most Orbit Transfer Vehicle (OTV) engine designs, has required an extension of the analytical methods used by the TDK computer program. Thus, the version of TDK that is described in this manual is in many respects different from the 1973 version of the program. This new material reflects the new capabilities of the TDK computer program, the most important of which are described

Perturbations of slowly rotating black holes: massive vector fields in the Kerr metric

We discuss a general method to study linear perturbations of slowly rotating black holes which is valid for any perturbation field, and particularly advantageous when the field equations are not separable. As an illustration of the method we investigate massive vector (Proca) perturbations in the Kerr metric, which do not appear to be separable in the standard Teukolsky formalism. Working in a perturbative scheme, we discuss two important effects induced by rotation: a Zeeman-like shift of nonaxisymmetric quasinormal modes and bound states with different azimuthal number m, and the coupling between axial and polar modes with different multipolar index l. We explicitly compute the perturbation equations up to second order in rotation, but in principle the method can be extended to any order. Working at first order in rotation we show that polar and axial Proca modes can be computed by solving two decoupled sets of equations, and we derive a single master equation describing axial perturbations of spin s=0 and s=+-1. By extending the calculation to second order we can study the superradiant regime of Proca perturbations in a self-consistent way. For the first time we show that Proca fields around Kerr black holes exhibit a superradiant instability, which is significantly stronger than for massive scalar fields. Because of this instability, astrophysical observations of spinning black holes provide the tightest upper limit on the mass of the photon: mv<4x10^-20 eV under our most conservative assumptions. Spin measurements for the largest black holes could reduce this bound to mv<10^-22 eV or lower.Comment: v1: 29 pages, 9 figures, 3 appendices. v2: References added and improved discussion. Matches the version to appear in Physical Review D. Mathematica notebooks available here http://blackholes.ist.utl.pt/?page=Files and http://www.phy.olemiss.edu/~berti/qnms.htm

Direct event location techniques based on Adams multistep methods for discontinuous ODEs

In this paper we consider numerical techniques to locate the event points of the differential system x′=f(x), where f is a discontinuous vector field along an event surface splitting the state space into two different regions R1 and R2 and f(x)=fi(x) when x∈Ri, for i=1,2 while f1(x)≠f2(x) when x∈Σ. Methods based on Adams multistep schemes which approach the event surface Σ from one side only and in a finite number of steps are proposed. Particularly, these techniques do not require the evaluation of the vector field f1 (respectively, f2) in the region R2 (respectively R1) and are based on the computation–at each step– of a new time ste

Birth/birth-death processes and their computable transition probabilities with biological applications

Birth-death processes track the size of a univariate population, but many biological systems involve interaction between populations, necessitating models for two or more populations simultaneously. A lack of efficient methods for evaluating finite-time transition probabilities of bivariate processes, however, has restricted statistical inference in these models. Researchers rely on computationally expensive methods such as matrix exponentiation or Monte Carlo approximation, restricting likelihood-based inference to small systems, or indirect methods such as approximate Bayesian computation. In this paper, we introduce the birth(death)/birth-death process, a tractable bivariate extension of the birth-death process. We develop an efficient and robust algorithm to calculate the transition probabilities of birth(death)/birth-death processes using a continued fraction representation of their Laplace transforms. Next, we identify several exemplary models arising in molecular epidemiology, macro-parasite evolution, and infectious disease modeling that fall within this class, and demonstrate advantages of our proposed method over existing approaches to inference in these models. Notably, the ubiquitous stochastic susceptible-infectious-removed (SIR) model falls within this class, and we emphasize that computable transition probabilities newly enable direct inference of parameters in the SIR model. We also propose a very fast method for approximating the transition probabilities under the SIR model via a novel branching process simplification, and compare it to the continued fraction representation method with application to the 17th century plague in Eyam. Although the two methods produce similar maximum a posteriori estimates, the branching process approximation fails to capture the correlation structure in the joint posterior distribution