Ordinary Differential Equation Solver Research Articles

Development of a fully implicit ODE-solver for containment analysis code

The thermal–hydraulic dynamics in containment are governed by a system of stiff ordinary differential equations (ODEs). A fully implicit discretization scheme is adopted to discretize these ODEs in order to mitigate the effects of stiffness. In comparison with explicit or semi-implicit discretization schemes that are subject to Courant limits on time steps, the fully implicit discretization scheme is more suitable for a containment analysis code that focuses on predicting both short-term and long-term thermal–hydraulic parameters after an accident. This study introduces a general-purpose ODE solver for the containment analysis code. The outline of the solver is as follows: The fully implicit discrete equations lead to a large set of nonlinear equations that need to be solved using Newton’s iterative method. The partial derivative components in the Jacobi matrix are calculated by the perturbation method using finite difference approximation, which avoids the complicated derivation of partial derivatives. The scaling modification technique is incorporated into this ODE solver to deal with significant differences in unknown variable magnitudes, and the line search method is introduced to address the difficulty of obtaining an accurate root estimate with Newton’s method when the initial guess is far from the actual root. This proposed ODE solver was applied to two typical stiff ODE problems to test its stiffness-suppressed ability and to demonstrate that this proposed solver can perform calculations with a very large time step. Then, the CASSIA code, a containment analysis code developed by China Nuclear Power Technology Research Institute Co., Ltd (CNPRI), equipped with this ODE solver, was applied to the CSNI (Committee on the Safety of Nuclear Installations) benchmark problem and the Carolinas Virginia Tube Reactor (CVTR) test 3 problem to preliminarily demonstrate that the proposed ODE solver can perform containment thermal–hydraulic analysis correctly. This study could provide references for the development of a home-made containment analysis code.

Determination of effective liquid-phase diffusion coefficients by tracer method in continuous packed columns

Many industrially relevant chemical processes are affected by the liquid-phase diffusion resistance in porous materials, for example in heterogeneous catalysts. Therefore, reliable experimental methods for the determination of effective diffusion coefficients in porous structures are highly desirable. In this work, a method was developed to study the internal diffusion effects in porous media by step response experiments with tracers injected to fixed beds filled with porous particles and inert material. The tracer responses were recorded at the outlet of the bed by UV–vis spectrometry. Ethanol was used as the tracer and water as the solvent. Porous aluminum oxide particles were selected as the model material. The step responses were described with a dynamic axial dispersion model, i.e. partial differential equations for the inert and catalyst sections of the bed. The partial differential equations describing the model were discretized with finite differences and the resulting ordinary differential equations were solved numerically with a solver for stiff ordinary differential equations. From the experimental data, the dispersion coefficients of the bed and the effective diffusion coefficient in the porous aluminum oxide particles were successfully determined by non-linear regression analysis. The values of the effective diffusion coefficents had a good agreement with the liquid-phase diffusion coefficients estimated from correlations. The method is applicable for any porous material.

LFT: Neural Ordinary Differential Equations With Learnable Final-Time.

Since the last decade, deep neural networks have shown remarkable capability in learning representations. The recently proposed neural ordinary differential equations (NODEs) can be viewed as the continuous-time equivalence of residual neural networks. It has been shown that NODEs have a tremendous advantage over the conventional counterparts in terms of spatial complexity for modeling continuous-time processes. However, existing NODEs methods entail their final time to be specified in advance, precluding the models from choosing a desirable final time and limiting their expressive capabilities. In this article, we propose learnable final-time (LFT) NODEs to overcome this limitation. LFT rebuilds the NODEs learning process as a final-time-free optimal control problem and employs the calculus of variations to derive the learning algorithm of NODEs. In contrast to existing NODEs methods, the new approach empowers the NODEs models to choose their suitable final time, thus being more flexible in adjusting the model depth for given tasks. Additionally, we analyze the gradient estimation errors caused by numerical ordinary differential equations (ODEs) solvers and employ checkpoint-based methods to obtain accurate gradients. We demonstrate the effectiveness of the proposed method with experimental results on continuous normalizing flows (CNFs) and feedforward models.

Physics-Informed Neural Networks for High-Frequency and Multi-Scale Problems Using Transfer Learning

Physics-Informed Neural Network (PINN) is a data-driven solver for partial and ordinary differential equations (ODEs/PDEs). It provides a unified framework to address both forward and inverse problems. However, the complexity of the objective function often leads to training failures. This issue is particularly prominent when solving high-frequency and multi-scale problems. We proposed using transfer learning to boost the robustness and convergence of training PINN, starting training from low-frequency problems and gradually approaching high-frequency problems through fine-tuning. Through two case studies, we discovered that transfer learning can effectively train PINNs to approximate solutions from low-frequency problems to high-frequency problems without increasing network parameters. Furthermore, it requires fewer data points and less training time. We compare the PINN results using direct differences and L2 relative error showing the advantage of using transfer learning techniques. We describe our training strategy in detail, including optimizer selection, and suggest guidelines for using transfer learning to train neural networks to solve more complex problems.

A hybrid kernel-based meshless method for numerical approximation of multidimensional Fisher’s equation

We propose and analyze a meshless method of lines by considering some hybrid radial kernels. These hybrid kernels are constructed by linearly combining infinite smooth radial functions to piecewise smooth radial functions; which are then used for spatial approximation on trial spaces spanned by translates of positive definite radial functions. After spatial approximation, a high-order ODE solver is invoked for efficient and stable time-integration of the resultant semi-discrete system of ordinary differential equations (ODEs). Unlike the mesh-based method of lines, the proposed method works for arbitrary scattered data points and is equally effective for problems over non-rectangular domains. The proposed method is tested on one-, two- and three-dimensional reaction–diffusion Fisher equation for its numerical stability, accuracy, and efficiency against the contemporary meshless and mesh-based methods. The economical computational cost, improved accuracy, eigenvalues stability, and well-conditioning of system matrices are observed against RBF collocation and RBF-PS methods.

Physics-Informed Neural Networks with skip connections for modeling and control of gas-lifted oil wells

Neural networks, while powerful, often lack interpretability. Physics-Informed Neural Networks (PINNs) address this limitation by incorporating physics laws into the loss function, making them applicable to solving Ordinary Differential Equations (ODEs) and Partial Differential Equations (PDEs). The recently introduced PINC framework extends PINNs to control applications, allowing for open-ended long-range prediction and control of dynamic systems. In this work, we enhance PINC for modeling highly nonlinear systems such as gas-lifted oil wells. By introducing skip connections in the PINC network and refining certain terms in the ODE, we achieve more accurate gradients during training, resulting in an effective modeling process for the oil well system. Our proposed improved PINC demonstrates superior performance, reducing the validation prediction error by an average of 67% in the oil well application and significantly enhancing gradient flow through the network layers, increasing its magnitude by four orders of magnitude compared to the original PINC. Furthermore, experiments showcase the efficacy of Model Predictive Control (MPC) in regulating the bottom-hole pressure of the oil well using the improved PINC model, even in the presence of noisy measurements.

Emulating the interstellar medium chemistry with neural operators

Context. The study of galaxy formation and evolution critically depends on our understanding of the complex photo-chemical processes that govern the evolution and thermodynamics of the interstellar medium (ISM). In a computational sense, resolving the chemistry is among the weightiest tasks in cosmological and astrophysical simulations. Aims. Astrophysical simulations can include photo-chemical models that allow for a wide range of densities (n), abundances of different species (ni/n) and temperature (T), and plausible evolution scenarios of the ISM under the action of a radiation field (F) with different spectral shapes and intensities. The evolution of such a non-equilibrium photo-chemical network relies on implicit, precise, computationally costly, ordinary differential equations (ODE) solvers. Here, we aim to substitute such procedural solvers with fast, pre-trained emulators based on neural operators. Methods. We emulated a non-equilibrium chemical network up to H2 formation (9 species, 52 reactions) by adopting the DeepONet formalism, namely: by splitting the ODE solver operator that maps the initial conditions and time evolution into a tensor product of two neural networks (named branch and trunk). We used KROME to generate a training set, spanning −2 < log(n/cm−3) ≤ 3.5, log(20) ≤ log(T/K) ≤ 5.5, −6 ≤ log(ni/n) < 0, and adopting an incident radiation field, F, sampled in 10 energy bins with a continuity prior. We separately trained the solver for T and each ni for ≃4.34 GPUhrs. Results. Compared with the reference solutions obtained by KROME for single-zone models, the typical precision obtained is of the order of 10−2, that is, it is 10 times better when using a training that is 40 times less costly, with respect to previous emulators that only considered a fixed F. DeepONet also performs well for T and ni outside the range of the training sample. Furthermore, the emulator aptly reproduces the ion and temperature profiles of photo dissociation regions as well; namely, by giving errors that are comparable to the typical difference between various photo-ionization codes. The present model achieves a speed-up of a factor of 128× with respect to stiff ODE solvers. Conclusions. Our neural emulator represents a significant leap forward in the modelling of ISM chemistry, offering a good balance of precision, versatility, and computational efficiency. Nevertheless, further work is required to address the challenges represented by the extrapolation beyond the training time domain and the removal of potential outliers.

Spatial-temporal graph neural ODE networks for skeleton-based action recognition.

In the field of skeleton-based action recognition, accurately recognizing human actions is crucial for applications such as virtual reality and motion analysis. However, this task faces challenges such intraindividual action differences and long-term temporal dependencies. To address these challenges, we propose an innovative model called spatial-temporal graph neural ordinary differential equations (STG-NODE). First, in the data preprocessing stage, the dynamic time warping (DTW) algorithm is used to normalize and calculate 3D skeleton data to facilitate the derivation of customized adjacency matrices for improving the influence of intraindividual action differences. Secondly, a custom ordinary differential equation (ODE) integrator is applied based on the initial conditions of the temporal features, producing a solution function that simulates the dynamic evolution trend of the events of interest. Finally, the outstanding ODE solver is used to numerically solve the time features based on the solution function to increase the influence of long-term dependencies on the recognition accuracy of the model and provide it with a more powerful temporal modeling ability. Through extensive experiments conducted on the NTU RGB+D 60 and Kinetics Skeleton 400 benchmark datasets, we demonstrate the superior performance of STG-NODE in the action recognition domain. The success of the STG-NODE model also provides new ideas and methods for the future development of the action recognition field.

Comparative Analysis of Adams-Bashforth-Moulton and Runge-Kutta Methods for Solving Ordinary Differential Equations Using MATLAB

Comparative Analysis of Adams-Bashforth-Moulton and Runge-Kutta Methods for Solving Ordinary Differential Equations Using MATLAB

On the Solution of Blasius Boundary Layer Equations of Prandtl-Eyring Fluid Flow Past a Stretching Sheet

Objective: This paper investigates velocity profile for two-dimensional, incompressible, laminar forced convection flow of the fluid model for Prandtl-Eyring fluid past a stretching sheet in the presence of fluid parameters. Methods: The governing partial differential equation for the flow was transformed into non-linear ordinary differential equation by using the deductive one parameter group theoretic method and numerical solution of non-linear ordinary differential equation (ODE) is solved by MATLAB bvp4c solver. Findings: The solution of velocity profile obtained as a function of parameter and . The effect of the fluid parameter was discussed graphically. Novelty: The main goal of this article is to analyze boundary layer flow of Prandtl-Eyring fluid over a stretching surface. The conservation equations of mass, momentum are converted into non-linear ordinary differential equations along with boundary conditions using deductive one parameter group theoretic method and solved by MATLAB ODE solver. Comparisons with previously published works are made, and results show a high level of agreement. This type of research is applicable to extrusion, paper production, fiber glass production, hot rolling, condensation process, crystal growing, polymer sheets etc. Keywords: Boundary layer, laminar flow, Deductive one parameter Group theoretic method, Absolute invariant, Stretching Sheet, Prandtl-Eyring fluid

Understanding the impact of numerical solvers on inference for differential equation models.

Most ordinary differential equation (ODE) models used to describe biological or physical systems must be solved approximately using numerical methods. Perniciously, even those solvers that seem sufficiently accurate for the forward problem, i.e. for obtaining an accurate simulation, might not be sufficiently accurate for the inverse problem, i.e. for inferring the model parameters from data. We show that for both fixed step and adaptive step ODE solvers, solving the forward problem with insufficient accuracy can distort likelihood surfaces, which might become jagged, causing inference algorithms to get stuck in local 'phantom' optima. We demonstrate that biases in inference arising from numerical approximation of ODEs are potentially most severe in systems involving low noise and rapid nonlinear dynamics. We reanalyse an ODE changepoint model previously fit to the COVID-19 outbreak in Germany and show the effect of the step size on simulation and inference results. We then fit a more complicated rainfallrun-off model to hydrological data and illustrate the importance of tuning solver tolerances to avoid distorted likelihood surfaces. Our results indicate that, when performing inference for ODE model parameters, adaptive step size solver tolerances must be set cautiously and likelihood surfaces should be inspected for characteristic signs of numerical issues.

Three-dimension holographic numerical simulation of gravity anomalies

In gravity-anomaly-based prospecting, it is difficult to achieve large-scale and high-precision inversion imaging of complex geological models. To address this issue, this paper proposes a three-dimensional holographic numerical simulation method for gravity anomalies. This method transforms the 3D partial differential equation of gravity potential into many independent 1D differential equations with different wavenumbers by performing a 2D Fourier transform along the horizontal direction. It decomposes a large-scale 3D numerical modeling problem into many 1D numerical modeling problems, which greatly reduces the computation and memory requirements of modeling, and each 1D differential equation is independent, so it has high parallelism. The vertical direction is reserved as the spatial domain, which has strict upper and lower boundary conditions; The horizontal 2D Fourier transform uses a holographic Fourier transform (Holo-FT) with a full interval of integration, consistent with the physical boundary in the horizontal direction, thus the physical information of the gravity potential described by the three-dimensional partial differential equation is fully and accurately simulated. Using the high parallelism of the algorithm, the CPU is used for parallel solving ordinary differential equations, while the GPU is used for parallel computing of Holo-FT, which implements the CPU-GPU parallel acceleration scheme and further improves the efficiency of this algorithm.

Solving Ordinary Differential Equations with Boundary Conditions Numerically

This article elucidates the method and implementation of finite differences, aiming to enhance clarity for students and provide detailed insights into its application, as many authors assume a pre-existing understanding of numerical analysis. It includes two case studies employing Dirichlet and Cauchy boundary conditions, demonstrating the precision and effectiveness of this numerical tool in obtaining results.

A Sparse Differential Algebraic Equation (DAE) and Stiff Ordinary Differential Equation (ODE) Solver in Maple

This paper implements efficient numerical methods in Maple to solve index-1 nonlinear Differential Algebraic Equations (DAEs) and stiff Ordinary Differential Equations (ODEs) systems. Single-step methods (like Trapezoid (TR), Implicit-mid point (IMP), Euler-backward (EB), Radau IIA (Rad) methods, TRBDF2, TRX2) and backward-difference formula of order 2 are implemented with adaptive time-stepping methods in Maple to solve index-1 nonlinear DAEs. Maple’s robust and efficient ability to search within a list/set is exploited to identify the sparsity pattern and automatically calculate the analytic Jacobian. The algorithm and implementation are robust and efficient for index-1 DAE problems and scale well for finite difference/finite element discretization of two-dimensional models with system size up to 10,000 nonlinear DAEs and solve the same in a few seconds. The computational efficiency of the proposed algorithm (provided as an open-access code) compares favorably with the commercial solver gPROMs, one of the most commonly used sparse DAE solvers in the industry.

Accurate error estimation for model reduction of nonlinear dynamical systems via data-enhanced error closure

Accurate error estimation is crucial in model order reduction, to obtain small reduced-order models as well as to certify their accuracy when deployed in downstream applications such as digital twins. In existing a posteriori error estimation approaches, knowledge about the time integration scheme is mandatory, e.g., the residual-based error estimators proposed for the reduced basis method. This poses a challenge when automatic ordinary differential equation solver libraries are used to perform the time integration. To address this, we present a data-enhanced approach for a posteriori error estimation. Our new formulation enables residual-based error estimators to be independent of any time integration method. To achieve this, we introduce a corrected reduced-order model which takes into account a data-driven closure term for improved accuracy. The closure term, subject to mild assumptions, is related to the local truncation error of the corresponding time integration scheme. We propose efficient computational schemes for approximating the closure term, at the cost of a modest amount of training data. Furthermore, the new error estimator is incorporated within a greedy process to obtain parametric reduced-order models. Numerical results on three different systems show the accuracy of the proposed error estimation approach and its ability to produce ROMs that generalize well.

On problematic practice of using normalization in self-modeling/multivariate curve resolution (S/MCR)

This short communication paper is briefly dealing with greater or lesser misused normalization in self-modeling/multivariate curve resolution (S/MCR) practice. The importance of the correct use of the ODE (ordinary differential equation) solvers and apt kinetic illustrations are elucidated. The new terms, external and internal normalizations are defined and interpreted. The problem of reducibility of a matrix is touched. Improper generalization/development of normalization-based methods are cited as examples. The position of the extreme values of the signal contribution function is clarified. An Executable Notebook with Matlab Live Editor was created for algorithmic explanations and depictions.

Hyperbolic modeling of gradient damage and one-dimensional finite volume simulations

This study provides a new formulation of gradient damage model which allows an efficient explicit numerical solution of dynamics problems. The proposed methodology is based on an ”extended Lagrangian approach” developed by one of the authors for the nondissipative and dispersive shallow water equation. By using this strategy, the global minimization problem commonly derived for gradient damage models is recast as a purely local hyperbolic one with source terms that can be easily solved using finite volumes. The numerical solution of the governing system is then based on a fractional-step method consisting of a classical Godunov-type scheme and an implicit Ordinary Differential Equation solver for the local source terms. Numerical results are presented on one-dimensional multi-fragmentation and spalling tests for illustration purpose.

A fast and efficient computation method for reflective diffraction simulations

We present a new computation method for simulating reflection high-energy electron diffraction and the total-reflection high-energy positron diffraction experiments. The two experiments are used commonly for the structural analysis of material surface. The present paper improves the conventional numerical method, the multi-slice method, for faster computation, since the present method avoids the matrix-eigenvalue solver for the computation of matrix exponentials and can adopt higher-order ordinary differential equation solvers. Moreover, we propose a high-performance implementation based on multi-thread parallelization and cache-reusable subroutines. In our tests, this new method performs up to 2,000 times faster than the conventional method without compromising accuracy.

Applications of particle swarm optimization for numerical simulation of Fisher’s equation using RBF

This research article focused on optimising the parameter related to radial basis function (RBF) by solving the Fisher’s equation numerically. To obtain the numerical solutions, Particle Swarm Optimisation (PSO) algorithm is used along with radial basis functions (RBFs). PSO has been used by a number of researchers to obtain the optimised results in various scientific applications. In this work, using RBF-pseudospectral method, partial differential equations (PDEs) are transformed into a system of ordinary differential equations (ODEs) then are solved in MATLAB using the ODE solver. Two numerical examples are presented to demonstrate the effectiveness and accuracy of the proposed method by computingL∞, L2, and Lrms error norms with their order of convergence. Also, it can be demonstrated from the tables that the obtained results are in conformity with the results available in the literature.

Mathematical solutions in internal dose assessment: A comparison of Python-based differential equation solvers in biokinetic modeling

In biokinetic modeling systems employed for radiation protection, biological retention and excretion have been modeled as a series of discretized compartments representing the organs and tissues of the human body. Fractional retention and excretion in these organ and tissue systems have been mathematically governed by a series of coupled first-order ordinary differential equations (ODEs). The coupled ODE systems comprising the biokinetic models are usually stiff due to the severe difference between rapid and slow transfers between compartments. In this study, the capabilities of solving a complex coupled system of ODEs for biokinetic modeling were evaluated by comparing different Python programming language solvers and solving methods with the motivation of establishing a framework that enables multi-level analysis. The stability of the solvers was analyzed to select the best performers for solving the biokinetic problems. A Python-based linear algebraic method was also explored to examine how the numerical methods deviated from an analytical or semi-analytical method. Results demonstrated that customized implicit methods resulted in an enhanced stable solution for the inhaled 60Co (Type M) and 131I (Type F) exposure scenarios for the inhalation pathway of the International Commission on Radiological Protection (ICRP) Publication 130 Human Respiratory Tract Model (HRTM). The customized implementation of the Python-based implicit solvers resulted in approximately consistent solutions with the Python-based matrix exponential method (expm). The differences generally observed between the implicit solvers and expm are attributable to numerical precision and the order of numerical approximation of the numerical solvers. This study provides the first analysis of a list of Python ODE solvers and methods by comparing their usage for solving biokinetic models using the ICRP Publication 130 HRTM and provides a framework for the selection of the most appropriate ODE solvers and methods in Python language to implement for modeling the distribution of internal radioactivity.