A variable stepsize implementation for stochastic differential equations

Stochastic differential equations (SDEs) arise from physical systems where the parameters describing the system can only be estimated or are subject to noise. Much work has been done recently on developing higher order Runge-Kutta methods for solving SDEs numerically. Fixed stepsize implementations of numerical methods have limitations when, for example, the SDE being solved is stiff as this forces the stepsize to be very small. This paper presents a completely general variable stepsize implementation of an embedded Runge Kutta pair for solving SDEs numerically; in this implementation, there is no restriction on the value used for the stepsize, and it is demonstrated that the integration remains on the correct Brownian path.

Adams-Type Methods for the Numerical Solution of Stochastic Ordinary Differential Equations

Stochastic differential equations (SDEs) arise fi om physical systems where the parameters describing the system can only be estimated or are subject to noise. There has been much work done recently on developing numerical methods for solving SDEs. This paper will focus on stability issues and variable stepsize implementation techniques for numerically solving SDEs effectively. (C) 2000 Elsevier Science B.V. All rights reserved.

A bound on the maximum strong order of stochastic Runge-Kutta methods for stochastic ordinary differential equations

In Burrage and Burrage [1] it was shown that by introducing a very general formulation for stochastic Runge-Kutta methods, the previous strong order barrier of order one could be broken without having to use higher derivative terms. In particular, methods of strong order 1.5 were developed in which a Stratonovich integral of order one and one of order two were present in the formulation. In this present paper, general order results are proven about the maximum attainable strong order of these stochastic Runge-Kutta methods (SRKs) in terms of the order of the Stratonovich integrals appearing in the Runge-Kutta formulation. In particular, it will be shown that if an s-stage SRK contains Stratonovich integrals up to order p then the strong order of the SRK cannot exceed min{(p + 1)/2, (s - 1)/2), p greater than or equal to 2, s greater than or equal to 3 or 1 if p = 1.

Efficient solution of multiple cracks in great number using eigen COD boundary integral equations with iteration procedure

A newly developed computational approach is proposed in the paper for the analysis of multiple crack problems based on the eigen crack opening displacement (COD) boundary integral equations. The eigen COD particularly refers to a crack in an infinite domain under fictitious traction acting on the crack surface. With the concept of eigen COD, the multiple cracks in great number can be solved by using the conventional displacement discontinuity boundary integral equations in an iterative fashion with a small size of system matrix to determine all the unknown CODs step by step. To deal with the interactions among cracks for multiple crack problems, all cracks in the problem are divided into two groups, namely the adjacent group and the far-field group, according to the distance to the current crack in consideration. The adjacent group contains cracks with relatively small distances but strong effects to the current crack, while the others, the cracks of far-field group are composed of those with relatively large distances. Correspondingly, the eigen COD of the current crack is computed in two parts. The first part is computed by using the fictitious tractions of adjacent cracks via the local Eshelby matrix derived from the traction boundary integral equations in discretized form, while the second part is computed by using those of far-field cracks so that the high computational efficiency can be achieved in the proposed approach. The numerical results of the proposed approach are compared not only with those using the dual boundary integral equations (D-BIE) and the BIE with numerical Green's functions (NGF) but also with those of the analytical solutions in literature. The effectiveness and the efficiency of the proposed approach is verified. Numerical examples are provided for the stress intensity factors of cracks, up to several thousands in number, in both the finite and infinite plates.

Efficient solution of two-sided nonlinear space-fractional diffusion equations using fast Poisson preconditioners

We develop a fast Poisson preconditioner for the efficient numerical solution of a class of two-sided nonlinear space fractional diffusion equations in one and two dimensions using the method of lines. Using the shifted Gr¨unwald finite difference formulas to approximate the two-sided(i.e. the left and right Riemann-Liouville) fractional derivatives, the resulting semi-discrete nonlinear systems have dense Jacobian matrices owing to the non-local property of fractional derivatives. We employ a modern initial value problem solver utilising backward differentiation formulas and Jacobian-free Newton-Krylov methods to solve these systems. For efficient performance of the Jacobianfree Newton-Krylov method it is essential to apply an effective preconditioner to accelerate the convergence of the linear iterative solver. The key contribution of our work is to generalise the fast Poisson preconditioner, widely used for integer-order diffusion equations, so that it applies to the two-sided space fractional diffusion equation. A number of numerical experiments are presented to demonstrate the effectiveness of the preconditioner and the overall solution strategy.

A preconditioned Lanczos method for space-fractional reaction-diffusion equations on two-dimensional unstructured meshes

We consider a two-dimensional space-fractional reaction diffusion equation with a fractional Laplacian operator and homogeneous Neumann boundary conditions. The finite volume method is used with the matrix transfer technique of Ilić et al. (2006) to discretise in space, yielding a system of equations that requires the action of a matrix function to solve at each timestep. Rather than form this matrix function explicitly, we use Krylov subspace techniques to approximate the action of this matrix function. Specifically, we apply the Lanczos method, after a suitable transformation of the problem to recover symmetry. To improve the convergence of this method, we utilise a preconditioner that deflates the smallest eigenvalues from the spectrum. We demonstrate the efficiency of our approach for a fractional Fisher’s equation on the unit disk.

Sediment-transport (wind) experiments in zero gravity

The carousel wind tunnel (CWT) can be a significant tool for the determination of the nature and magnitude of interparticlar forces at threshold of motion. By altering particle and drum surface electrical properties and/or by applying electric potential difference across the inner and outer drums, it should be possible to separate electrostatic effects from other forces of cohesion. Besides particle trajectory and bedform analyses, suggestions for research include particle aggregation in zero and sub-gravity environments, effect of suspension-saltation ratio on soil abrasion, and the effects of shear and shear free turbulence on particle aggregation as applied to evolution of solar nebula.

Gravity-driven fingering simulations for a thin liquid film flowing down the outside of a vertical cylinder

A numerical study is presented to examine the fingering instability of a gravity-driven thin liquid film flowing down the outer wall of a vertical cylinder. The lubrication approximation is employed to derive an evolution equation for the height of the film, which is dependent on a single parameter, the dimensionless cylinder radius. This equation is identified as a special case of that which describes thin film flow down an inclined plane. Fully three-dimensional simulations of the film depict a fingering pattern at the advancing contact line. We find the number of fingers observed in our simulations to be in excellent agreement with experimental observations and a linear stability analysis reported recently by Smolka & SeGall (Phys Fluids 23, 092103 (2011)). As the radius of the cylinder decreases, the modes of perturbation have an increased growth rate, thus increasing cylinder curvature partially acts to encourage the contact line instability. In direct competition with this behaviour, a decrease in cylinder radius means that fewer fingers are able to form around the circumference of the cylinder. Indeed, for a sufficiently small radius, a transition is observed, at which point the contact line is stable to transverse perturbations of all wavenumbers. In this regime, free surface instabilities lead to the development of wave patterns in the axial direction, and the flow features become perfectly analogous to the two-dimensional flow of a thin film down an inverted plane as studied by Lin & Kondic (Phys Fluids 22, 052105 (2010)). Finally, we simulate the flow of a single drop down the outside of the cylinder. Our results show that for drops with low volume, the cylinder curvature has the effect of increasing drop speed and hence promoting the phenomenon of pearling. On the other hand, drops with much larger volume evolve to form single long rivulets with a similar shape to a finger formed in the aforementioned simulations.

Numerical solutions for thin film flow down the outside and inside of a vertical cylinder

We consider a model for thin film flow down the outside and inside of a vertical cylinder. Our focus is to study the effect that the curvature of the cylinder has on the gravity-driven instability of the advancing contact line and to simulate the resulting fingering patterns that form due to this instability. The governing partial differential equation is fourth order with a nonlinear degenerate diffusion term that represents the stabilising effect of surface tension. We present numerical solutions obtained by implementing an efficient alternating direction implicit scheme. When compared to the problem of flow down a vertical plane, we find that increasing substrate curvature tends to increase the fingering instability for flow down the outside of the cylinder, whereas flow down the inside of the cylinder substrate curvature has the opposite effect. Further, we demonstrate the existence of nontrivial travelling wave solutions which describe fingering patterns that propagate down the inside of a cylinder at constant speed without changing form. These solutions are perfectly analogous to those found previously for thin film flow down an inclined plane.

Synthesis and characterisation of titanium sol-gels in varied gravity environments

Sol-gel synthesis in varied gravity is only a relatively new topic in the literature and further investigation is required to explore its full potential as a method to synthesise novel materials. Although trialled for systems such as silica, the specific application of varied gravity synthesis to other sol-gel systems such as titanium has not previously been undertaken. Current literature methods for the synthesis of sol-gel material in reduced gravity could not be applied to titanium sol-gel processing, thus a new strategy had to be developed in this study. To successfully conduct experiments in varied gravity a refined titanium sol-gel chemical precursor had to be developed which allowed the single solution precursor to remain un-reactive at temperatures up to 50oC and only begin to react when exposed to a pressure decrease from a vacuum. Due to the new nature of this precursor, a thorough characterisation of the reaction precursors was subsequently undertaken with the use of techniques such as Nuclear Magnetic Resonance, Infra-red and UV-Vis spectroscopy in order to achieve sufficient understanding of precursor chemistry and kinetic stability. This understanding was then used to propose gelation reaction mechanisms under varied gravity conditions. Two unique reactor systems were designed and built with the specific purpose to allow the effects of varied gravity (high, normal, reduced) during synthesis of titanium sol-gels to be studied. The first system was a centrifuge capable of providing high gravity environments of up to 70 g’s for extended periods, whilst applying a 100 mbar vacuum and a temperature of 40-50oC to the reaction chambers. The second system to be used in the QUT Microgravity Drop Tower Facility was also required to provide the same thermal and vacuum conditions used in the centrifuge, but had to operate autonomously during free fall. Through the use of post synthesis characterisation techniques such as Raman Spectroscopy, X-Ray diffraction (XRD) and N2 adsorption, it was found that increased gravity levels during synthesis, had the greatest effect on the final products. Samples produced in reduced and normal gravity appeared to form amorphous gels containing very small particles with moderate surface areas. Whereas crystalline anatase (TiO2), was found to form in samples synthesised above 5 g with significant increases in crystallinity, particle size and surface area observed when samples were produced at gravity levels up to 70 g. It is proposed that for samples produced in higher gravity, an increased concentration gradient of water is forms at the bottom of the reacting film due to forced convection. The particles formed in higher gravity diffuse downward towards this excess of water, which favours the condensation reaction of remaining sol gel precursors with the particles promoting increased particle growth. Due to the removal of downward convection in reduced gravity, particle growth due to condensation reaction processes are physically hindered hydrolysis reactions favoured instead. Another significant finding from this work was that anatase could be produced at relatively low temperatures of 40-50oC instead of the conventional method of calcination above 450oC solely through sol-gel synthesis at higher gravity levels. It is hoped that the outcomes of this research will lead to an increased understanding of the effects of gravity on chemical synthesis of titanium sol-gel, potentially leading to the development of improved products suitable for diverse applications such as semiconductor or catalyst materials as well as significantly reducing production and energy costs through manufacturing these materials at significantly lower temperatures.

Segmental torso masses and gravity-induced coronal plane joint moments in adolescent idiopathic scoliosis

Introduction: Calculating segmental (vertebral level-by-level) torso masses in Adolescent Idiopathic Scoliosis (AIS) patients allows the gravitational loading on the scoliotic spine during relaxed standing to be estimated. This study used supine CT scans of AIS patients to measure segmental torso masses and explored the joint moments in the coronal plane, particularly at the apex of a scoliotic major curve. Methods: Existing low dose CT data from the Paediatric Spine Research Group was used to calculate vertebral level-by-level torso masses and joint moments occurring in the spine for a group of 20 female AIS patients with right sided thoracic curves. The mean age was 15.0 ± 2.7 years and all curves were classified Lenke Type 1 with a mean Cobb angle 52 ± 5.9°. Image processing software, ImageJ (v1.45 NIH USA) was used to create reformatted coronal plane images, reconstruct vertebral level-by-level torso segments and subsequently measure the torso volume corresponding to each vertebral level. Segment mass was then determined by assuming a tissue density of 1.04x103 kg/m3. Body segment masses for the head, neck and arms were taken from published anthropometric data (Winter 2009). Intervertebral joint moments in the coronal plane at each vertebral level were found from the position of the centroid of the segment masses relative to the joint centres with the segmental body mass data. Results and Discussion: The magnitude of the torso masses from T1-L5 increased inferiorly, with a 150% increase in mean segmental torso mass from 0.6kg at T1 to 1.5kg at L5. The magnitudes of the calculated coronal plane joint moments during relaxed standing were typically 5-7 Nm at the apex of the curve, with the highest apex joint torque of 7Nm. The CT scans were performed in the supine position and curve magnitudes are known to be 7-10° smaller than those measured in standing, due to the absence of gravity acting on the spine. Hence, it can be expected that the moments produced by gravity in the standing individual will be greater than those calculated here.

Morphology of osteocyte cells in normal and hyper-gravity

Osteocyte cells are the most abundant cells in human bone tissue. Due to their unique morphology and location, osteocyte cells are thought to act as regulators in the bone remodelling process, and are believed to play an important role in astronauts’ bone mass loss after long-term space missions. There is increasing evidence showing that an osteocyte’s functions are highly affected by its morphology. However, changes in an osteocyte’s morphology under an altered gravity environment are still not well documented. Several in vitro studies have been recently conducted to investigate the morphological response of osteocyte cells to the microgravity environment, where osteocyte cells were cultured on a two-dimensional flat surface for at least 24 hours before microgravity experiments. Morphology changes of osteocyte cells in microgravity were then studied by comparing the cell area to 1g control cells. However, osteocyte cells found in vivo are with a more 3D morphology, and both cell body and dendritic processes are found sensitive to mechanical loadings. A round shape osteocyte’s cells support a less stiff cytoskeleton and are more sensitive to mechanical stimulations compared with flat cellular morphology. Thus, the relative flat and spread shape of isolated osteocytes in 2D culture may greatly hamper their sensitivity to a mechanical stimulus, and the lack of knowledge on the osteocyte’s morphological characteristics in culture may lead to subjective and noncomprehensive conclusions of how altered gravity impacts on an osteocyte’s morphology. Through this work empirical models were developed to quantitatively predicate the changes of morphology in osteocyte cell lines (MLO-Y4) in culture, and the response of osteocyte cells, which are relatively round in shape, to hyper-gravity stimulation has also been investigated. The morphology changes of MLO-Y4 cells in culture were quantified by measuring cell area and three dimensionless shape features including aspect ratio, circularity and solidity by using widely accepted image analysis software (ImageJTM). MLO-Y4 cells were cultured at low density (5×103 per well) and the changes in morphology were recorded over 10 hours. Based on the data obtained from the imaging analysis, empirical models were developed using the non-linear regression method. The developed empirical models accurately predict the morphology of MLO-Y4 cells for different culture times and can, therefore, be used as a reference model for analysing MLO-Y4 cell morphology changes within various biological/mechanical studies, as necessary. The morphological response of MLO-Y4 cells with a relatively round morphology to hyper-gravity environment has been investigated using a centrifuge. After 2 hours culture, MLO-Y4 cells were exposed to 20g for 30mins. Changes in the morphology of MLO-Y4 cells are quantitatively analysed by measuring the average value of cell area and dimensionless shape factors such as aspect ratio, solidity and circularity. In this study, no significant morphology changes were detected in MLO-Y4 cells under a hyper-gravity environment (20g for 30 mins) compared with 1g control cells.

A simplified approach for calculating the moments of action for linear reaction-diffusion equations

The mean action time is the mean of a probability density function that can be interpreted as a critical time, which is a finite estimate of the time taken for the transient solution of a reaction-diffusion equation to effectively reach steady state. For high-variance distributions, the mean action time under-approximates the critical time since it neglects to account for the spread about the mean. We can improve our estimate of the critical time by calculating the higher moments of the probability density function, called the moments of action, which provide additional information regarding the spread about the mean. Existing methods for calculating the nth moment of action require the solution of n nonhomogeneous boundary value problems which can be difficult and tedious to solve exactly. Here we present a simplified approach using Laplace transforms which allows us to calculate the nth moment of action without solving this family of boundary value problems and also without solving for the transient solution of the underlying reaction-diffusion problem. We demonstrate the generality of our method by calculating exact expressions for the moments of action for three problems from the biophysics literature. While the first problem we consider can be solved using existing methods, the second problem, which is readily solved using our approach, is intractable using previous techniques. The third problem illustrates how the Laplace transform approach can be used to study coupled linear reaction-diffusion equations.

Structure-preserving Runge-Kutta methods for stochastic Hamiltonian equations with additive noise

There has been considerable recent work on the development of energy conserving one-step methods that are not symplectic. Here we extend these ideas to stochastic Hamiltonian problems with additive noise and show that there are classes of Runge-Kutta methods that are very effective in preserving the expectation of the Hamiltonian, but care has to be taken in how the Wiener increments are sampled at each timestep. Some numerical simulations illustrate the performance of these methods.

Development and validation of anthropometric prediction equations for estimation of lean body mass and appendicular lean soft tissue in Indian men and women

Lean body mass (LBM) and muscle mass remains difficult to quantify in large epidemiological studies due to non-availability of inexpensive methods. We therefore developed anthropometric prediction equations to estimate the LBM and appendicular lean soft tissue (ALST) using dual energy X-ray absorptiometry (DXA) as a reference method. Healthy volunteers (n= 2220; 36% females; age 18-79 y) representing a wide range of body mass index (14-44 kg/m2) participated in this study. Their LBM including ALST was assessed by DXA along with anthropometric measurements. The sample was divided into prediction (60%) and validation (40%) sets. In the prediction set, a number of prediction models were constructed using DXA measured LBM and ALST estimates as dependent variables and a combination of anthropometric indices as independent variables. These equations were cross-validated in the validation set. Simple equations using age, height and weight explained > 90% variation in the LBM and ALST in both men and women. Additional variables (hip and limb circumferences and sum of SFTs) increased the explained variation by 5-8% in the fully adjusted models predicting LBM and ALST. More complex equations using all the above anthropometric variables could predict the DXA measured LBM and ALST accurately as indicated by low standard error of the estimate (LBM: 1.47 kg and 1.63 kg for men and women, respectively) as well as good agreement by Bland Altman analyses. These equations could be a valuable tool in large epidemiological studies assessing these body compartments in Indians and other population groups with similar body composition.