Integration Of Ordinary Differential Equations Research Articles

Exploring Inductive Linearization for simulation and estimation with an application to the Michaelis–Menten model

Nonlinear ordinary differential equations (ODEs) are common in pharmacokinetic–pharmacodynamic systems. Although their exact solutions cannot generally be determined via algebraic methods, their rapid and accurate solutions are desirable. Thus, numerical methods have a critical role. Inductive Linearization was proposed as a method to solve systems of nonlinear ODEs. It is an iterative approach that converts a nonlinear ODE into a linear time-varying (LTV) ODE, for which a range of standard integration techniques can then be used to solve (e.g., eigenvalue decomposition [EVD]). This study explores the properties of Inductive Linearization when coupled with EVD for integration of the LTV ODE and illustrates how the efficiency of the method can be improved. Improvements were based on three approaches, (1) incorporation of a convergence criterion for the iterative linearization process (for simulation and estimation), (2) creating more efficient step sizes for EVD (for simulation and estimation), and (3) updating the initial conditions of the Inductive Linearization (for estimation). The performance of these improvements were evaluated using single subject stochastic simulation-estimation with an application to a simple pharmacokinetic model with Michaelis–Menten elimination. The reference comparison was a standard non-stiff Runge–Kutta method with variable step size (ode45, MATLAB). Each of the approaches improved the speed of the Inductive Linearization technique without diminishing accuracy which, in this simple case, was faster than ode45 with comparable accuracy in the parameter estimates. The methods described here can easily be implemented in standard software programme such as R or MATLAB. Further work is needed to explore this technique for estimation in a population approach setting.

An Algorithm for Numerical Integration of ODE with Sampled Unknown Functional Factors

The problem of having ordinary differential equations (ODE) whose coefficients are unknown functions is frequent in several fields. Sometimes, it is possible to obtain samples of the values of these functions in different instants or spatial points. The present paper presents a methodology for the numeric solving of these ODE. There are approximations to the problem for specific cases of equations, especially in the case where the parameters correspond to constants. Other studies focus on the case in which the functions under consideration are linear or meet a certain condition. There are two main advantages of the proposed algorithm. First, it does not impose any condition over the data or the subsequent function from where these sample data are derived. Additionally, the methodology used in the functions modeling can control the possibility of overfitting in the function modeling. This is a crucial point in order to limit the influence of model biases in the numerical solution of the ordinary differential equation under study.

A numerical verification method to specify homoclinic orbits as application of local Lyapunov functions

AbstractWe propose a verification method for specification of homoclinic orbits as application of our previous work for constructing local Lyapunov functions by verified numerics. Our goal is to specify parameters appeared in the given systems of ordinary differential equations (ODEs) which admit homoclinic orbits to equilibria. Here we restrict ourselves to cases that each equilibrium is independent of parameters. The feature of our methods consists of Lyapunov functions, integration of ODEs by verified numerics, and Brouwer’s coincidence theorem on continuous mappings. Several techniques for constructing continuous mappings from a domain of parameter vectors to a region of the phase space are shown. We present numerical examples for problems in 3 and 4-dimensional cases.

Implementation of a practical Markov chain Monte Carlo sampling algorithm in PyBioNetFit.

Bayesian inference in biological modeling commonly relies on Markov chain Monte Carlo (MCMC) sampling of a multidimensional and non-Gaussian posterior distribution that is not analytically tractable. Here, we present the implementation of a practical MCMC method in the open-source software package PyBioNetFit (PyBNF), which is designed to support parameterization of mathematical models for biological systems. The new MCMC method, am, incorporates an adaptive move proposal distribution. For warm starts, sampling can be initiated at a specified location in parameter space and with a multivariate Gaussian proposal distribution defined initially by a specified covariance matrix. Multiple chains can be generated in parallel using a computer cluster. We demonstrate that am can be used to successfully solve real-world Bayesian inference problems, including forecasting of new Coronavirus Disease 2019 case detection with Bayesian quantification of forecast uncertainty. PyBNF version 1.1.9, the first stable release with am, is available at PyPI and can be installed using the pip package-management system on platforms that have a working installation of Python 3. PyBNF relies on libRoadRunner and BioNetGen for simulations (e.g. numerical integration of ordinary differential equations defined in SBML or BNGL files) and Dask.Distributed for task scheduling on Linux computer clusters. The Python source code can be freely downloaded/cloned from GitHub and used and modified under terms of the BSD-3 license (https://github.com/lanl/pybnf). Online documentation covering installation/usage is available (https://pybnf.readthedocs.io/en/latest/). A tutorial video is available on YouTube (https://www.youtube.com/watch?v=2aRqpqFOiS4&t=63s). Supplementary data are available at Bioinformatics online.

АСИМПТОТИКА РЕШЕНИЯ ПЕРВОЙ КРАЕВОЙ ЗАДАЧИ ДЛЯ СИНГУЛЯРНО ВОЗМУЩЕННОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ В ЧАСТНЫХ ПРОИЗВОДНЫХ ВТОРОГО ПОРЯДКА ПАРАБОЛИЧЕСКОГО ТИПА

The article constructs a complete uniform asymptotic expansion in a small parameter of the solution of the first boundary value problem. The first boundary value problem is posed for a singularly perturbed linear inhomogeneous second-order partial differential equation with two independent variables of parabolic type. The problem is investigated on a rectangle. The peculiarities of the problem are the presence of a small parameter in front of the heat conduction operator, the existence of corner boundary layers at the lower corners of the rectangle. It is required to construct a uniform asymptotic expansion of the solution of the first boundary value problem on a rectangle, with any degree of accuracy, as the small parameter tends to zero. The asymptotic expansion of the solution in terms of a small parameter is constructed by the Vishik–Lyusternik method. When solving the problem, we use: methods of integration of ordinary differential equations, the classical method of a small parameter, the Vishik–Lyusternik boundary function method, and the maximum principle. As usual, the problem is solved in two stages: in the first stage, a formal expansion of the solution of the first boundary value problem is constructed; and in the second stage, the remainder of the resulting expansion is estimated and this proves that the resulting expansion is indeed asymptotic over the entire rectangle. In the first stage, a formal asymptotic solution is sought in the form of a sum of six functions (solutions): an external solution defined on the entire rectangle; boundary layer solution in a small neighborhood of the lower side of the rectangle; two lateral boundary layer solutions in a small neighborhood of the lateral sides of the rectangle and two corner boundary layer solutions in the neighborhood of the lower vertices of the rectangle. All these boundary layer solutions exponentially decrease outside the boundary layers.

Stable linear multistep methods with off-step points for the solution of ordinary differential equations

Of recent, stability has become an important concept and a qualitative property in any numerical integration scheme. In this work, we propose two stable linear multistep methods with off-step points for the numerical integration of ordinary differential equations whose development is collocation and interpolation based. The boundary locus techniques show that the proposed schemes are zero-stable, A-stable and -stable for some step number and are found suitable for stiff differential equations. Numerical results obtained compare favourably with some existing methods in literature.

The Double Phospho/Dephosphorylation Cycle as a Benchmark to Validate an Effective Taylor Series Method to Integrate Ordinary Differential Equations

The double phosphorylation/dephosphorylation cycle consists of a symmetric network of biochemical reactions of paramount importance in many intracellular mechanisms. From a network perspective, they consist of four enzymatic reactions interconnected in a specular way. The general approach to model enzymatic reactions in a deterministic fashion is by means of stiff Ordinary Differential Equations (ODEs) that are usually hard to integrate according to biologically meaningful parameter settings. Indeed, the quest for model simplification started more than one century ago with the seminal works by Michaelis and Menten, and their Quasi Steady-State Approximation methods are still matter of investigation nowadays. This work proposes an effective algorithm based on Taylor series methods that manages to overcome the problems arising in the integration of stiff ODEs, without settling for model approximations. The double phosphorylation/dephosphorylation cycle is exploited as a benchmark to validate the methodology from a numerical viewpoint.

Learning Dynamic Systems Using Gaussian Process Regression with Analytic Ordinary Differential Equations as Prior Information

Recently the data-driven learning of dynamic systems has become a promising approach because no physical knowledge is needed. Pure machine learning approaches such as Gaussian process regression (GPR) learns a dynamic model from data, with all physical knowledge about the system discarded. This goes from one extreme, namely methods based on optimizing parametric physical models derived from physical laws, to the other. GPR has high flexibility and is able to model any dynamics as long as they are locally smooth, but can not generalize well to unexplored areas with little or no training data. The analytic physical model derived under assumptions is an abstract approximation of the true system, but has global generalization ability. Hence the optimal learning strategy is to combine GPR with the analytic physical model. This paper proposes a method to learn dynamic systems using GPR with analytic ordinary differential equations (ODEs) as prior information. The one-time-step integration of analytic ODEs is used as the mean function of the Gaussian process prior. The total parameters to be trained include physical parameters of analytic ODEs and parameters of GPR. A novel method is proposed to simultaneously learn all parameters, which is realized by the fully Bayesian GPR and more promising to learn an optimal model. The standard Gaussian process regression, the ODE method and the existing method in the literature are chosen as baselines to verify the benefit of the proposed method. The predictive performance is evaluated by both one-time-step prediction and long-term prediction. By simulation of the cart-pole system, it is demonstrated that the proposed method has better predictive performances.

DLBFoam: An open-source dynamic load balancing model for fast reacting flow simulations in OpenFOAM

Computational load imbalance is a well-known performance issue in multiprocessor reacting flow simulations utilizing directly integrated chemical kinetics. We introduce an open-source dynamic load balancing model named DLBFoam to address this issue within OpenFOAM, an open-source C++ library for Computational Fluid Dynamics (CFD). Due to the commonly applied operator splitting practice in reactive flow solvers, chemistry can be treated as an independent stiff ordinary differential equation (ODE) system within each computational cell. As a result of the highly non-linear characteristics of chemical kinetics, a large variation in the convergence rates of the ODE integrator may occur, leading to a high load imbalance across multiprocessor configurations. However, the independent nature of chemistry ODE systems leads to a problem that can be parallelized easily (called an embarrassingly parallel problem in the literature) during the flow solution. The presented model takes advantage of this feature and balances the chemistry load across available resources. Additionally, a reference mapping model is utilized to further speed-up the simulations. When DLBFoam it utilized with both these features enabled, a speed-up by a factor of 10 is reported for reactive flow benchmark cases. To the best of our knowledge, this model is the first open-source implementation of chemistry load balancing in the literature. Program summaryProgram Title: DLBFoamCPC Library link to program files:https://doi.org/10.17632/bb9zjfzcmm.1Developer's repository link:https://github.com/blttkgl/DLBFoamLicensing provisions: GPLv3Programming language: C++Nature of problem: Solution of chemical kinetics in parallel reacting flow solvers raises a computational imbalance across multiprocessor architectures. DLBFoam balances the load distribution evenly, providing significant speed-up in reacting CFD applications.Solution method: The dynamic load balancing is implemented by distributing the point-wise chemistry problems from most loaded processes to less loaded ones using MPI communication protocol.Additional comments including restrictions and unusual features: The present model is designed to work with the standard chemistry model class available in OpenFOAM (versions 7 and 8). For the time being, the model does not support derived combustion models such as “TDAC” and covers gas-phase reaction kinetics only. In addition, the boundary surface chemistry problems are neglected by the model.

Applying the Method of Integration of Ordinary Differential Equations Based on the Chebyshev Series to the Restricted Plane Circular Three-Body Problem

Applying the Method of Integration of Ordinary Differential Equations Based on the Chebyshev Series to the Restricted Plane Circular Three-Body Problem

Some Remarks on the Shifner-Erougin-Salakhova-Chebotarev Type Differential Equations

This article studies a type of integrable linear ordinary differential equations considered by Shifner [17] and generalised by Erougin [5], Salakhova and Chebotarev [15]. We will consider a slight generalisation of their results.

Enhancing Accuracy of Runge–Kutta-Type Collocation Methods for Solving ODEs

In this paper, a new class of Runge–Kutta-type collocation methods for the numerical integration of ordinary differential equations (ODEs) is presented. Its derivation is based on the integral form of the differential equation. The approach enables enhancing the accuracy of the established collocation Runge–Kutta methods while retaining the same number of stages. We demonstrate that, with the proposed approach, the Gauss–Legendre and Lobatto IIIA methods can be derived and that their accuracy can be improved for the same number of method coefficients. We expressed the methods in the form of tables similar to Butcher tableaus. The performance of the new methods is investigated on some well-known stiff, oscillatory, and nonlinear ODEs from the literature.

Common dynamic estimation via structured low-rank approximation with multiple rank constraints

Abstract We consider the problem of detecting the common dynamic among several observed signals. It has been shown in (Markovsky et al., 2019) that the problem is equivalent to a generalization of the classical Hankel low-rank approximation to the case of multiple rank constraints. We propose an optimization method based on the integration of ordinary differential equations describing a descent dynamic for a suitable functional to be minimized. We show how the proposed algorithm improves the numerical solutions computed by existing subspace methods which solve the same problem.

Superconvergence of Galerkin variational integrators

We study the order of convergence of Galerkin variational integrators for ordinary differential equations. Galerkin variational integrators approximate a variational (Lagrangian) problem by restricting the space of curves to the set of polynomials of degree at most $s$ and approximating the action integral using a quadrature rule. We show that, if the quadrature rule is sufficiently accurate, the order of the integrators thus obtained is $2s$.

АСИМПТОТИКА РЕШЕНИЯ ДВУХЗОННОЙ ДВУХТОЧЕЧНОЙ КРАЕВОЙ ЗАДАЧИ

The article investigates the asymptotic behavior of the solution of a two-point boundary value problem on an interval for a linear inhomogeneous ordinary differential equation of the second order with a small parameter at the highest derivative. The essential features of the problem are the presence of a small parameter in front of the second-order derivative of the desired function, the existence of a twodimensional boundary layer at the left end of the segment at x = 0, and the non-smoothness of the solution to the corresponding unperturbed boundary value problem. It is required to construct a uniform asymptotic expansion of the solution to a two-zone two-point boundary value problem on a unit interval, with any degree of accuracy, as the small parameter tends to zero. Due to the second and third features of the problem, it is not easy to construct an asymptotic solution expansion with respect to the small parameter using the known asymptotic methods. When solving the problem, the following methods are used: methods of integration of ordinary differential equations; the method of a small parameter; the classical method of boundary functions; and the generalized method of boundary functions and the maximum principle. The problem is solved in two stages: in the first stage, a formal expansion of the solution to the two-point boundary value problem is constructed, and in the second stage, the justification of this expansion is given, i.e. the remainder term of the expansion is estimated. In the first stage, a formal asymptotic solution is sought in the form of the sum of three solutions: a smooth outer solution on the entire segment; classical boundary layer solution in the vicinity of x = 0, which exponentially decreases outside the boundary layer; and an intermediate boundary layer solution at x = 0, which decreases in power mode outside the boundary layer. The constructed asymptotic expansion of the solution to the two-point boundary value problem is asymptotic in the sense of Erdei.

Integrable Abel equation and asymptotics of symmetry solutions of Korteweg-de Vries equation

We provide a general solution for a first order ordinary differential equation with a rational right-hand side, which arises in constructing asymptotics for large time of simultaneous solutions of the Korteweg-de Vries equation and the stationary part of its higher non-autonomous symmetry. This symmetry is determined by a linear combination of the first higher autonomous symmetry of the Korteweg-de Vries equation and of its classical Galileo symmetry. This general solution depends on an arbitrary parameter. By the implicit function theorem, locally it is determined by the first integral explicitly written in terms of hypergeometric functions. A particular case of the general solution defines self-similar solutions of the Whitham equations, found earlier by G.V. Potemin in 1988. In the well-known works by A.V. Gurevich and L.P. Pitaevsky in early 1970s, it was established that these solutions of the Whitham equations describe the origination in the leading term of non-damping oscillating waves in a wide range of problems with a small dispersion. The result of this article supports once again an empirical rule saying that under various passages to the limits, integrable equations can produce only integrable, in certain sense, equations. We propose a general conjecture: integrable ordinary differential equations similar to that considered in the present paper should also arise in describing the asymptotics at large times for other symmetry solutions to evolution equations admitting the application of the method of inverse scattering problem.

Numerical integrators for Lagrangian oceanography

Abstract. A common task in Lagrangian oceanography is to calculate a large number of drifter trajectories from a velocity field precalculated with an ocean model. Mathematically, this is simply numerical integration of an ordinary differential equation (ODE), for which a wide range of different methods exist. However, the discrete nature of the modelled ocean currents requires interpolation of the velocity field in both space and time, and the choice of interpolation scheme has implications for the accuracy and efficiency of the different numerical ODE methods. We investigate trajectory calculation in modelled ocean currents with 800 m, 4 km, and 20 km horizontal resolution, in combination with linear, cubic and quintic spline interpolation. We use fixed-step Runge–Kutta integrators of orders 1–4, as well as three variable-step Runge–Kutta methods (Bogacki–Shampine 3(2), Dormand–Prince 5(4) and 8(7)). Additionally, we design and test modified special-purpose variants of the three variable-step integrators, which are better able to handle discontinuous derivatives in an interpolated velocity field. Our results show that the optimal choice of ODE integrator depends on the resolution of the ocean model, the degree of interpolation, and the desired accuracy. For cubic interpolation, the commonly used Dormand–Prince 5(4) is rarely the most efficient choice. We find that in many cases, our special-purpose integrators can improve accuracy by many orders of magnitude over their standard counterparts, with no increase in computational effort. Equivalently, the special-purpose integrators can provide the same accuracy as standard methods at a reduced computational cost. The best results are seen for coarser resolutions (4 and 20 km), thus the special-purpose integrators are particularly advantageous for research using regional to global ocean models to compute large numbers of trajectories. Our results are also applicable to trajectory computations on data from atmospheric models.

A note on fitting a generalised Moody diagram for wall modelled large-eddy simulations

ABSTRACT Motivated by the needs of wall modelled Large Eddy Simulation (LES), we introduce fits to numerical solutions of the Reynolds Averaged Navier–Stokes equations in their simplest near-wall, boundary layer approximation including a mixing-length model. We formulate the problem such that independent dimensionless variables are those directly available in LES. We provide practical fits for the dependent variable, fits that encompass a smooth transition between the viscous sublayer and inertial logarithmic layer, and then progress first considering moderate pressure gradients as well as roughness effects under the assumption that the mixing-length is not affected by the pressure gradient. An alternative fit based on the empirical wall model of Nickels (Inner scaling for wall-bounded flows subject to large pressure gradients. J Fluid Mech. 2004;521:217–239) is also provided, taking into account possible effects of pressure gradient on turbulence near-wall structure. We then consider the case of general pressure gradients, both favourable and adverse, up to conditions of separation, for both smooth and rough surfaces. The proposed fitting functions constitute a generalised Moody chart, comply with analytical solutions valid in various asymptotic regimes, and obviate the need for numerical iterative solution methods or near-wall numerical integration of ordinary differential equations during LES.

A tutorial on solving ordinary differential equations using Python and hybrid physics-informed neural network

We present a tutorial on how to directly implement integration of ordinary differential equations through recurrent neural networks using Python. In order to simplify the implementation, we leveraged modern machine learning frameworks such as TensorFlow and Keras. Besides, offering implementation of basic models (such as multilayer perceptrons and recurrent neural networks) and optimization methods, these frameworks offer powerful automatic differentiation. With all that, the main advantage of our approach is that one can implement hybrid models combining physics-informed and data-driven kernels, where data-driven kernels are used to reduce the gap between predictions and observations. Alternatively, we can also perform model parameter identification. In order to illustrate our approach, we used two case studies. The first one consisted of performing fatigue crack growth integration through Euler’s forward method using a hybrid model combining a data-driven stress intensity range model with a physics-based crack length increment model. The second case study consisted of performing model parameter identification of a dynamic two-degree-of-freedom system through Runge–Kutta integration. The examples presented here as well as source codes are all open-source under the GitHub repository https://github.com/PML-UCF/pinn_code_tutorial.

Nested Implicit Runge–Kutta Pairs of Gauss and Lobatto Types with Local and Global Error Controls for Stiff Ordinary Differential Equations

The problem of efficient global error estimation and control is studied in embedded nested implicit Runge–Kutta pairs of Gauss and Lobatto types as applied to stiff ordinary differential equations (ODEs). Stiff problems may arise in many areas of engineering, and their accurate numerical solution is an important issue of computational and applied mathematics. A cheap global error estimation technique designed recently for the mentioned Runge–Kutta pairs can severely overestimate the global error when applied to stiff ODEs and, hence, this reduces the efficiency of those solvers. In the present paper, we explain the cause of that error overestimation and show how to improve the mentioned computation techniques for stiff problems. Such modifications not only boost the efficiency of numerical integration of stiff ODEs, but also make the embedded nested implicit Runge–Kutta pairs with scaled modified local and global error controls superior to stiff built-in MATLAB ODE solvers with only local error control when applied to important test examples.