Numerical Recipes Research Articles

Functional central limit theorems for rough volatility

The non-Markovian nature of rough volatility makes Monte Carlo methods challenging, and it is in fact a major challenge to develop fast and accurate simulation algorithms. We provide an efficient one for stochastic Volterra processes, based on an extension of Donsker’s approximation of Brownian motion to the fractional Brownian case with arbitrary Hurst exponent H∈(0,1)\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$H \\in (0,1)$\\end{document}. Some of the most relevant consequences of this ‘rough Donsker (rDonsker) theorem’ are functional weak convergence results in Skorokhod space for discrete approximations of a large class of rough stochastic volatility models. This justifies the validity of simple and easy-to-implement Monte Carlo methods, for which we provide detailed numerical recipes. We test these against the current benchmark hybrid scheme and find remarkable agreement (for a large range of values of H\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$H$\\end{document}). Our rDonsker theorem further provides a weak convergence proof for the hybrid scheme itself and allows constructing binomial trees for rough volatility models, the first available scheme (in the rough volatility context) for early exercise options such as American or Bermudan options.

Generation of Scale-Free Assortative Networks via Newman Rewiring for Simulation of Diffusion Phenomena

By collecting and expanding several numerical recipes developed in previous work, we implement an object-oriented Python code, based on the networkX library, for the realization of the configuration model and Newman rewiring. The software can be applied to any kind of network and “target” correlations, but it is tested with focus on scale-free networks and assortative correlations. In order to generate the degree sequence we use the method of “random hubs”, which gives networks with minimal fluctuations. For the assortative rewiring we use the simple Vazquez-Weigt matrix as a test in the case of random networks; since it does not appear to be effective in the case of scale-free networks, we subsequently turn to another recipe which generates matrices with decreasing off-diagonal elements. The rewiring procedure is also important at the theoretical level, in order to test which types of statistically acceptable correlations can actually be realized in concrete networks. From the point of view of applications, its main use is in the construction of correlated networks for the solution of dynamical or diffusion processes through an analysis of the evolution of single nodes, i.e., beyond the Heterogeneous Mean Field approximation. As an example, we report on an application to the Bass diffusion model, with calculations of the time tmax of the diffusion peak. The same networks can additionally be exported in environments for agent-based simulations like NetLogo.

Quantum chemical “Aufbau” principles: how to estimate the shape of highly flexible (bio-)polymers? A recursively extendable “chemion picture” of Euler-Hückel-type

An outline is given of how to split the n-dimensional space of torsion angles occurring in flexible (bio-)polymers (like alkanes, nucleic acids, or proteins, for instance) into n one-dimensional potential curves. Forthcoming applications will focus on the “protein folding problem,” beginning with polyglycine.ContextIn accordance with Euler’s topology rules, molecules are considered to be composed of “vertices” (atoms, ligands, bonding sites, functional groups, and bigger fragments). Following Hückel, each vertex is represented by only one basis function. Starting from the “monofocal” hydrids CH4\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\ ext {CH}_{4}$$\\end{document}, NH3\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\ ext {NH}_{3}$$\\end{document}, OH2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\ ext {OH}_{2}$$\\end{document}, FH, and SiH4\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\ ext {SiH}_{4}$$\\end{document}, PH3\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\ ext {PH}_{3}$$\\end{document}, SH2\\documentclass[12pt]{minimal} \\usepackage{amsmath} \\usepackage{wasysym} \\usepackage{amsfonts} \\usepackage{amssymb} \\usepackage{amsbsy} \\usepackage{mathrsfs} \\usepackage{upgreek} \\setlength{\\oddsidemargin}{-69pt} \\begin{document}$$\ ext {SH}_{2}$$\\end{document}, ClH as anchor units, “chemionic” Hamiltonians (of individual “chemion ensembles” and proportional nuclear charges) are constructed recursively, together with an appropriate basis set for the first five (normal) alkanes and some related oligomers like primary alcohols, alkyl amines, and alkyl chlorides.MethodsStandard methods (“Restricted Hartree-Fock RHF” and “Full Configuration Interaction FCI”) are used to solve the various stationary Schrödinger equations. Two software packages are indispensable: “SMILES” for integral evaluations over Slater-type orbitals (STO), and “Numerical Recipes” for matrix diagonalizations and inversions. While managing with only two-center repulsion integrals, “implicit multi-center integrations” lead us to the non-empirical fundament of Hoffmann’s “Extended-Hückel Theory.”

Numerical recipes of virtual element method for phase field modeling of brittle fracture

In this work, a new and efficient virtual element formulation for non-standard phase field model of brittle fracture is presented. A multi-pass alternative minimization solution scheme based on algorithm operator splitting is utilized, which decouples the whole problem into two parts, namely, mechanical and damage sub-problems. The former is treated as elasto-static problem, while the latter one is treated as Poisson-type of reaction-diffusion equation subjected to bounded and irreversibility constraint. To demonstrate the performance of proposed formulation, several benchmark problems are studied and results are in good agreement with corresponding finite element calculations and experimental studies.

TEMPus VoLA: The timed Epstein multi-pressure vessel at low accelerations.

The field of planetary system formation relies extensively on our understanding of the aerodynamic interaction between gas and dust in protoplanetary disks. Of particular importance are the mechanisms triggering fluid instabilities and clumping of dust particles into aggregates, and their subsequent inclusion into planetesimals. We introduce the timed Epstein multi-pressure vessel at low accelerations, which is an experimental apparatus for the study of particle dynamics and rarefied gas under micro-gravity conditions. This facility contains three experiments dedicated to studying aerodynamic processes: (i) the development of pressure gradients due to collective particle-gas interaction, (ii) the drag coefficients of dust aggregates with variable particle-gas velocity, and (iii) the effect of dust on the profile of a shear flow and resultant onset of turbulence. The approach is innovative with respect to previous experiments because we access an untouched parameter space in terms of dust particle packing fraction, and Knudsen, Stokes, and Reynolds numbers. The mechanisms investigated are also relevant for our understanding of the emission of dust from active surfaces, such as cometary nuclei, and new experimental data will help interpreting previous datasets (Rosetta) and prepare future spacecraft observations (Comet Interceptor). We report on the performance of the experiments, which has been tested over the course of multiple flight campaigns. The project is now ready to benefit from additional flight campaigns, to cover a wide parameter space. The outcome will be a comprehensive framework to test models and numerical recipes for studying collective dust particle aerodynamics under space-like conditions.

A Collection of Numerical Recipes Useful for Building Scalable Psychometric Applications

This article is concerned with a subset of numerically stable and scalable algorithms useful to support computationally complex psychometric models in the era of machine learning and massive data. The subset selected here is a core set of numerical methods that should be familiar to computational psychometricians and considers whitening transforms for dealing with correlated data, computational concepts for linear models, multivariable integration, and optimization techniques.

CARAMEL-gas: A Step toward Fast Empirical Models of the Broad-line-emitting Gas

Fast empirical models of the broad emission line region (BLR) are a powerful tool to interpret velocity-resolved reverberation mapping (RM) data, estimate the mass of the supermassive black holes, and gain insight into its geometry and kinematics. Much of the effort so far has been devoted to describing the emissivity of one emission line at a time. We present here an alternative approach aimed at describing the underlying BLR gas distribution, by exploiting simple numerical recipes to connect it with emissivity. This approach is a step toward describing multiple emission lines originating from the same gas and allows us to clarify some issues related to the interpretation of RM data. We illustrate this approach—implemented in the code CARAMEL-gas—using three data sets covering the Hβ emission line (Mrk 50, Mrk 1511, Arp 151) that have been modeled using the emissivity-based version of the code. As expected, we find differences in the parameters describing the BLR gas and emissivity distribution, but the emissivity-weighted lag measurements and all other model parameters including black hole mass and overall BLR morphology and kinematics are consistent with the previous measurements. We also model the Hα emission line for Arp 151 using both the gas- and emissivity-based BLR models. We find ionization stratification in the BLR with Hα arising at larger radii than Hβ, while all other model parameters are consistent within the uncertainties.

Synchronization of dynamical systems on Riemannian manifolds by an extended PID-type control theory: Numerical evaluation

<p style='text-indent:20px;'>The present document outlines a non-linear control theory, based on the PID regulation scheme, to synchronize two second-order dynamical systems insisting on a Riemannian manifold. The devised extended PID scheme, referred to as M-PID, includes an unconventional component, termed 'canceling component', whose purpose is to cancel the natural dynamics of a system and to replace it with a desired dynamics. In addition, this document presents numerical recipes to implement such systems, as well as the devised control scheme, on a computing platform and a large number of numerical simulation results focused on the synchronization of Duffing-like non-linear oscillators on the unit sphere. Detailed numerical evaluations show that the canceling contribution of the M-PID control scheme is not critical to the synchronization of two oscillators, however, it possesses the beneficial effect of speeding up their synchronization. Simulation results obtained in non-ideal conditions, namely in the presence of additive disturbances and delays, reveal that the devised synchronization scheme is robust against high-frequency additive disturbances as well as against observation delays.</p>

Numerical recipes for faster MAS-DNP simulations

Numerical simulations of Magic Angle Spinning Dynamic Nuclear Polarization (MAS-DNP) have transformed the way the DNP process is understood in rotating samples. In 2012, two methods were concomitantly developed to simulate small spin systems (< 4 spin-1/2). The development of new polarizing agents, including those containing metal centers with S > 1/2, makes it necessary to further expand the numerical tools with minimal approximations that will help rationalize the experimental observations and build approximate models. In this paper, three strategies developed in the past five years are presented: an adaptive integration scheme, a hybrid Hilbert/Liouville formalism, and a method to truncate the Liouville space basis for periodic Hamiltonian. Each of these methods enable time savings ranging from a factor of 3 to > 100. We illustrate the code performance by reporting for the first time the MAS-DNP field profiles for “AMUPol”, in which the couplings to the nitrogen nuclei are explicitly considered, as well as Cross-Effect MAS-DNP field profiles with two electrons spin 5/2 interacting with a nuclear spin 1/2.

A practical guide to writing a radiative transfer code

Using our decades-long experience in radiative transfer (RT) code development for Earth science, we endeavor to reduce the knowledge gap of bringing RT from theory to code quickly. Despite numerous classic and recent literature, it is still hard to develop an RT code from scratch within a few weeks. It is equally hard to understand, not to mention modify, an existing “monster” RT code, for which the developer is either located remotely or has retired. Following the format of “Numerical Recipes” by Press et al., we collocate in this paper small pieces of necessary theory with corresponding small pieces of RT code. These are arranged in an order that is natural for code development, which is often opposite of the natural order for laying out the theoretical basis. We focus on the transfer of unpolarized monochromatic solar radiation in a plane-parallel atmosphere over a reflecting surface. Both the surface and the atmosphere are homogeneous (uniform) at all directions. The multiple scattering is numerically solved using the deterministic method of Gauss-Seidel iterations. Except for the presented Python-Numba open-source RT code gsit, the paper does not report any new scientific results, but rather serves as an academic demonstration. If development time is an issue or the reader is familiar with basic concepts of RT theory, we recommend proceeding directly to Sec. 3 “RT code development”. Program summaryProgram title: gsit (pronounced “jeezit”)CPC Library link to program files:https://doi.org/10.17632/d3zt5zhx49.1Developer's repository link:https://github.com/korkins/gsitLicensing provisions: MITProgramming language: Python 3Nature of problem: We present a tutorial code in Python for deterministic (non-stochastic) numerical simulation of multiple scattering of monochromatic solar light in a plane-parallel Earth atmosphere bounded from below by a reflecting surface. The problem is solved in a simplified form (i.e., uniform atmosphere, no polarization, uniform surface reflectance, etc.) to better explain numerical features, rather than physics, of propagation of light in the atmosphere.Solution method: The method of Gauss-Seidel iterations. It relies on the Fourier decomposition of the Radiative Transfer Equation over azimuth, Gauss quadrature for numerical integration over the zenith and iterative process for integration over height (optical depth) with analytical (hence known) single scattering approximation being the starting point. The method is relatively simple to code and does not require any external libraries.

Solving Schrödinger Bridges via Maximum Likelihood.

The Schrödinger bridge problem (SBP) finds the most likely stochastic evolution between two probability distributions given a prior stochastic evolution. As well as applications in the natural sciences, problems of this kind have important applications in machine learning such as dataset alignment and hypothesis testing. Whilst the theory behind this problem is relatively mature, scalable numerical recipes to estimate the Schrödinger bridge remain an active area of research. Our main contribution is the proof of equivalence between solving the SBP and an autoregressive maximum likelihood estimation objective. This formulation circumvents many of the challenges of density estimation and enables direct application of successful machine learning techniques. We propose a numerical procedure to estimate SBPs using Gaussian process and demonstrate the practical usage of our approach in numerical simulations and experiments.

Generic Models for Disk-resolved and Disk-integrated Phase-dependent Linear Polarization of Light Reflected from Exoplanets

Abstract Similar to the case of solar system planets, reflected starlight from exoplanets is expected to be polarized due to atmospheric scattering, and the net disk-integrated polarization should be nonzero owing to the asymmetrical illumination of the planetary disk. The computation of the disk-integrated reflected flux and its state of polarization involves techniques for the calculation of the local reflection matrices, as well as the numerical recipes for integration over the planetary disks. In this paper, we present a novel approach to calculate the azimuth-dependent reflected intensity vectors at each location on the planetary disk divided into grids. We achieve this by solving the vector radiative transfer equations that describe linear polarization. Our calculations incorporate self-consistent atmospheric models of exoplanets over a wide range of equilibrium temperature, surface gravity, atmospheric composition, and cloud structure. A comparison of the flux and the amount of polarization calculated by considering both single and multiple scattering exhibits the effect of depolarization due to multiple scattering of light depending on the scattering albedo of the atmosphere. We have benchmarked our basic calculations against some of the existing models. We have also presented our models for the hot Jupiter HD 189733 b, indicating the level of precision required by future observations to detect the polarization of this planet in the optical and near-infrared wavelength region. The generic nature and the accuracy offered by our models make them an effective tool for modeling future observations of the polarized light reflected from exoplanets.

Electrostatic potential energy stored in a hemispherical surface with uniform surface charge distribution

We report an exact result for the electrostatic potential energy stored in a hemispherical surface with uniform surface charge density. This system lacks the spherical symmetry of a uniformly charged spherical surface or its solid sphere counterpart. Consequently, a hemispherical surface with all its charge uniformly spread on its surface represents an interesting case study in electrostatics. In this work we combine the fact that the body has axial symmetry with clever mathematical tools to obtain an outcome that is remarkably simple. The final exact expression obtained for the electrostatic potential energy stored in a hemispherical surface with uniform surface charge density answers a long-standing question. Though the bulk of the treatment is analytical, the accuracy of the findings can also be ascertained by using symbolic computation software or numerical recipes. Hence, the present analytical results can be compared to ones obtained with various powerful numerical methods which nowadays are routinely used to check the correctness of analytical and semi-analytical calculations.

Real Gas Models in Coupled Algorithms Numerical Recipes and Thermophysical Relations

In the majority of compressible flow CFD simulations, the standard ideal gas state equation is accurate enough. However, there is a range of applications where the deviations from the ideal gas behaviour is significant enough that performance predictions are no longer valid and more accurate models are needed. While a considerable amount of the literature has been written about the application of real gas state equations in CFD simulations, there is much less information on the numerical issues involved in the actual implementation of such models. The aim of this article is to present a robust implementation of real gas flow physics in an in-house, coupled, pressure-based solver, and highlight the main difference that arises as compared to standard ideal gas model. The consistency of the developed iterative procedures is demonstrated by first comparing against results obtained with a framework using perfect gas simplifications. The generality of the developed framework is tested by using the parameters from two different real gas state equations, namely the IAPWS-97 and the cubic state equations state equations. The highly polynomial IAPWS-97 formulation for water is applied to a transonic nozzle case where steam is expanded at transonic conditions until phase transition occurs. The cubic state equations are applied to a two stage radial compressor setup. Results are compared in terms of accuracy with a commercial code and measurement data. Results are also compared against simulations using the ideal gas model, highlighting the limitations of the later model. Finally, the effects of the real gas formulations on computational time are compared with results obtained using the ideal gas model.

Numerical recipes for elastodynamic virtual element methods with explicit time integration

SummaryWe present a general framework to solve elastodynamic problems by means of the virtual element method (VEM) with explicit time integration. In particular, the VEM is extended to analyze nearly incompressible solids using the B‐bar method. We show that, to establish a B‐bar formulation in the VEM setting, one simply needs to modify the stability term to stabilize only the deviatoric part of the stiffness matrix, which requires no additional computational effort. Convergence of the numerical solution is addressed in relation to stability, mass lumping scheme, element size, and distortion of arbitrary elements, either convex or nonconvex. For the estimation of the critical time step, two approaches are presented, ie, the maximum eigenvalue of a system of mass and stiffness matrices and an effective element length. Computational results demonstrate that small edges on convex polygonal elements do not significantly affect the critical time step, whereas convergence of the VEM solution is observed regardless of the stability term and the element shape in both two and three dimensions. This extensive investigation provides numerical recipes for elastodynamic VEMs with explicit time integration and related problems.

Finite element based Green's function integral equation for modelling light scattering.

We revisit the Green's function integral equation for modelling light scattering with discretization strategies as well as numerical integration recipes borrowed from finite element method. The finite element based Green's function integral equation is implemented by introducing auxiliary variables, which are used to discretize the Green's function integral equation. The merits of introducing finite element techniques into Green's function integral equation are apparent. Firstly, the finite element discretization provides a better geometric approximation of the scatterers, compared with that of the conventional discretization method using staircase approximation. Secondly, the accuracy of numerical integral inside one element associated with Green's function integral equations can be improved by using more quadrature points, where the singular terms confined inside each triangle can be approximated analytically. We then illustrate the advantages of our finite element based Green's function integral equation method via a few concrete examples in modelling light scattering by optically large and complex scatterers in 2-dimensional scenarios.

Near-Wall Flow in Cerebral Aneurysms

The region where the vascular lumen meets the surrounding endothelium cell layer, hence the interface region between haemodynamics and cell tissue, is of primary importance in the physiological functions of the cardiovascular system. The functions include mass transport to/from the blood and tissue, and signalling via mechanotransduction, which are primary functions of the cardiovascular system and abnormalities in these functions are known to affect disease formation and vascular remodelling. This region is denoted by the near-wall region in the present work, and we outline simple yet effective numerical recipes to analyse the near-wall flow field. Computational haemodynamics solutions are presented for six patient specific cerebral aneurysms, at three instances in the cardiac cycle: peak systole, end systole (taken as dicrotic notch) and end diastole. A sensitivity study, based on Newtonian and non-Newtonian rheological models, and different flow rate profiles, is effected for a selection of aneurysm cases. The near-wall flow field is described by the wall shear stress (WSS) and the divergence of wall shear stress (WSSdiv), as descriptors of tangential and normal velocity components, respectively, as well as the wall shear stress critical points. Relations between near-wall and free-stream flow fields are discussed.

Numerical methods as an introduction to calculus

This paper develops the ideas of The Pedagogical Power of Context: Iterative Calculus Methods and the Epidemiology of Eyam (French et al 2018 J. Phys. Educ.), where we considered the application of the Euler method to solve epidemiological problems. Our purpose was to convey some examples of school level work in numerical methods that can provide useful visualisations to aid in comprehension and to incentivise the acquisition of skills in calculus. We recommend that the modern pre-university student of the physical sciences would do well to develop the precursors to coding, and learn the basics of modelling using spreadsheets in tandem with analytical methods that dominate contemporary schemes of work. In the present paper, we consider how the first-order Euler (Press et al 2003 Numerical Recipes in C++: the Art of Scientific Computing 2nd edn (Cambridge: Cambridge University Press) pp 712–55) iterative method may be extended to provide a useful platform to embark on introductory work in more accurate Verlet (https://en.wikipedia.org/wiki/Verlet_integration (Accessed: 2 September 2018)) techniques, particularly with applications to the solution of second order ordinary differential equations.

Numerical recipes for successfully modeling the phase transitions in thermal energy storage adsorption systems

AbstractStrong convergence problems may arise during the dynamic simulation of the phase transitions of adsorption systems subjected to larger and/or random fluctuations of the boundary conditions, which do not follow the conventional cycle sequence (e.g., adsorptive thermal storage systems). This sequence could even be reversed at some instants and its stages not completed. Such convergence problems preclude solving those phase transitions, and thus the successful dynamic simulation of the adsorption systems. A set of numerical procedures for solving those problems is proposed, consisting of a small group of new conditions to be added to the conventional numerical algorithms. Its successful use is demonstrated in the dynamic simulation of a thermal energy storage system provided with an adsorption module, modeled in two distinct but coupled platforms (TRNSYS and MATLAB), for which the convergence difficulties are identified and described. Results indicate successful convergence by applying the proposed methodology, in contrast with the original algorithm, attesting for its capability to overcome the convergence problems that occur during the system's dynamic simulation. This set of numerical recipes can be adapted for other systems' algorithms, or to serve as a basis for its successful dynamic simulation, thus showing a promising and enlarged range of use.

Asymptotics for Volatility Derivatives in Multi-Factor Rough Volatility Models

We present small-time implied volatility asymptotics for Realised Variance (RV) and VIX options for a number of (rough) stochastic volatility models via large deviations principle. We provide numerical results along with efficient and robust numerical recipes to compute the rate function; the backbone of our theoretical framework. Based on our results, we further develop approximation schemes for the density of RV, which in turn allows to express the volatility swap in close-form. Lastly, we investigate different constructions of multi-factor models and how each of them affects the convexity of the implied volatility smile. Interestingly, we identify the class of models that generate non-linear smiles around-the-money.