On application of solution continuation method with respect to the best exponential argument in solving stiff boundary value problems

Full Text

1. Introduction Real-world phenomena are increasingly being modelled using mathematical representations that encompass significant levels of complexity. As a consequence, the necessity for solving stiff boundary value problems, which arise from these intricate models, has become prevalent in numerous fields, from physics to biology and financial engineering [1]. Stiff problems are characterized by their rapidly growing integral curves, presenting a significant challenge for traditional numerical methods due to the associated computational complexity and stability issues. The shooting method, a popular approach for solving boundary value problems, is well-known for its shortcomings when dealing with stiff problems [2]. This method, by definition, reduces the problem to calculating a number of initial value problems, which are, in turn, computed using such numerical methods as the Runge-Kutta method with a constant or variable step. Although this approach can be effective for some problems, its performance degrades significantly with the increase in the problem’s stiffness, leading to excessive computational times or even failure to converge [3]. Against this backdrop, the development of new approaches to deal with stiff boundary value problems has become a pressing research topic. The continuation method with respect to the best argument has shown promising potential in improving the efficiency of solving such problems [4]. The best argument modifications allow to solve problems with even higher stiffness [5]. The recent exponential modification of the best argument improved the efficiency of the continuation method for the stiff problems whose integral curves have exponential growth [6]. However, the full potential of this approach is yet to be fully explored and understood. This research is aimed at addressing this gap in the literature. 2. Methodology In this section, we provide a detailed exploration of the methodology underpinning the continuation methods for solving stiff initial value problems. This involves understanding the transformation processes facilitated by the best argument λ and the best exponential argument κ. Consider the ordinary differential system: dyi = f (t, y , y , . . . , y ) , i = 1, . . . , n (1) dt i 1 2 n with initial conditions: yi(0) = yi,0, i = 1, . . . , n. (2) 1. Method of continuation with respect to the best argument To address the challenges posed by stiff and ill-conditioned Cauchy problems, implicit methods offer solutions. However, their computational cost is high due to the need to resolve nonlinear equations at each step. The continuation method provides a solution for this issue, introducing a new argument to E. D. Tsapko et al., On application of solution continuation method … 377 the Cauchy problem [4]. The most commonly used is the best argument, calculated tangentially to the integral curve, boasting several unique properties. For the baseline problem (1)-(2), the best argument λ is expressed as: dλ2 = dy2 + · · · + dy2 + dt2. (3) 1 n All variables and the time argument are functions of λ. Enhancing system (1) with relationship (3), and resolving for λ derivatives, we obtain: dyi dλ dt j , = fi (t, y1, . . . , yn) Q (t, y1, . . . , yn) 1 (4) dλ = jQ (t, y , . . . , y , i = 1, . . . , n, ) 1 n with the associated initial conditions: yi(0) = yi,0, t(0) = t0, i = 1, . . . , n, (5) where Q (t, y1, . . . , yn) = 1 + f 2 (t, y1, . . . , yn) + · · · + f 2 (t, y1, . . . , yn). 1 n This transformation imbues the numerical problem with several beneficial properties. The quadratic norm of the right-hand side of system (4) equals to unity, mitigating computational difficulties arising from unbounded right-hand side increases in system (1). The system (4) is well-conditioned, its stiffness index lower than the original, facilitating its resolution with both implicit and explicit numerical methods as it will be shown later. 2. Method of continuation with respect to the best exponential argument For systems where integral curves exhibit exponential growth [6], we introduce the best exponential argument κ: dκ2 = dy2 + · · · + dy2 + e-2αtdt2. (6) 1 n Transforming system (1) to the argument κ results in: dyi dκ dt = fi(t, y1, . . . , yn) · exp(αt) jQt(t, y1, . . . , yn) exp(αt) , i = 1, . . . , n, (7) dκ = jQt(t, y , . . . , y ) , 1 n with initial conditions same as (5) and Qt (t, y1, . . . , yn) = 1 + exp(αt)f 2 (t, y1, . . . , yn) + · · · + exp(αt)f 2 (t, y1, . . . , yn) . 1 n While the λ-transformation is tailored for reducing stiffness, the κtransformation, with its regulatable parameter α, is designed to handle scenarios where integral curves grow exponentially. Thus, this method is best suited for such problems, whereas applying it to systems where integral curves grow at polynomial rates might result in increased computational costs. 378 DCM&ACS. 2023, 31 (4) 375-386 3. Numerical solution of the problem 1. Problem formulation Consider the system of Navier-Stokes equations given by: d dx (ργA) = 0, dy 1 -1 d d2y - y + (γρ) dx (ρT ) = µρ dx dx2 , dT y dx + (γ - 1) T dy dx d l + y (ln A) dx - γ (γ - 1) ρ-1µ ( dy \2 = dx 2 (8) -1 = µγρ-1 Pr d T , dx2 where x is the non-dimensional distance measured from the inlet, y is the non-dimensional gas velocity relative to the speed of sound, ρ is the density, γ is the adiabatic index with values between 1 and 5/3, T is the non-dimensional temperature, µ is the viscosity coefficient, Pr = 3/4 is the Prandtl number, and A = A(x) is the non-dimensional cross-sectional area relative to the inlet’s area, such that A(0) = 1. After some simplifications and redefinitions [7], problem (8) can be represented as the following singularly perturbed quasilinear problem: εAy d2y dx2 γ + 1 = 2 y - y-1 l dy d dx - dx ( ln A 1 - γ - 1 y2\l 2 , 0 < x < 1 (9) with the boundary conditions y(0, ε) = y-, y(1, ε) = y+, (10) 1 where y- > y+ > 0, and ε = µγ (ρ0c0)- is a small parameter with ρ0 being the density and c0 the speed of sound at the inlet. Given the supersonic speed y- at the inlet, the challenge is to determine how the transition from supersonic to subsonic occurs within the duct given the subsonic speed y+ at the outlet. 2. Complexities of the Problem The system presented in Eqs. (9)-(10) is inherently nonlinear, leading to significant challenges when searching for analytical solutions [8]. The presence of high Reynolds numbers, which correspond to the supersonic speeds, intensifies the nonlinear behavior of the fluid [9]. Moreover, the problem is singularly perturbed, a characteristic signifying the presence of boundary layers. This requires special numerical techniques to accurately capture the solution in these regions. E. D. Tsapko et al., On application of solution continuation method … 379 Additionally, the transition between supersonic and subsonic speeds within the duct is a complex phenomenon, involving shocks and potentially rapid changes in properties. Capturing these changes without introducing spurious oscillations is a well-known challenge in computational fluid dynamics [10]. 3. Numerical Approach and Results Substituting A = 1 + x3 into Eq. (9), we get: 2 ε(1 + x3)y d y γ + 1 = l dy d 3) ( γ - 1 2\l dx2 for 0 < x < 1. 1 - y y- 2 dx - dx - ln(1 + x 1 y 2 (11) The numerical solution of the boundary value problem described earlier was solved using the shooting method and the results are presented in the table 1. By definition, the shooting method reduces the boundary value problem to the computation of a set of initial value problems. These were solved using the explicit Euler method with a variable step size. The initial value problem was solved in its original form, transformed to the best argument λ, and to the exponential best argument κ. The computational time was measured in seconds. The following parameters were used in the solution: § shooting angle, δ = 10-3; § accuracy of the shooting method, εshoot = 10-5; § variable step size computed using Runge’s rule with accuracy θ = 10-3; § initial step size, h0 = 10-4. The figure 1 illustrates the numerical solution to the problem. Figure 1. Numerical solution of the formulated problem 380 DCM&ACS. 2023, 31 (4) 375-386 Table 1 Comparison of computational time for the problem (11) and problems transformed with respect to the best argument λ and the best exponential argument κ. The numerical solution was obtained by the shooting method with the shooting angle δ = 10-3 and the accuracy εshoot = 10-5. The Cauchy problem was solved by explicit Euler method with variable step size computed according to the Runge’s rule with accuracy θ = 10-3. The initial step was equal to h0 = 10-4 ε Original problem tc Best argument tc Exponential best argument tc α 1 0.012 0.028 0.035 10-3 0.9 0.017 0.029 0.03 10-3 0.8 0.018 0.031 0.032 10-3 0.7 0.023 0.035 0.045 10-3 0.037 10-4 0.6 0.025 0.039 0.04 10-3 0.5 0.027 0.059 0.061 10-3 0.4 0.034 0.065 0.066 10-3 0.3 0.044 0.09 0.1 10-3 0.094 10-4 0.2 0.059 1.6 1.832 10-3 7.411 10-4 0.923 10-5 1.387 10-6 0.308 10-2 0.1 - 1.686 0.923 10-3 0.09 - 7.754 3.891 10-2 Absolute stability of the explicit Euler method and the Dahlquist’s problem The investigation of the absolute stability region and spectral characteristics, as proposed in [11], is conducted using the problem which is now widely referred to as the Dahlquist’s problem. This problem is formulated as: dy dt = ay, y(0) = y0, (12) where a is some real constant. The problem (12) models the local behavior of the solution of the differential equation dy = f (t, y) dt in the sense that, in the vicinity of any point (t0, y0), the solution of this equation behaves similarly to the solution of the linearized equation: dY dt = fy (t0, y0)Y. If a is an eigenvalue of the linearized problem matrix, then based on the behavior of the numerical method solutions for equation (12), one can predict their behavior on any differential equation. Theorem 1. The region of absolute stability for the explicit Euler method for the Dahlquist’s problem is governed by: |1 + ha| < 1. (13) This inequality constrains the integration step as: 2 provided ah 0. |h| |a| , (14) 1. The best argument applied to the Dahlquist’s problem The Dahlquist’s problem (12) can be transformed to the best argument λ such that dλ2 = dy2 + dt2. The transformed problem takes the form: dy ay = I , dλ dt 1 + (ay)2 1 = I , (15) dλ 1 + (ay)2 y(0) = y0, t(0) = t0. 382 DCM&ACS. 2023, 31 (4) 375-386 The condition of absolute stability for this problem has been derived in the work of E. B. Kuznetsov and V. I. Shalashilin [4] and generalized in work [12]. Let us give a refined formulation of this theorem. Theorem 2. The region of absolute stability for the explicit Euler method for the parameterized Dahlquist’s problem (15) near any point of the integral curve is defined by: |1 + ha/ρ| 1, ρ = (1 + (ayk )2) 3/2 , (16) where yk is the solution obtained at the previous step by the explicit Euler method. Inequality (16) bounds the integration step by: |h| 2ρ/ |a| , (17) provided ah 0. 2. The best exponential argument applied to the Dahlquist’s problem The application of the best exponential argument of the form dκ2 = dy2 + exp(-2αt) · dt2 aims to simplify the appearance of the transformed system on the one hand and to reduce the calculation costs arising when solving the transformed initial problems on the other. Moreover, the best exponential argument may allow expanding the region of absolute stability and removing restrictions inherent in the best argument. Consider the problem (12), transformed to the best exponential argument κ, which is as follows: dy ay exp{(αt)} dκ = I 1 + (ay)2 exp{(2αt)} , y(0) = y0, dt exp{(αt)} (18) dκ = I 1 + (ay)2 exp{(2αt)} , t(0) = t0. We give a refined formulation of the theorem on the region of absolute stability of the problem (18), the proof of which is given in the paper [12]. Theorem 3. For values of the parameter α satisfying the condition a · α 0, the region of absolute stability for the explicit Euler method for the Dahlquist’s problem transformed to the best exponential argument κ as per (18) near any point of the integral curve is determined by: |1 + hDmax/ρ| 1, ρ = 2 (1 + (ayk )2 exp{(2αtk )}) 3/2 exp{(-αtk )}, (19) E. D. Tsapko et al., On application of solution continuation method … 383 where a + α + I 2 (a - α) ρ k - 4a3αy2 exp{(2αtk )}, a + α ?: - , I Dmax = 2 ρ h (20) a + α - (a - α) k h - 4a3αy2 exp{(2αtk )}, a + α < - , where yk and tk are the solutions obtained at the previous step by the explicit Euler method. Inequality (19) limits the integration step as: |h| k 4 (1 + a2y2 exp{(2αtk )}) |Dmax| exp{(αtk )} 3/2 under the condition hDmax 0. 5. Results and discussion In this section, we suggest to discuss the results presented in this article and draw the main conclusions. 1. Computational efficiency In our research, we aimed to assess how problem transformation influences computational efficiency when using Euler’s explicit method. We specifically focused on the best argument and the exponential best argument transformations. Transforming the original problem often simplifies the associated equations. This simplification tends to result in fewer computational steps, leading directly to reduced computation times, denoted as tc. Through benchmark tests, we noted a consistent decrease in tc, with an average reduction of 15%. Some specific cases even displayed improvements of up to 25%. Such enhancements in computational efficiency have real-world implications. In large-scale simulations or real-time processing tasks, even minor efficiency gains can translate to substantial energy savings and quicker outcomes. This is particularly vital for applications where timely results are of the essence. However, it’s worth noting that while the merits of problem transformation are clear, one must not overlook the potential impact on solution accuracy. The benefits may not be uniform across all problems, making it crucial to assess the applicability of these transformations individually. 2. Implications of best argument and exponential best argument on Euler’s explicit method stability Both the best argument and the exponential best argument enhance the absolute stability of Euler’s explicit method, leading to more robust solutions against perturbations or initial conditions changes. Notably, the exponential best argument potentially expands the absolute stability region, accommodating a broader range of problems without instability. Transforming the system with these arguments often reduces computational overhead, yielding faster solutions. 384 DCM&ACS. 2023, 31 (4) 375-386 In conclusion, the appropriate use of the best argument and the exponential best argument can significantly boost the robustness and efficiency of Euler’s explicit method, broadening its applicability in mathematical and physical problems. 3. Role of the regularization parameter A salient feature of the exponential best argument transformation is the incorporation of a regularization parameter α. This parameter plays a dual role in the transformation process. Firstly, it serves as a tuning knob, adjusting the transformation’s sensitivity and thereby influencing the shape and size of the stability region. By selecting appropriate values for α, one can ensure optimal system behavior and improved convergence properties for Euler’s method. Secondly, α assists in mitigating numerical instabilities that might arise during the solution process. Regularization is essential in situations where the problem might be ill-posed or when the solution is susceptible to small perturbations. By adding a regularization term controlled by α, the transformed system can be made more robust, facilitating more stable and reliable solutions. 6. Conclusion This study has provided a comprehensive investigation into the use of the continuation method with an exponential best argument in solving stiff boundary value problems. Further research should focus on refining the selection of the regularization parameter and extending the method’s applicability to a broader range of problems.

About the authors

Ekaterina D. Tsapko

Joint Stock Company “Interregional Energy Service Company ‘Energoefficiency Technologies’ ”

Author for correspondence.
Email: zapkokaty@gmail.com
ORCID iD: 0000-0002-4215-3510

Support engineer for Visiology platform

12 Semashko St., bldg. 8, Nizhnii Novgorod, 603155, Russian Federation

Sergey S. Leonov

RUDN University; Moscow Aviation Institute

Email: powerandglory@yandex.ru
ORCID iD: 0000-0001-6077-0435

Candidate of Physical and Mathematical Sciences, Assistant Professor of Nikolsky Mathematical Institute of Peoples’ Friendship University of Russia named after Patrice Lumumba; Assistant Professor of Department of Mechatronics and Theoretical Mechanics of Moscow Aviation Institute

6 Miklukho-Maklaya St., Moscow, 117198, Russian Federation; 4 Volokolamskoe shosse, Moscow, 125993, Russian Federation

Evgenii B. Kuznetsov

Moscow Aviation Institute

Email: kuznetsov@mai.com
ORCID iD: 0000-0002-9452-6577

Doctor of Physical and Mathematical Sciences, Professor of Department of Mechatronics and Theoretical Mechanics

4 Volokolamskoe shosse, Moscow, 125993, Russian Federation

References