Ordinary Differential Equation Solver Research Articles

Solving Ordinary Differential Equation Using Parallel Fourth Order Runge-Kutta Method With Three Processors

This paper presents a derivation of the Runge-Kutta or fourth method with six stages suitable for parallel implementation. Development of a parallel model based on the sparsity structure of the fourth type Runge-Kutta which is divided into three processors. The calculation of the parallel computation model and the sequential model from the accurate side shows that the sequential model is better. However, generally, the parallel method will end the analytic solution by increasing the number of iterations. In terms of execution time, parallel method has advantages over sequential method.

Preliminary Flow Modeling by Hybrid Automata Alternating Continuous Reaction and Discrete Transit for Pharmacokinetics.

In traditional pharmacokinetic models, blood flow or liquid transit is often expressed as first-order kinetics. When the flow expression by first-order kinetics is used for dynamic simulation, the flow velocity illogically depends on the step size of a solver of ordinary differential equations. In this study, we propose flow modeling using hybrid automata that combine ordinary differential equations and recursive equations, and we have preliminarily applied the constructed models to several examples. The blood concentration-time profiles of p-aminohippurate and propranolol after intravenous administration were successfully reproduced by simple hybrid automata. The simulation results of one-dimensional tube flow have demonstrated that the fluid velocity in the hybrid automata was independent of the step size of the ordinary differential equation solver. A body fluid model coordinated various flows in a human body with scheduled daily activities and could be used as a drug container to describe formulation-dependent disposition of 5-aminosalicylic acid and enterohepatic circulation of a virtual drug. These findings suggested that flow modeling using hybrid automata could avoid the logical inconsistency in the traditional pharmacokinetic modeling and that the hybrid automata have high versatility and a wide range of applicability to pharmacokinetic analysis. Because our method can define various intervals for multiple recursive equations, the resolution of a specific part of a model can be adjusted relatively freely while the whole body is being roughly modeled, which would be beneficial to refine a coarse model into a fine model in future. SIGNIFICANCE STATEMENT: There is a logical inconsistency in flow expression by first-order kinetics in ordinary differential equations used in traditional pharmacokinetic modeling. It is difficult to model a whole human body using flow models in partial differential equations because of the excessive calculation costs. Our simulations on tube flow and body fluids have demonstrated that the flow modeling using hybrid automata could avoid the problems. The preliminary applications of hybrid automata to several examples highlighted its high versatility in pharmacokinetic analysis.

Discussion of Numerical Methods used in Positive Displacement Comprehensive Mechanistic Models: Case Study using the Z-Compressor

Comprehensive mechanistic models have been widely used for simulating positive displacement (PD) compressors or expanders. Such models are based on a thermodynamic analysis of generic control volumes that describe the working chambers. The thermodynamic analysis results in a set of non-linear ordinary differential equations (ODEs) that describe conservation of mass and energy. The derived ODEs are used to predict the evolution of the thermodynamic properties of the working fluid with respect to crank angle or time. The solution scheme of a mechanistic model involves different layers of solvers that handle the step-by-step integration, the cycle-to-cycle continuity as well as the overall energy balance. Since most of the PD machines feature multiple control volumes that merge or split (e.g. scroll-type) or wrap-around a rotor (e.g. screw-type) the integration scheme faces challenges due to the stiffness of the equations or discontinuities in the control volume definition.In this paper, a discussion of different integration schemes and numerical considerations are proposed that are applied to a novel compressor called the z-compressor which presents a high-degree of intrinsic stiffness. A simplified comprehensive model is developed which is validated against experimental mass flow rate data for a prototype compressor with an average model-predicted error of less than 2%. This model is used to explore the differences in two ODE solvers from the ODE suite in MATLAB, (ODE15s, and ODE113). These results are further compared against an existing, simplified, Z-compressor model developed in PDSim which uses an explicit, lower order, solver the Adaptive Runge-Kutta (RKF). The results suggest that the semi-implicit solver (ODE15s) provides an efficient means for stepping through the compression process ODEs compared to the explicit solvers.

A unified and rigorous solution for quasi-static cylindrical cavity expansion in plasticity constitutive models

This paper presents a unified and rigorous solution for the quasi-static cylindrical cavity expansion problem in all existing plasticity constitutive models. Rigorous three-dimensional definitions of p’ and q are used as proposed by Chen & Abousleiman. Then, the p-q stress space is transformed to the three stress components in a cylindrical coordinate system. The governing partial differential equations (PDEs) for cylindrical cavity expansion, including the stress equilibrium equation, the constitutive equation, the consistency equation, the continuity equation, and the drainage conditions, are written with respect to the three stress components in a cylindrical coordinate system. The PDEs are subsequently reduced to ordinary differential equations (ODEs), which can be summarized in unified matrix form by means of a similarity transformation. A rigorous semi-analytical solution can immediately be obtained by numerically solving the ODEs using commercial ODE solvers, such as MATLAB. The proposed solution procedure is general and can be applied to both the drained and undrained conditions in all plasticity constitutive models. Finally, the partially drained effect was discussed through numerical analysis, and empirical equations considering the partially drained effect were proposed for calculating the limit effective radial stress and excess pore pressure at the cavity wall.

From SIR to SEAIRD: A novel data-driven modeling approach based on the Grey-box System Theory to predict the dynamics of COVID-19.

Common compartmental modeling for COVID-19 is based on a priori knowledge and numerous assumptions. Additionally, they do not systematically incorporate asymptomatic cases. Our study aimed at providing a framework for data-driven approaches, by leveraging the strengths of the grey-box system theory or grey-box identification, known for its robustness in problem solving under partial, incomplete, or uncertain data. Empirical data on confirmed cases and deaths, extracted from an open source repository were used to develop the SEAIRD compartment model. Adjustments were made to fit current knowledge on the COVID-19 behavior. The model was implemented and solved using an Ordinary Differential Equation solver and an optimization tool. A cross-validation technique was applied, and the coefficient of determination R2 was computed in order to evaluate the goodness-of-fit of the model. Key epidemiological parameters were finally estimated and we provided the rationale for the construction of SEAIRD model. When applied to Brazil’s cases, SEAIRD produced an excellent agreement to the data, with an R2 ≥ 90%. The probability of COVID-19 transmission was generally high (≥ 95%). On the basis of a 20-day modeling data, the incidence rate of COVID-19 was as low as 3 infected cases per 100,000 exposed persons in Brazil and France. Within the same time frame, the fatality rate of COVID-19 was the highest in France (16.4%) followed by Brazil (6.9%), and the lowest in Russia (≤ 1%). SEAIRD represents an asset for modeling infectious diseases in their dynamical stable phase, especially for new viruses when pathophysiology knowledge is very limited.Supplementary InformationThe online version contains supplementary material available at 10.1007/s10489-021-02379-2.

A grid-free approach for simulating sweep and cyclic voltammetry.

We present a computational approach to simulate linear sweep and cyclic voltammetry experiments that does not require a discretized grid in space to quantify diffusion. By using a Green's function solution coupled to a standard implicit ordinary differential equation solver, we are able to simulate current and redox species concentrations using only a small grid in time. As a result, where benchmarking is possible, we find that the current method is faster than (and quantitatively identical to) established techniques. The present algorithm should help open the door for studying adsorption effects in inner sphere electrochemistry.

Least Square Methods Based on Control Points of Said Ball Curves for Solving Ordinary Differential Equations

This paper presents the use of Said Ball curve’s control points to approximate the solutions of linear ordinary differential equations (ODEs). Least squares methods (LSM) is proposed to find the control points of Said Ball curves by minimizing the error of residual function.The residual error is measured by taking the sum of squares of the Said Ball curve’s control points of the residual function. Then the approximate solution of ODEs is obtained by minimizing residual error.Two numerical examples are given in term of error and compared with the exact solution to demonstrate the efficiency of the proposed method.

Modeling the effect of application-specific program transformations on energy and performance improvements of parallel ODE solvers

Ordinary differential equations (ODEs) are important for modelling many problems from science and engineering and efficient ODE solvers are required, for example when solving time-dependent partial differential equations (PDEs) with the method of lines. Since an ODE solver may perform a large number of iteration steps, the execution time for solving an ODE problem might be quite large. Thus, a reduction of the execution time is desirable and should affect each iteration step of the simulation. Programming techniques to reduce the execution time of ODE solver are parallelism and modification of the memory access structure such that the memory access time decreases. In this article, we investigate multithreaded solution methods for ODEs with different memory access behavior and their influence on the performance. Additionally the energy consumption is considered. The parallelism is implemented as shared memory program for multicore processors. The memory access behavior is investigated using different program variants which result from application-specific program transformations changing the memory access order while guaranteeing the numerical correctness. For the investigation of the performance, experimental data have been gathered on five different recent multicore processors. Additionally, an analytical power and energy model for modeling the performance and energy consumption is introduced. As ODE solver, the popular embedded Runge-Kutta methods with error correction is used. The simulation problems are two different ODEs resulting from discretized PDEs. The experimental data give insight into the quite diverse performance behavior of the ODE solver variants solving the same problem on different platforms.

Comparison of various zero‐dimensional models to analyze the brachial artery blood flow of an arterio‐venous fistula during hemodialysis procedure

AbstractOne of the most recommended vascular accesses for hemodialysis is an autogenous arterio‐venous fistula (AVF), which is a surgically created conduit between an artery and a vein. AVFs can be compromised by stenosis‐induced arterial or venous thrombosis. An effective and efficient mathematical fluid mechanical model can provide insight into the mechanisms for stenosis and thrombosis. In an earlier paper, a laminar entry flow model (Shah) was applied to a fifth order lumped parameter zero dimensional (0D) model to predict the inlet flow rate in the brachial artery of an AVF. ODE45, a built in MATLAB ordinary differential equation (ODE) solver, was used to produce the solution. The Shah model showed closer correlation with a computational fluid dynamics model, performed in COMSOL, than did the Poiseuille flow model. In this study, an outlet flow rate in the brachial artery of an AVF is generated with third, fourth, and fifth order lumped parameter zero dimensional models with laminar and turbulent flow conditions within different arterial segments. Butcher's sixth order Runge‐Kutta (RK) method was used as an ODE solver and the results were compared with the MATLAB ODE solvers (ODE45 and ODE113). The lumped parameter model solved with Butcher's sixth order RK method generated less oscillation with a more stable and consistent outlet flow rate in comparison to the other solvers. Brachial artery blood flow, modeled with Shah's entry flow equation in the inlet segment and Blasius turbulent equation modeled at the outlet segment prior to the beginning of the anastomosis, generated an average outlet flow rate similar to that of clinical data. The dependence of brachial artery on the flow equations indicates the importance of wisely selecting flow resistance models and the effect of blood pressure and peripheral resistance on fistula maturation. The average and peak flow rate also depended on the ODE solver. Third and fourth order lumped parameter models generated results similar to the fifth order model only at normal peripheral resistance. Therefore, this study can help to understand the AVF performance in vitro and guide the application of an appropriate lumped parameter model for blood flow prediction.

Topology optimization employing a condensation method for nonlinear structural frames with supplemental mass

AbstractA topology optimization schema employing a condensation method for nonlinear frame structures with supplemental mass subjected to time‐varying excitation is presented. In the context of the design of structural frames, in certain applications, the supplemental mass can be order(s) of magnitude larger than the mass of the system itself. Thus, condensing the system of governing equations to only those associated with the supplemental mass, reduces the complexity and computational cost of the dynamic analysis and thus the optimization process. In addition to considering material nonlinearity, distributed plasticity, and multiaxial interactions, the hysteretic beam finite element (FE) model employed in this study has constant elastic stiffness and hysteretic matrices facilitating a Guyan‐type condensation of the nonlinear dynamic equations. Furthermore, the nonlinear dynamic equations and element evolution equations of hysteretic FE modeling are concisely presented as a system of first‐order nonlinear ordinary differential equations (ODEs) that can be solved using general ODE solvers. The schema is demonstrated on various design problems considering different arrangements of supplemental mass, where the objective is to minimize the volume of the structure subject to maximum displacement constraint using a single p‐norm.

A Numerical Method for Solving Ordinary Differential Equations by Converting Them into the Form of a Shannon

A Numerical Method for Solving Ordinary Differential Equations by Converting Them into the Form of a Shannon

Nelder-Mead Algorithm in Solving Ordinary Differential Equations Whose Solutions Possess Singularities

This research considers Initial Value Problems (IVPs) in Ordinary Differential Equations (ODEs) whose solutions possess singularities. Here, we represent the theoretical solution by a rational function as it is more convenient representing a function close to a singularity by a rational function. The process of transforming the IVP to a constrained optimization problem and application of Nelder-Mead algorithm in obtaining approximate solution is presented in this work. Accuracy and efficiency of this scheme is demonstrated on two numerical examples. The proposed approach produced better results compared with existing methods discussed in literature.

Inference of neutrino flavor evolution through data assimilation and neural differential equations

The evolution of neutrino flavor in dense environments such as core-collapse supernovae and binary compact object mergers constitutes an important and unsolved problem. Its solution has potential implications for the dynamics and heavy-element nucleosynthesis in these environments. In this paper, we build upon recent work to explore inference-based techniques for estimation of model parameters and neutrino flavor evolution histories. We combine data assimilation, ordinary differential equation solvers, and neural networks to craft an inference approach tailored for non-linear dynamical systems. Using this architecture, and a simple two-neutrino, two-flavor model, we test various optimization algorithms with the help of four experimental setups. We find that employing this new architecture, together with evolutionary optimization algorithms, accurately captures flavor histories in the four experiments. This work provides more options for extending inference techniques to large numbers of neutrinos.

Analytical evaluation of transformation matrix between functional spaces and its application in obtaining exact solutions of thermo-mechanical stresses in disks

Thermo-mechanical stresses developed in a disk, inner boundary of which is subject to the nonstationary thermo-mechanical loadings due to variation in temperature and pressure inside the disk is considered in the presented paper. Problems of theoretical analysis are divided and described as follows: The first problem implies analytical solution of the heat conduction equation with nonhomogeneous boundary conditions of Robin type, and the second problem deals with obtaining analytical solution of dynamic equation of displacement for nonisothermal elastic body (Navier equation) with traction boundary conditions in the plane strain case of uncoupled thermo-elasticity. Finite Hankel transform is used to solve the problems. Two sets of orthogonal basis functions necessary for integral transform are obtained. The first one corresponds to the boundary conditions written for the heat conduction equation, and the second fits the boundary conditions for Navier equation. Integrals defining elements of transformation matrix between the two sets of basis functions are evaluated analytically. In general, convenient method of obtaining analytical solution of linear hyperbolic equation with nonhomogeneous term, which itself is solution of nonhomogeneous parabolic equation, is elaborated within the presented work. Using finite difference method, discrete counterparts of the described problems are constructed and, via the built-in ordinary differential equation (ODE) solvers in MATLAB, numerical solutions are obtained, which exhibit excellent accordance with analytical ones.

A Novel Robust Strategy for Discontinuous Galerkin Methods in Computational Fluid Mechanics: Why? When? What? Where?

In this paper we will review a recent emerging paradigm shift in the construction and analysis of high order Discontinuous Galerkin methods applied to approximate solutions of hyperbolic or mixed hyperbolic-parabolic partial differential equations (PDEs) in computational physics. There is a long history using DG methods to approximate the solution of partial differential equations in computational physics with successful applications in linear wave propagation, like those governed by Maxwell’s equations, incompressible and compressible fluid and plasma dynamics governed by the Navier-Stokes and the Magnetohydrodynamics equations, or as a solver for ordinary differential equations (ODEs), e.g., in structural mechanics. The DG method amalgamates ideas from several existing methods such as the Finite Element Galerkin method (FEM) and the Finite Volume method (FVM) and is specifically applied to problems with advection dominated properties, such as fast moving fluids or wave propagation. In the numerics community, DG methods are infamous for being computationally complex and, due to their high order nature, as having issues with robustness, i.e., these methods are sometimes prone to crashing easily. In this article we will focus on efficient nodal versions of the DG scheme and present recent ideas to restore its robustness, its connections to and influence by other sectors of the numerical community, such as the finite difference community, and further discuss this young, but rapidly developing research topic by highlighting the main contributions and a closing discussion about possible next lines of research.

Benchmarking of numerical integration methods for ODE models of biological systems

Ordinary differential equation (ODE) models are a key tool to understand complex mechanisms in systems biology. These models are studied using various approaches, including stability and bifurcation analysis, but most frequently by numerical simulations. The number of required simulations is often large, e.g., when unknown parameters need to be inferred. This renders efficient and reliable numerical integration methods essential. However, these methods depend on various hyperparameters, which strongly impact the ODE solution. Despite this, and although hundreds of published ODE models are freely available in public databases, a thorough study that quantifies the impact of hyperparameters on the ODE solver in terms of accuracy and computation time is still missing. In this manuscript, we investigate which choices of algorithms and hyperparameters are generally favorable when dealing with ODE models arising from biological processes. To ensure a representative evaluation, we considered 142 published models. Our study provides evidence that most ODEs in computational biology are stiff, and we give guidelines for the choice of algorithms and hyperparameters. We anticipate that our results will help researchers in systems biology to choose appropriate numerical methods when dealing with ODE models.

Topology optimization of structural frames considering material nonlinearity and time-varying excitation

An approach for the topology optimization of structures composed of nonlinear beam elements under time-varying excitation is presented. Central to this approach is a hysteretic beam finite element model that accounts for distributed plasticity and axial-moment interaction through appropriate hysteretic interpolation functions and yield/capacity function, respectively. Nonlinearity is represented via the hysteretic variables for curvature and axial deformations that evolve according to first order nonlinear ordinary differential equations (ODEs), referred to as evolution equations, and the yield function. Hence, the governing dynamic equilibrium equations and hysteretic evolution equations can thus be concisely presented as a system of first-order nonlinear ODEs that can be solved using general ODE solvers without the need for linearization. The approach is applied for the design of frame structures with an objective to minimize the total volume in the domain, such that the maximum displacement at specified node(s) satisfies a specified constraint (i.e., drift limit) for the given excitation. The maximum displacement is approximated using the p-norm and thus permits the completion of the analytical sensitivities required for gradient-based updating. Several numerical examples are presented to demonstrate the approach for the design of structural frames subjected to pulse, harmonic, and seismic base excitation. Topologies obtained using the suggested, nonlinear approach are compared to solutions obtained from topology optimization problems assuming linear-elastic material behavior. These comparisons show that although similarities between the designs exist, in general the nonlinear designs differ in composition and, importantly, outperform the linear designs when assessed by nonlinear dynamic analysis.

High Order Block Method for Third Order ODEs

Many initial value problems are difficult to be solved using ordinary, explicit step-by-step methods because most of these problems are considered stiff. Certain implicit methods, however, are capable of solving stiff ordinary differential equations (ODEs) usually found in most applied problems. This study aims to develop a new numerical method, namely the high order variable step variable order block backward differentiation formula (VSVO-HOBBDF) for the main purpose of approximating the solutions of third order ODEs. The computational work of the VSVO-HOBBDF method was carried out using the strategy of varying the step size and order in a single code. The order of the proposed method was then discussed in detail. The advancement of this strategy is intended to enhance the efficiency of the proposed method to approximate solutions effectively. In order to confirm the efficiency of the VSVO-HOBBDF method over the two ODE solvers in MATLAB, particularly ode15s and ode23s, a numerical experiment was conducted on a set of stiff problems. The numerical results prove that for this particular set of problem, the use of the proposed method is more efficient than the comparable methods. VSVO-HOBBDF method is thus recommended as a reliable alternative solver for the third order ODEs.

Theoretical Analysis of Bulletproof Capability of Multilayer Ceramic Composites Subjected to Impact by an Armor Piercing Projectile

Inspired by the theories of Tate and Zaera, a theoretical analysis model including the erosion of the projectile, the cracking of ceramic composites, and the deformation of metal backplate was established in this study to investigate the bulletproof capability of the ceramic composites under impact by an armor piecing projectile (AP). The analysis results were verified by ballistic tests. As for the ceramic composites, the volume of the cracked ceramic conoid and the change in the compressive strength were included. Regarding the deformation of the metal backplate, the plastic deformation work, the external work, and the conservation of kinetics were considered. Based on the thickness of the target plate, failure modes were separated into the plug type and the petal type. The ordinary differential equation solver of MATLAB, ode45, was adopted to solve relevant ordinary differential equations. In this study, the powder metallurgy was used to produce the Al2O3/ZrO2 multilayered ceramic composites of three layers; each layer was 3 mm in thickness. The ceramic composites were paired with a backplate made of 6061‐T6 aluminum alloy with a thickness of either 1 mm or 4 mm. The ballistic tests were executed by using 0.30″ AP projectiles to impact the specimens. The results from theoretical model and ballistic tests were compared and shown consistent in the field of residual velocity, residual bullet mass, and the failure modes of the metal backplate.

A Novel Method for Solving Ordinary Differential Equations with Artificial Neural Networks

This research work investigates the use of Artificial Neural Network (ANN) based on models for solving first and second order linear constant coefficient ordinary differential equations with initial conditions. In particular, we employ a feed-forward Multilayer Perceptron Neural Network (MLPNN), but bypass the standard back-propagation algorithm for updating the intrinsic weights. A trial solution of the differential equation is written as a sum of two parts. The first part satisfies the initial or boundary conditions and contains no adjustable parameters. The second part involves a feed-forward neural network to be trained to satisfy the differential equation. Numerous works have appeared in recent times regarding the solution of differential equations using ANN, however majority of these employed a single hidden layer perceptron model, incorporating a back-propagation algorithm for weight updation. For the homogeneous case, we assume a solution in exponential form and compute a polynomial approximation using statistical regression. From here we pick the unknown coefficients as the weights from input layer to hidden layer of the associated neural network trial solution. To get the weights from hidden layer to the output layer, we form algebraic equations incorporating the default sign of the differential equations. We then apply the Gaussian Radial Basis function (GRBF) approximation model to achieve our objective. The weights obtained in this manner need not be adjusted. We proceed to develop a Neural Network algorithm using MathCAD software, which enables us to slightly adjust the intrinsic biases. We compare the convergence and the accuracy of our results with analytic solutions, as well as well-known numerical methods and obtain satisfactory results for our example ODE problems.