Analytic solution of nonlinear singularly perturbed initial value problems through iteration

This paper is concerned with singularly perturbed initial value problems for systems of ordinary differential equations. Here our emphasis will be on nonlinear phenomena and properties, particularly those with physical relevance. Since very few nonlinear systems can be solved explicitly, one must typically rely on a numerical scheme to accurately approximate the solution. However, numerical schemes do not always give accurate results, and we discuss the class of stiff differential equations, which present a more serious challenge to numerical analysts. In this paper, we derive in closed from, analytic solution of stiff nonlinear initial value problems, through iteration. The obtained sequence of iterates is based on the use of Lagrange multipliers. Moreover, the illustrative examples shows the efficiency of the method.     

Dynamic reduction of a CH4/air chemical mechanism appropriate for investigating vortex–flame interactions

This paper describes two methods, piecewise reusable implementation of solution mapping (PRISM) and dynamic steady‐state approximation (DYSSA), in which chemistry is reduced dynamically to reduce the computational burden in combustion simulations. Each method utilizes the large range in species timescales to reduce the dimensionality to the number of species with slow timescales. The methods are applied within a framework that uses hypercubes to partition multidimensional chemical composition space, where each chemical species concentration, plus temperature, is represented by an axis in space. The dimensionality of the problem is reduced uniquely in each hypercube, but the dimensionality of chemical composition space is not reduced. The dimensionality reduction is dynamic and is different for different hypercubes, thereby escaping the restrictions of global methods in which reductions must be valid for all chemical mixtures. PRISM constructs polynomial equations in each hypercube, replacing the chemical kinetic ordinary differential equation (ODE) system with a set of quadratic polynomials with terms related to the number of species with slow timescales. Earlier versions of PRISM were applied to smaller chemical mechanisms and used all chemical species concentrations as terms. DYSSA is a dynamic treatment of the steady‐state approximation and uses the fast–slow timescale separation to determine the set of steady‐state species in each hypercube. A reduced number of chemical kinetic ODEs are integrated rather than the original full set. PRISM and DYSSA are evaluated for simulations of a pair of counterrotating vortices interacting with a premixed CH₄/air laminar flame. DYSSA is sufficiently accurate for use in combustion simulations, and when relative errors are less than 1.0%, speedups on the order of 3 are observed. PRISM does not perform as well as DYSSA with respect to accuracy and efficiency. Although the polynomial evaluation that replaces the ODE solver is sufficiently fast, polynomials are not reused sufficiently to enable their construction cost to be recovered. © 2007 Wiley Periodicals, Inc. 39: 204–220, 2007     

Exact solutions to magnetogasdynamic equations in Lagrangian coordinates

The Lie group of point transformations, which leave the equations for a simplified model of one dimensional ideal gas in magnetogasdynamics invariant, are used to obtain some exact solutions for the governing system of hyperbolic partial differential equations (PDEs). Similarity variables which reduces the governing system of PDEs into system of ordinary differential equations (ODEs) are determined through the transformations. The resulting ODEs are solved analytically to obtain some exact solutions that exhibits space-time dependence. Further, we study the propagation of weak discontinuity through a state characterized by one of the solutions.     

Integration of chemical stiff ODEs using exponential propagation method

In this paper, we study the numerical long time integration of large stiff systems of differential equations arising from chemical reactions by exponential propagation methods. These methods, which typically converge faster, use matrix-vector products with the exponential or other related function of the Jacobian that can be effectively approximated by Krylov sabspace methods. We equip these methods to an automatic stepsize control technique and apply the method of order 4 for numerical integration of some famous stiff chemical problems such as Belousov-Zhabotinskii reaction, the Chapman atmosphere, Hydrogen chemistry, chemical Akzo-Nobel problem and air pollution problem.     

Numerical treatment of Casson nanofluid Bioconvectional flow with heat transfer due to stretching cylinder/plate: Variable physical properties

PurposeThe purpose of the current framework is to scrutinize the two-dimensional flow and heat transfer of Casson nanofluid over cylinder/plate along with impacts of thermophoresis and Brownian motion effects. Also, the effects of exponential thermal sink/source, bioconvection, and motile microorganisms are taken.Methodology/ApproachThe resulting non-linear equations (PDEs) are reformed into nonlinear ODEs by using appropriate similarity variables. The resultant non-linear (ODEs) were numerically evaluated by the use of the Bvp4c package in the mathematical solver MATLAB.FindingsThe numerical and graphical illustration regarding outcomes represents the performance of flow-involved physical parameters on velocity, temperature, concentration, and microorganism profiles. Additionally, the skin friction coefficient, local Nusselt number, local Sherwood number, and local microorganism density number are computed numerically for the current presented system. We noted that the velocity profile diminishes for the rising estimations of magnetic and mixed convection parameters. The Prandtl number corresponds with the declining performance of the temperature profile observed. The enhancement in the values of the Solutal Biot number and Brownian motion parameter increased in the concentration profile.OriginalityIn specific, this framework focuses on the rising heat transfer of Casson nanofluid with bioconvection by using a shooting mathematical model. The novel approach of the presented study is the use of motile microorganisms with exponential thermal sink/source in a Casson nano-fluid through a cylinder/plate. A presented study performed first time in the author’s opinion. Understanding the flow characteristics and behaviors of these nanofluids is crucial for the scientific community in the developing subject of nanofluids.     

On identifiability for chemical systems from measurable variables

The dynamics of the composition of chemical species in reacting systems can be characterized by a set of autonomous differential equations derived from mass conservation principles and some elementary hypothesis related to chemical reactivity. These sets of ordinary differential equations (ODEs) are basically non-linear, their complexity grows as much increases the number of substances present in the reacting media and can be characterized by a set of phenomenological constants (kinetic rate constants) which contains all the relevant information about the physical system. The determination of these kinetic constants is critical for the design or control of chemical systems from a technological point of view but the non-linear nature of the ODEs implies that there are hidden correlations between the parameters which maybe can be revealed with a identifiability analysis.     

A combined explicit-implicit method for high accuracy reaction path integration

We present the use of an optimal combined explicit-implicit method for following the reaction path to high accuracy. This is in contrast to most purely implicit reaction path integration algorithms, which are only efficient on stiff ordinary differential equations. The defining equation for the reaction path is considered to be stiff, however, we show here that the reaction path is not uniformly stiff and instead is only stiff near stationary points. The optimal algorithm developed in this work is a combination of explicit and implicit methods with a simple criterion to switch between the two. Using three different chemical reactions, we combine and compare three different integration methods: the implicit trapezoidal method, an explicit stabilized third order algorithm implemented in the code DUMKA3 and the traditional explicit fourth order Runge-Kutta method written in the code RKSUITE. The results for high accuracy show that when the implicit trapezoidal method is combined with either explicit method the number of energy and gradient calculations can potentially be reduced by almost a half compared with integrating either method alone. Finally, to explain the improvements of the combined method we expand on the concepts of stability and stiffness and relate them to the efficiency of integration methods.     

Locally implicit solution of a reaction-diffusion system with stiff kinetics

The Tyson-Fife reaction-diffusion equations are solved numerically using a locally implicit approach. Since the variables evolve at very different time scales, the resulting system of equations is stiff. The reaction term is responsible for the stiffness and the time step is increased by using an implicit method. The diffusion operator is evaluated explicitly and the system of implicit nonlinear equations is decoupled. The method is particularly useful for parameter values in which the equations are very stiff, such as the values obtained directly from the experimental reaction rate constants. Previous efforts modified the parameters on the equations to avoid stiffness. The equations then become a simplified model of excitable media and, for those cases, the locally implicit method gives a faster although less accurate solution. Nevertheless, since the modified equations no longer represent a particular chemical system an accurate solution is not as important. The algorithm is applied to observe the transition from simple motion to compound motion of a spiral tip.     

A new technique for the early detection of stiffness in coupled differential equations and application to standard Runge-Kutta algorithms

Higher-order Runge-Kutta (RK) algorithms employing local truncation error (LTE) estimates have had very limited success in solving stiff differential equations. These LTEs do not recognize stiffness until the region of instability has been crossed after which no correction is possible. A new technique has been designed, using the local stiffness function (LSF), which can detect stiffness very early before instability occurs. The LSF is a normalized dimensionless ratio which is essentially based on the product of the step size and the geometric mean of all the slopes. It is exceedingly sensitive to the onset of stiffness. Together, the LSF and the LTE form a complementary pair which can cooperate to help solve some mildly stiff equations which were previously intractable to RK algorithms alone. Examples are given of implementation and LSF performance. Received: 18 September 1997 / Accepted: 13 February 1998 / Published online: 17 June 1998     

Time scales and other problems in linking simulations of simple chemical systems to more complex ones

We shall start with very small systems like H₂ and H₃, computed with very accurate methods (Hylleraas–CI ) or atomic systems up to Zn with accurate methods (CI ), then move to more complex ones, like C₆₀, but now with somewhat less accurate methods, specifically Hartree–Fock with density functionals, the latter for the correlation energy but not for the exchange energy. For even more complex tasks like geometry optimization of C₆₀, we have resorted to even simpler and parametrized methods, like local density functionals. Then, we could use quantum mechanics either to provide interaction potentials for classical molecular dynamics or to directly solve dynamical systems, in a quantum molecular dynamics approximation. Having demonstrated that we can use the computational output from small systems as input to larger ones, we discuss in detail a new model for liquid water, which is borne out entirely from ab initio methods and nicely links spectroscopic, thermodynamics, and other physicochemical data. Concerning time scales, we use classical molecular dynamics to determine friction coefficients, and with these we perform stochastic dynamic simulations. The use of simulation results from smaller systems to provide inputs for larger system simulations is the “global simulation” approach, which, today, with the easily available computers, is becoming more and more feasible. Projections on simulations in the 1996–1998 period are discussed, new computational areas are outlined, and a N⁴ complexity algorithm is compared to density functional approaches. © 1993 John Wiley & Sons, Inc.     

Using ordinary differential equations system to solve isoconversional problems in non-isothermal kinetic analysis

The mathematical evaluation of the activation energy, E, of non-isothermal degradation reactions is usually made using the Ozawa/Flynn–Wall isoconversion principle and involves the numerical resolution of a set of integrals without closed form solution, which are solved by polynomial approximation or by numeric integration. In the present work, the isoconversion principle, originally described and maintained until now as an algebraic problem, was written as a set of ordinary differential equations (ODEs). The individual ODEs obtained are integrated by numeric methods and are used to estimate the activation energy of simulated examples. A least square error (LSE) objective function using the introduced ODEs was written to deal with multiple heating rate CaCO₃ thermal decomposition TG experiments.     

The ORAC Assay: Mathematical Analysis of the Rate Equations and Some Practical Considerations

ORAC (oxygen radical absorbance capacity), a method widely used for measuring the total antioxidant capacity of biological samples, can also be used for the determination of the relative reactivity of an antioxidant compound (XH) by examining the dependence of the rate of consumption of the probe (PH) on the concentration of XH; initial conditions are chosen in such a way that the rate of consumption of the starting reactants may be assumed to follow a drastically simplified kinetic scheme, and the steady‐state approximation for the concentration of the azo compound peroxyl (ROO^•) radical is invoked to simplify the analysis. Here we first attempted to find an analytical solution to the coupled first‐order ordinary differential equations (ODEs) of the minimal ORAC kinetic system, applying Lie symmetry group theory without any precondition. However, the Lie symmetry transformations applied to the Chini equation, which appeared after mathematical transformations, showed that the form of the coefficients of the Chini equation precluded the analytical solution of the minimal ORAC kinetic system through symmetry reduction. Consequently, an approximate analytical solution was sought, valid for the case when the bimolecular rate constant of XH with ROO^• (i.e., k_x ) was much larger than that of PH with ROO^• (i.e., k_p ). Using numerical solutions of the original set of ODEs of the ORAC kinetic system, the quality of the approximate solution was inspected under conditions that correspond to those employed in several ORAC methods together with a low initial concentration of the azo compound radical initiator. The simulations allowed us to conclude that the approximate analytical solution of the ODEs of the minimal ORAC kinetic system was not entirely devoid of academic interest, but its applicability was restricted to conditions where both k_x ≫ k_p and the initial concentration of XH was higher than that of PH.     

An approach for simultaneous estimation of reaction kinetics and curve resolution from process and spectral data

Parameter estimation of reaction kinetics from spectroscopic data remains an important and challenging problem. This study describes a unified framework to address this challenge. The presented framework is based on maximum likelihood principles, nonlinear optimization techniques, and the use of collocation methods to solve the differential equations involved. To solve the overall parameter estimation problem, we first develop an iterative optimization‐based procedure to estimate the variances of the noise in system variables (eg, concentrations) and spectral measurements. Once these variances are estimated, we then determine the concentration profiles and kinetic parameters simultaneously. From the properties of the nonlinear programming solver and solution sensitivity, we also obtain the covariance matrix and standard deviations for the estimated kinetic parameters. Our proposed approach is demonstrated on 7 case studies that include simulated data as well as actual experimental data. Moreover, our numerical results compare well with the multivariate curve resolution alternating least squares approach.     

A quick and efficient method for consistent initialization of battery models

Secondary batteries are usually modeled as a system of coupled nonlinear partial differential equations. These models are typically solved by applying finite differences or other discretization techniques in the spatial directions and solving the resulting system of differential algebraic equations (DAEs) numerically in time. These DAEs are very difficult to solve even using popular DAE solvers. The complications arise partly due to the difficulty in obtaining consistent, or closely consistent, initial conditions for the DAEs. In this paper, a shooting method is proposed as an effective and rapid technique for the initialization of battery models. This method is built on a regionwise shooting approach with initial guess at one end of the electrode and physics based shooting criterion on the other end that can ultimately satisfy all the required conditions in a battery unit. Notably, the computation time required for the proposed method is only milliseconds in a FORTRAN environment for the case of initializing a standard physics based lithium-ion battery model. Also the initial values obtained are exact and can readily be fed into any DAE solver for achieving accurate solutions without solver failure. This rapid method will help in simulating batteries in hybrid environments in real-time (milliseconds).     

Two reliable wavelet methods to Fitzhugh–Nagumo (FN) and fractional FN equations

Fractional reaction–diffusion equations serve as more relevant models for studying complex patterns in several fields of nonlinear sciences. In this paper, we have developed the wavelet methods to find the approximate solutions for the Fitzhugh–Nagumo (FN) and fractional FN equations. The proposed method techniques provide the solutions in rapid convergence series with computable terms. To the best of our knowledge, until now there is no rigorous wavelet solutions have been reported for the FN and fractional FN equations arising in gene propagation and model. With the help of Laplace operator and Legendre wavelets operational matrices, the FN equation is converted into an algebraic system. Finally, we have given some numerical examples to demonstrate the validity and applicability of the wavelet methods. The power of the manageable method is confirmed. Moreover, the use of the wavelet methods is found to be accurate, efficient, simple, low computation costs and computationally attractive.     

Stiffness and energy conservation in molecular dynamics: An improved integrator

Molecular dynamics is the integration of a set of coupled differential equations describing the motion of atoms over time. These equations exhibit the unfortunate property of stiffness, that is, terms of the equations (the forces on the atoms) are defined on several scales—ranging from tens of kcal/mol/Å to thousands of kcal/mol/Å. Additional nonconservative and stiff effects occur when a distance cutoff is used for the electrostatics and nonbonded potentials. Because the first derivative at the cutoff is essentially infinite, small variations in positions will cause large variations in energy and violate conservation of energy. The effects are demonstrated in a small system of 125 isolated water molecules. It is possible to greatly reduce and nearly eliminate the stiff integration effects with an improved integrator. The nonconservative effects of the distance cutoff cannot be removed by changing the integrator. © John Wiley & Sons, Inc.     

Wavelet formulation of the polarizable continuum model

The first implementation of a wavelet discretization of the Integral Equation Formalism (IEF) for the Polarizable Continuum Model (PCM) is presented here. The method is based on the application of a general purpose wavelet solver on the cavity boundary to solve the integral equations of the IEF‐PCM problem. Wavelet methods provide attractive properties for the solution of the electrostatic problem at the cavity boundary: the system matrix is highly sparse and iterative solution schemes can be applied efficiently; the accuracy of the solver can be increased systematically and arbitrarily; for a given system, discretization error accuracy is achieved at a computational expense that scales linearly with the number of unknowns. The scaling of the computational time with the number of atoms N is formally quadratic but a N^1.5 scaling has been observed in practice. The current bottleneck is the evaluation of the potential integrals at the cavity boundary which scales linearly with the system size. To reduce this overhead, interpolation of the potential integrals on the cavity surface has been successfully used. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Generalized CC-TDSCF and LCSA: The system-energy representation

Typical (sub)system-bath quantum dynamical problems are often investigated by means of (approximate) reduced equations of motion. Wavepacket approaches to the dynamics of the whole system have gained momentum in recent years and there is hope that properly designed approximations to the wavefunction will allow one to correctly describe the subsystem evolution. The continuous-configuration time-dependent self-consistent field (CC-TDSCF) and local coherent-state approximation (LCSA) methods, for instance, use a simple Hartree product of bath single-particle-functions for each discrete variable representation (DVR) state introduced in the Hilbert space of the subsystem. Here we focus on the above two methods and replace the DVR states with the eigenstates of the subsystem Hamiltonian, i.e., we adopt an energy-local representation for the subsystem. We find that stable and semiquantitative results are obtained for a number of dissipative problems, at the same (small) computational cost of the original methods. Furthermore, we find that both methods give very similar results, thus suggesting that coherent-states are well suited to describe (local) bath states. As a whole, present results highlight the importance of the system basis-set in the selected-multiconfiguration expansion of the wavefunction. They suggest that accurate and yet computationally cheap methods may be simply obtained from CC-TDSCF/LCSA by letting the subsystem states be variationally optimized.     

A family of numerical multistep methods with three distinct schemes: explicit advanced step-point (EAS) methods and the EAS1 approach

In this work, for the first time in an article, we present in a comprehensive way the explicit advanced step-point (EAS) methods. The EAS methods is a family of methods designed for the numerical solution of non-stiff and mildly stiff initial value problems (IVPs) and comprises three distinct schemes: EAS1, EAS2 and EAS3. A thorough theoretical analysis of the EAS family of predictor–corrector methods is presented in terms of their accuracy and stability characteristics and requirements, as well as the rationale for creating the three distinct schemes mentioned above. In this paper we also examine in detail one of the three schemes, the EAS1 methods. EAS1 are assessed for the very first time, are meticulously studied and their superior regions of absolute stability are presented. Furthermore the computational efficiency of EAS1 is examined and comparative numerical results are presented with the use of a variable step, variable order EAS1 code. The numerical results provide good evidence that EAS1 could be seen as superior to the well established Adams methods for the numerical solution of mildly stiff initial value problems.