On a Variational Method for Stiff Differential Equations Arising from Chemistry Kinetics

: For the approximation of stiff systems of ODEs arising from chemistry kinetics, implicit integrators emerge as good candidates. This paper proposes a variational approach for this type of systems. In addition to introducing the technique, we present its most basic properties and test its numerical performance through some experiments. The main advantage with respect to other implicit methods is that our approach has a global convergence. The other approaches need to ensure convergence of the iterative scheme used to approximate the associated nonlinear equations that appear for the implicitness. Notice that these iterative methods, for these nonlinear equations, have bounded basins of attraction.


Introduction
We start with the following stiff differential problem: Due to its stiffness, it is convenient to use implicit Runge-Kutta methods for the time integration step.The error analysis and the existence and uniqueness of the solution obtained through the Runge-Kutta method were analyzed in [1,2].However, a difficulty appeared related to the approximation of the nonlinear Runge-Kutta equations by iterative methods [2], since we need to find good initial guesses.We will point out these difficulties in the Motivation Section.
The variational approach that we will use is based on the error functional: that will be minimized among the absolutely-continuous path x : (0, T) → R N with x(0) = x 0 .Note that if E(x) is finite for one such path x, then automatically, x is square integrable.The error functional in (1) is associated with the original problem in a natural way: it only has a local minimum that is the solution of the problem.This property will imply the global convergence of our approach based on optimality conditions and minimization schemes like (steepest) descent methods.We only need to approximate linear problems; therefore, we can use the well-developed theory of the convergence of Runge-Kutta methods for stiff linear problems.
We will focus our attention on kinetic chemical reactions modeled through stiff differential equations.We will provide mathematical expressions for the reaction rate for problems where this rate is widely varying [3,4].The main goal of this article is to test the variational approach exposed before to solve this kind of problems.

Motivation
One of the simplest stiff problems is the following linear test differential equation, where λ is a complex number with Re λ << 0. This equation is used to test the A-stability and imposes restrictions on the step size of explicit methods.In order to point out the problems with the initial guesses in the use of implicit methods, we will consider a nonlinear modification of the problem in (2).We consider the following nonlinear variation, The explicit solution is given by the expression, In order to solve this problem, we can consider the classical trapezoid rule, The approximate solution of the problem in (2) obtained using this classical implementation is shown in Figure 1.To obtain these results, we have used λ = −100 (left) and λ = −10,000 (right).We can see that the method provides a bad approximation when the stiffness of the problem increases, as shown in Figure 1 on the right.This is due to the fact that the initial guess is outside of the basin of attraction of Newton's method [1].However, looking at Figure 2, we can see that the variational method converges in the two considered cases.

Some Theoretical Analysis
One initial assumption that we must make about f : R N → R N is its smoothness, so that ∇ f : R N → R N×N is continuous.We will also assume global Lipschitzianity with Lipschitz constant M > 0 (|∇ f | ≤ M).Moreover, we will suppose that for every positive C > 0 and small > 0, there is a constant D C, > 0 so that, This regularity comes from the relation of our approach with optimality.On the other hand, these properties are verified for most of the important problems in chemistry kinetics.
The next theorem is a local existence result.Fortunately, if the function f verifies the hypotheses globally, applying this theorem to successively smaller intervals, we can derive a global strong convergence result.
Theorem 1.Our iterative procedure x (j) = x (j−1) + y (j) , starting from any initial approximation x (0) , where: converges strongly in L ∞ (0, T) and in H 1 (0, T) to the unique solution of the original problem.
All the mathematical models considered in this work will verify the theoretical hypotheses since they will have both gradients and Hessians bounded.

Sketch of the proof:
We start with any initial approximation x (0) (≡ x), for instance x (0) = x 0 for all t, or x (0) (t) = x 0 + t f (x 0 ).In order to reduce the error finding a descent direction, we compute the Gâteaux derivative of E in (1) at a given x in the direction y with y(0) = 0. Namely, In particular, we can propose to choose y such that, This way it is clear that E (x)y = −2E(x), and the (local) decrease of the error is of the size E(x).Finding y requires solving the above linear problem.In some sense, this is like a Newton method.
Through an induction process, we obtain the proof of the theorem.The complete proof of Theorem 1 can be found in [5].

Some Numerical Examples
In this section, we introduce some well-known chemical problems and apply our iterative scheme in order to obtain an approximation of the solution and to show the efficiency of the new method.The models are: the Robertson problem, the chemical Akzo Nobel problem, and the Hires problem.For the three models, our approach is convergent.Indeed, we have: Corollary 1.The three models (Robertson, chemical Akzo Nobel, and Hires) verify the hypotheses of Theorem 1.
Proof.The three models considered in this section verify the theoretical hypotheses of the variational approach.Two of them are quadratic polynomials that have a constant Hessian, and the third is based on concentrations (all positive and bounded), in monomials and square roots.Therefore, all are sufficiently differentiable with bounded gradients and Hessians.
In particular, the user can consider the variational method as a black box with theoretical guarantees.This is not the case of other approaches, like implicit Runge-Kutta, where only numerical evidence is presented.For a more general context, of the chemistry kinetics and its applications, we refer to [6].
The choice of initial conditions when integrating systems of ODEs is usually not trivial, especially when the equations are stiff, and hence, the solution is not easily predictable.Therefore, the introduction of a new method that, under reasonable hypotheses, allows for a global convergence is rather important in our opinion.The figures correspond to the approximation of our method using a stopping criterion when the error functional is smaller than 10 −5 , and the global convergence of our approach ensures the good approximation to the solution.In fact, our methods can be used to check the existence of a solution, since the method only diverges basically when the problem has no solution.In any case, the values and the plots are those expected for these problems.

Robertson Problem
The following chemical reaction process was introduced by Robertson in 1966 [7]: The problem consists of three equations where k 1 , k 2 , and k 3 represent the rate constants and A, B, and C are the chemical species involved.
Under some conditions and applying the mass action law for the rate functions, the following mathematical model, which consists of a stiff system of three non-linear ordinary differential equations, is obtained, ( In the previous system, y 1 (t), y 2 (t), and y 3 (t) were the concentrations of the chemical species A, B, and C, respectively.The initial values at time t = 0 can be represented by (y 01 , y 02 , y 03 ) T .The values of the constants used in the test were k 1 = 0.04, k 2 = 3.10 7 , and k 3 = 10 4 , and the initial concentrations were y 01 = 1, y 02 = 0, and y 03 = 0.The time interval considered was 0 ≤ t ≤ 40.The stiffness of this problem is due to the large difference between the kinetic constants k i , which produce a very quick initial transient, as can be seen in Figure 3.The convergence of the method has been theoretically proven in [5].In Figure 3, it is shown that the method converges and approximates well the solution of the Robertson problem.The numerical value of the solution for the three components y 1 , y 2 , y 3 of the Robertson problem at t = 40, using the variational method, can be seen in Table 1.

Chemical Akzo Nobel Problem
This problem arose in Akzo Nobel Central Research (The Netherlands) [8].The problem describes a chemical process, in which two species, FLBand ZLU, are mixed, while carbon dioxide is continuously added.The resulting species of importance is ZLA.The names of the different chemical species are fictitious due to the commercial competition in this field.The chemical reaction scheme is the following, FLB + ZHU FLBT.ZHU K s represents the constant of equilibrium of the last equation and has the expression, Brackets represent concentrations.The velocities of the previous chemical equations are given respectively by: F in defines the inflow of oxygen per unit volume and is given by: where kl A denotes the mass transfer coefficient, H is the Henry constant, and p (O 2 ) represents the partial pressure of oxygen and is considered independent of [O 2 ]. k 1 , k 2 , k 3 , k 4 , K, kl A, H, and p (O 2 ) are constant parameters.The mathematical model consists of a stiff system of six non-linear ordinary differential equations.
where the auxiliary variables r i and F in are: The values of the parameters were k 1 = 18.7, k 4 = 0.42, K s = 115.83k 2 = 0.58, K = 34.4,p (O 2 ) k 3 = 0.009, kl A = 3.3, H = 737.Initially, 0.44 mol/liter of FLB were mixed with 0.007 mol/liter of ZHU.At the beginning, the carbon dioxide had a concentration of 0.00123 mol/liter.The time of simulation was 180 minutes.The relation between chemical and mathematical formulation was that the concentration of the different components in chemical formulation [ZLA], [FLB.ZHU] was respectively y 1 , . . ., y 6 in the mathematical formulation.The method was able to approximate the solution of this stiff problem.The convergence of the method can be seen in Figure 4. Table 2 shows the numerical solution of the components of the chemical Akzo Nobel problem y 1 , y 2 , y 3 , y 4 , y 5 , y 6 , at t = 180, when the variational method is used.

Hires Problem
Schäfer introduced this problem in 1975 [9].It represents the high irradiance response (HIRES) of photomorphogenesis on the basis of the phytochrome.The chemical reactions have eight reactants.The problem can be formulated mathematically through a stiff system of eight non-linear differential equations.The initial values were given by the vector y = (1, 0, 0, 0, 0, 0, 0, 0.0057) T .The end points of the time interval for each component have been arbitrarily chosen, as can be seen in the numerical results (Figure 5).For the components y 1 , y 2 , y 3 , y 4 , y 7 , and y 8 , we have chosen the interval [0, 5], and for the components y 5 and y 6 , we have chosen the interval [0, 350].The problem has eight components, and it may be considered as a large number of components.In this case, the method also converges to the solution sought; see Figure 5.The numerical value of the solution for the eight components y 1 , y 2 , y 3 , y 4 , y 7 y 8 at t = 5 and y 5 , y 6 at t = 350 of the Hires problem, using the variational method, can be seen in Table 3.

Conclusions
In stiff systems of ODEs, the restriction of the time step becomes severe, particularly when the simulation time is large.In this paper, the new variational method proposed has been successfully applied to several initial value problems arising from chemical reactions composed of large systems of stiff ODEs.Moreover, the variational method has the advantage of never getting stuck at local minima.It always converges to the solution regardless of the initialization.This is an important improvement over the methods that need an initialization that is close to the solution of the problem like the implicit Runge-Kutta methods.

Figure 3 .
Figure 3. Numerical solution of the concentrations y 1 (top left), y 2 (top right), and y 3 (bottom left) and a zoom of the transient around the origin of y 2 (bottom right) of the Robertson problem (3).

Figure 5 .
Figure 5. From left to right and top to bottom, we present the numerical solution of the HIRES problem for t ∈ [0, 5] min for the concentrations y 1 , y 2 , y 3 , y 4 , y 7 , y 8 and for t ∈ [0, 350] min for y 5 and y 6 .

Table 1 .
Results of the Robertson problem at t = 40.

Table 3 .
Results of the Hires problem at t = 5 for y 1 , y 2 , y 3 , y 4 , y 7 , y 8 and at t = 350 for y 5 and y 6 .