28 Analysis of Multilayer Methods of Building Approximate Solutions of Differential Equations in the Context of Solving a Homogeneous Duffing Equation* Serafim K. Boyarsky [0000-0003-0093-3719], Tatiana V. Lazovskaya [0000-0002-3324-6213] and Dmitry A. Tarkhov [0000-0002-9431-8241] Peter the Great St. Petersburg Polytechnic University dtarkhov@gmail.com Abstract. Multilayer methods are an alternative approach to building the approx- imate analytical solution of differential equations. This paper presents the study of the results obtained by the implementation of our modifications of acknowl- edged implicit and explicit numerical methods. The homogeneous Duffing equa- tion is of practical interest for modeling nonlinear oscillations and considered as a model equation. The accuracy of the obtained solutions is compared. It is shown that moving an initial point can significantly increase the accuracy of the solu- tions. Keywords: differential equations, numerical methods, analytical methods, mul- tilayer models, Duffing equation. 1 Introduction Modeling the behavior of many real objects often reduces to initial or boundary value problems for differential equations. In practice, the analytical solution of differential equations usually cannot be built, therefore, a numerical approximate solution is often sought. But such solutions are not clear enough. It is complicated to use it for studying the effect of changing the parameters of the original problem or adjust it to the behavior of the simulated object using the results of observation. Another well-known approach is building the approximate analytical solution. A lot of different approaches to finding it has been developed. There are various asymptotic methods, series expansion meth- ods, etc. [1]. Classic perturbation methods [2] are quite versatile but, like other non- linear analytical methods, they have significant limitations and restrictions. The quality of the solution may directly depend on the choice of the parameter by the researcher. In recent decades, new methods have appeared and old ones have been improved [3]. * Copyright 2021 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). 29 But often the higher-order analytical approximations require analytical solving of sys- tems of complex nonlinear algebraic equations or additional restrictions on the param- eters or the function [4-7]. Other approaches imply building a functional approximation in the form of broken lines or splines based on points of the numerical solution. In this article, we consider methods for building multilayer models that allow us to obtain an approximate analytical solution based on classical numerical methods. We compare the solutions built according to recurrence formulas of various numerical methods and in- vestigate the possibility of increasing the accuracy of the obtained solutions on the base of the initial point moving. 2 General Description of Multilayer Methods The essence of our approach is to apply the well-known recurrence formulas for the numerical integration of differential equations to an interval with a variable right end- point. [8-12]. The result is an approximate analytical solution in the form of a function of this endpoint. Consider the Cauchy problem for a system of ordinary differential equations 𝑦 β€² = 𝑓 (π‘₯, 𝑦) { π‘₯ ∈ 𝑅, 𝑦 ∈ π‘…π‘š (1) 𝑦(π‘₯0 ) = 𝑦0 The search for a solution is carried out on the interval [π‘₯0 , π‘₯0 + π‘Ž]. According to the main idea of our approach, we use well-known methods for the numerical solution of the problem (1) to an interval with a variable right endpoint π‘₯ ∈ [π‘₯0 , π‘₯0 + π‘Ž]. We choose a partition 𝑇𝑛 (π‘₯) of the interval [π‘₯0 , π‘₯ ] into n subintervals π‘₯0 < π‘₯1 < β‹― < π‘₯π‘˜ < β‹― < π‘₯𝑛 = π‘₯, β„Žπ‘˜ = π‘₯π‘˜+1 βˆ’ π‘₯π‘˜ . By applying the formula π‘¦π‘˜+1 = π‘¦π‘˜ + 𝐹(𝑓, β„Žπ‘˜ , π‘₯π‘˜ , π‘¦π‘˜+1 , π‘¦π‘˜ , π‘¦π‘˜βˆ’1 , … , 𝑦0 ) (2) n times, we obtain an approximate solution 𝑦𝑛 (π‘₯). The operator 𝐹 defines a specific method we modify as described above. The result is a function defined on the interval [π‘₯0 , π‘₯0 + π‘Ž]. We can replace an interval [π‘₯0 , π‘₯0 + π‘Ž] with [π‘₯0 βˆ’ 𝑏, π‘₯0 + π‘Ž]. 3 Multilayer Methods in the Context of Solving the Duffing Equation As the model task, we consider the homogeneous Duffing equation with constant coef- ficients [4-5, 13-16]. 𝑦 β€²β€² + π‘Žπ‘¦ + 𝑏𝑦 3 = 0 (3) 𝑦(0) = 𝑦0 , 𝑦 β€² (0) = 𝑦1 Higher-order differential equations can always be reduced to the form (1) by increas- ing the dimension of the system. In our case, it is easy to make a replacement 𝑦 = 𝑣, 𝑦 β€² = 𝑒: 30 𝑣′ = 𝑒 𝑒′ = βˆ’π‘Žπ‘£ βˆ’ 𝑏𝑣 3 (4) 𝑣 (0) = 𝑦0 , 𝑒(0) = 𝑦1 The search for a solution is carried out on the interval [0,3], the initial conditions and parameters of the equation are as follows π‘Ž = 1, 𝑏 = 1, 𝑦0 = 1, 𝑦1 = 1 For simplicity, the partition 𝑇𝑛 (π‘₯) is considered uniform for each method, namely β„Žπ‘˜ = (π‘₯ βˆ’ π‘₯0 )⁄𝑛. 3.1 Explicit Methods Euler's method. The simplest numerical method for which the operator F has the form 𝐹 (𝑓, β„Žπ‘˜ , π‘₯π‘˜ , π‘¦π‘˜ ) = β„Žπ‘˜ 𝑓 (π‘₯π‘˜ , π‘¦π‘˜ ) Refined Euler's method. Another well-known numerical method in which we used the formula 𝐹 (𝑓, β„Žπ‘˜ , π‘₯π‘˜ , π‘¦π‘˜ , π‘¦π‘˜βˆ’1 ) = 2β„Žπ‘˜ 𝑓(π‘₯π‘˜ , π‘¦π‘˜ ) + π‘¦π‘˜βˆ’1 βˆ’ π‘¦π‘˜ . In this case, to start the algorithm, the expression β„Ž1 β„Ž1 𝑦1 = 𝑦0 + β„Ž1 𝑓 (π‘₯0 + , 𝑦0 + 𝑓 (π‘₯0 , 𝑦0 )) 2 2 is used. y10 x yx y10 x yx y y 3 3 2 2 1 1 x x 0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0 1 1 2 2 3 3 Fig. 1. The plot of the exact solution of Fig. 2. The plot of the exact solution of problem (3) and the approximate solution problem (3) and the approximate solution built by our modification of the Euler built by our modification of the refined method in the case of n=10 Euler method in the case of n=10 The modified Euler method. This method works under the formula β„Žπ‘˜ β€² 𝐹 (𝑓, β„Žπ‘˜ , π‘₯π‘˜ , π‘¦π‘˜ ) = β„Žπ‘˜ [𝑓(π‘₯π‘˜ , π‘¦π‘˜ ) + (𝑓 (π‘₯ , 𝑦 ) + 𝑓𝑦′ (π‘₯π‘˜ , π‘¦π‘˜ )𝑓(π‘₯π‘˜ , π‘¦π‘˜ ))]. 2 π‘₯ π‘˜ π‘˜ 31 Second-order Runge-Kutta method: β„Žπ‘˜ 𝐹 (𝑓, β„Žπ‘˜ , π‘₯π‘˜ , π‘¦π‘˜ ) = β„Žπ‘˜ 𝑓 (π‘₯π‘˜ , π‘¦π‘˜ + 𝑓(π‘₯π‘˜ , π‘¦π‘˜ )). 2 y10 x yx y10 x yx y y 3 3 2 2 1 1 x x 0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0 1 1 2 2 3 3 Fig. 3. The plot of the exact solution of Fig. 4. The plot of the exact solution of prob- problem (3) and the approximate solution built lem (3) and the approximate solution built by by our modification of the modified Euler our modification of the Runge-Kutta method method in the case of n=10 in the case of n=10 StΓΆrmer Method. Since initially, the Duffing equation is a second-order differential equation, we can apply the StΓΆrmer method. In this case π‘¦π‘˜+1 = 2π‘¦π‘˜ βˆ’ π‘¦π‘˜βˆ’1 + β„Žπ‘˜2 𝑓 (π‘₯π‘˜ , π‘¦π‘˜ ). This method requires 𝑦1 . We calculated it approximately by the Taylor formula 𝑦 β€² (π‘₯ 0 ) 𝑦 β€²β€² (π‘₯0 ) 2 𝑦′′′(π‘₯0 ) 3 𝑦1 = 𝑦0 + β„Ž1 + β„Ž1 + β„Ž1 , 1! 2! 3! where 𝑦 β€²β€² (π‘₯0 ) and 𝑦′′′(π‘₯0 ) are easily obtained from differential equation (3). y5 x yx y10 x yx y y 3 3 2 2 1 1 x x 0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0 1 1 2 2 3 3 Fig. 5. The plot of the exact solution of prob- Fig. 6. The plot of the exact solution of prob- lem (3) and the approximate solution built by lem (3) and the approximate solution built by our modification of the StΓΆrmer Method in the our modification of the StΓΆrmer Method in the case of n=5 case of n=10 32 3.2 Implicit Methods Implicit methods are applicable if the equation π‘¦π‘˜+1 = π‘¦π‘˜ + 𝐹(𝑓, β„Žπ‘˜ , π‘₯π‘˜ , π‘¦π‘˜+1 , π‘¦π‘˜ , π‘¦π‘˜βˆ’1 , … , 𝑦0 ) is solvable for π‘¦π‘˜+1 . About our problem, this means solving the cubic equation at each step. In some cases, it is advisable not to look for the exact solution of this cubic equation but to use the approximate solving methods. We use one step of the Newton method. The justification of this approach is demonstrated below. The implicit Euler method. The operator F for an implicit method is as follows 𝐹 (𝑓, β„Žπ‘˜ , π‘₯π‘˜ , π‘¦π‘˜+1 ) = β„Žπ‘˜ 𝑓 (π‘₯π‘˜+1 , π‘¦π‘˜+1 ). Substituting this expression in (2) we obtain the system π‘£π‘˜+1 = π‘£π‘˜ + β„Žπ‘˜ π‘’π‘˜+1 { 3 ). (5) π‘’π‘˜+1 = π‘’π‘˜ + β„Žπ‘˜ (βˆ’π‘Žπ‘£π‘˜+1 βˆ’ π‘π‘£π‘˜+1 This system allows the exact expression π‘£π‘˜+1 , π‘’π‘˜+1 in terms of π‘£π‘˜ , π‘’π‘˜ , β„Žπ‘˜ : 𝑣 = 𝐡1 (π‘£π‘˜ , π‘’π‘˜ , β„Žπ‘˜ ) { π‘˜+1 . π‘’π‘˜+1 = 𝐡2 (π‘£π‘˜ , π‘’π‘˜ , β„Žπ‘˜ ) Then we can use formula π‘£π‘˜+1 π‘£π‘˜ 𝐡 (𝑣 , 𝑒 , β„Ž ) (𝑒 ) = (𝑒 ) + ( 1 π‘˜ π‘˜ π‘˜ ). π‘˜+1 π‘˜ 𝐡2 (π‘£π‘˜ , π‘’π‘˜ , β„Žπ‘˜ ) to perform computations. The result of such solving is presented below. On the other hand, the solution of system (5) can be obtained by using one step of the Newton method. In this case, we obtain another expression, 𝑣 = 𝑁1 (π‘£π‘˜ , π‘’π‘˜ , β„Žπ‘˜ ) { π‘˜+1 . π‘’π‘˜+1 = 𝑁2 (π‘£π‘˜ , π‘’π‘˜ , β„Žπ‘˜ ) y5 x yx y5 x yx y y 3 3 2 2 1 1 x x 0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0 1 1 2 2 3 3 Fig. 7. The plot of the exact solution of prob- Fig. 8. The plot of the exact solution of prob- lem (3) and the approximate solution built by lem (3) and the approximate solution built by our modification of the implicit Euler method our modification of the implicit Euler method in the case of the exact solution of (5) and n=5 in the case of the approximate solution of (5) (one step of the Newton method) and n=5 33 The maximum error in the first case is equal to 0.89163, in the second case it is 0.96993. As we can see, there is no significant change in the plot behavior. But when using the Newton method an expression is easier and therefore the complexity and time of calculations are lower than if we solve system (5) analytically. The solutions below are obtained using the Newton method. One-step Adams method. Another implicit method, the equation has the form β„Ž π‘¦π‘˜+1 = π‘¦π‘˜ + 2π‘˜ (𝑓(π‘₯π‘˜ , π‘¦π‘˜ ) + 𝑓 (π‘₯π‘˜+1 , π‘¦π‘˜+1 )). y5 x yx y5 x yx y y 3 3 2 2 1 1 x x 0.5 1.0 1.5 2.0 2.5 3.0 0.5 1.0 1.5 2.0 2.5 3.0 1 1 2 2 3 3 Fig. 9. The plot of the exact solution of prob- Fig. 20. The plot of the exact solution of prob- lem (3) and the approximate solution built by lem (3) and the approximate solution built by our modification of the One-step Adams our modification of the One-step Adams method in the case of the exact solution of the method in the case of the approximate solution cubic equation and n=5 of the cubic equation (one step of the Newton method) and n=5 The maximum error in the case of the exact solution according to the cubic equation is equal to 1.11217 and in the case of the approximate solution, it is 1.14777. 3.3 Computational results The computational results of all methods are presented in Table1. Each method corre- sponds to its maximum error in the specified interval. We compare the error on the intervals [0,1.5] and [0,3]. Table 1. The maximum errors of the studied methods with the number of layers n=2 and n=5. Number of iterations n 5 10 Interval [0, 1.5] [0, 3] [0, 1.5] [0, 3] Euler method 0.17474 7.11176 0.14765 1.19331 Refined Euler method 0.10981 11.70914 0.022633 0.19536 Modified Euler method 0.1148 0.90283 0.026272 0.13579 Runge-Kutta method 0.13924 1.90271 0.032397 0.19168 StΓΆrmer Method 0.023506 0.18407 0.0058227 0.034254 Implicit Euler method 0.30049 0.96993 0.16892 0.62051 One-step Adams method 0.77332 1.14778 0.19285 0.67787 34 Implicit methods have no advantages over explicit methods for this task. Of the explicit methods we examined, the most accurate is the StΓΆrmer method. 4 Initial Point Moving To improve the accuracy of the model, we investigated the following approach. Using the methods described above, we build a solution to the Cauchy problem (1) starting from the point π‘₯1 ∈ [π‘₯0 , π‘₯0 + π‘Ž], other than π‘₯0. The unknown initial condition 𝑦(π‘₯1 ) = 𝑦1 in this case is the parameter of the resulting solution 𝑦𝑛 (π‘₯, 𝑦|1). This parameter can be determined from the equation 𝑦𝑛 (π‘₯0 , 𝑦|1) = 𝑦0 . From the computational results be- low, it follows that by moving an initial point in this way, it is possible to improve the solution on the interval. As an example, we consider the StΓΆrmer method, as the most accurate method of the previously considered. The table below shows the maximum error of the solution ob- tained by this method on the interval [0,1.5] in the case of moving the initial point from 1 zero in increments of . The number of layers n we took equal to 2 and 5. 10 Table 2. The maximum errors in the interval [0,1.5] in the case of moving the initial point. π‘₯1 n=2 n=5 π‘₯1 n=2 n=5 0 0.11273 0.023506 0.8 0.03919 0.0077468 0.1 0.10936 0.022624 0.9 0.062919 0.0058778 0.2 0.15743 0.024725 1 0.11639 0.012437 0.3 0.20216 0.027498 1.1 0.16437 0.019892 0.4 0.21231 0.028058 1.2 0.20177 0.027324 0.5 0.18527 0.024551 1.3 0.23007 0.034529 0.6 0.13004 0.017213 1.4 0.25373 0.042016 0.7 0.063503 0.0091472 1.5 0.27807 0.050909 As we can see, solutions with moving the initial point have significantly better accuracy on the interval than a conventional solution built starting from zero. y2 x yx y2 x yx y y 2.0 2.0 1.5 1.5 1.0 1.0 0.5 0.5 x x 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.5 0.5 1.0 1.0 Fig. 11. The plot of the exact solution of prob- Fig. 32. The plot of the exact solution of prob- lem (3) and the approximate solution built by lem (3) and the approximate solution built by our modification of the StΓΆrmer Method in the our modification of the StΓΆrmer Method in the case of n=2 and initial point π‘₯1 = 0.8 case of n=2 and initial point π‘₯1 = 0 35 Next, we select new starting points from the gap between the previous best results in 1 increments of 100. 1/100. From this list of models, the most accurate model can be cho- sen. Table 3. The maximum errors in the interval [0,1.5] in the case of moving the initial point in the interval [0.8,0.9] π‘₯1 n=2 n=5 0.8 0.03919 0.0077468 0.81 0.036768 0.007522 0.82 0.034308 0.0072874 0.83 0.031835 0.0070447 0.84 0.034117 0.0067956 0.85 0.038495 0.006542 0.86 0.043071 0.0062857 0.87 0.047822 0.0060285 0.88 0.052729 0.0057722 0.89 0.057768 0.0055186 0.9 0.062919 0.0058778 Thus, the approach with moving an initial point allowed us to significantly reduce the deviation of the solution on the interval. The plots below illustrate the difference be- tween models with a selected initial point and models built starting from zero. y2 x y x y2 x y x y y 0.14 0.14 0.12 0.12 0.10 0.10 0.08 0.08 0.06 0.06 0.04 0.04 0.02 0.02 x x 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Fig. 43. The plot of the approximate solution Fig. 54. The plot of the approximate solution error module in the case of 𝑛 = 2 and initial error module in the case of 𝑛 = 2 and initial point π‘₯1 = 0.83 point π‘₯0 = 0 36 y5 x y x y5 x y x y y 0.025 0.025 0.020 0.020 0.015 0.015 0.010 0.010 0.005 0.005 x x 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Fig. 15. The plot of the approximate solution Fig. 16. The plot of the approximate solution error module in the case of 𝑛 = 5 and initial error module in the case of 𝑛 = 5 and initial point π‘₯1 = 0.89 point π‘₯0 = 0 5 Conclusions and Discussion The study allows us to draw the following conclusions: 1. The methods we have proposed allow us to construct an approximate solution of the Duffing equation in the form of a function with the required accuracy. 2. For model equations with parameters considered, implicit methods do not have sig- nificant advantages over explicit methods. Implicit methods make sense for those parameters when the task becomes stiff. 3. Moving an initial point lets us obtain an approximate analytical solution of the model task which is several times more accurate than a solution obtained without moving an initial point. Wherein an accuracy increases with the number of layers. 4. Our methods, without requiring additional assumptions, allow us to build parametric approximate analytical solutions [9] concerning the parameters of the original task. Acknowledgment This paper is based on research carried out with the financial support of the grant of the Russian Scientific Foundation (project β„–18-19-00474). References 1. Hairer E., Norsett S. P., Wanner G. Solving Ordinary Differential Equations I: Nonstiff Problem, Springer-Verlag, Berlin, 480 p. 2. Nayfeh, A. H. Problems in perturbation, John Wiley & Sons, 556 p. 3. He J.-H. Ξ‘ Review on Some New Recently Developed Nonlinear Analytical Techniques. International Journal of Nonlinear Sciences and Numerical Simulation 1, 51-70. 4. Wu B.S. et al. An analytical approximate technique for a class of strongly non-linear oscil- lators. International Journal of Non-Linear Mechanics 41, 766-774. 37 5. Farzaneh Y., Tootoonchi A. A. Global Error Minimization method for solving strongly non- linear oscillator differential equations, Computers and Mathematics with Applications, 59, 2887–2895. 6. Ismail G. M. An analytical coupled homotopy-variational approach for solving strongly non- linear differential equation Journal of the Egyptian Mathematical Society Volume 25, Issue 4, 434-437. 7. Rahman M. S., Hasan A. S. M. Z. Modified harmonic balance method for the solution of nonlinear jerk equations Results in Physics, Volume 8, 893-897. 8. Lazovskaya T. and Tarkhov D. Multilayer neural network models based on grid methods // Journal of Physics Conference Series: Materials Science and Engineering, Volume 158, Number 1. 9. Tarkhov D., Vasilyev А. Semi-empirical Neural Network Modeling and Digital Twins De- velopment Academic Press, Elsevier, 288 p. 10. Kaverzneva T. T., Malykhina G. F., Tarkhov D. A. From Differential Equations to Multi- layer Neural Network Models, International Symposium on Neural Networks ISNN 2019: Advances in Neural Networks – ISNN 2019 Lecture Notes in Computer Science book series (LNCS, volume 11554) 19-27. 11. Kuznetsov E. B., Leonov S. S., Tarkhov D. A. and Vasilyev A. N. Multilayer method for solving a problem of metals rupture under creep conditions, Thermal Science: vol. 23, suppl. 2, S575-S582. 12. Shemyakina T. Tarkhov D., Vasilyev A. and Velichko Y. Comparison of two neural network approaches to modeling processes in a chemical reactor, Thermal Science: vol. 23, suppl. 2, S583-S589. 13. Moon F. C. Chaotic Vibrations: An Introduction for Applied Scientists and Engineers, John Wiley & Sons, 316 p. 14. Astrakhov V.V., Koblyanskiy S.A., Shabunin A.V. Ostsilyator Duffinga: Uchebnoye po- sobiye dlya studentov vuzov (in Russian), 52 p. 15. Lakshmanan M., Murali K. Chaos in Nonlinear Oscillators: Controlling and Synchroniza- tion, World Scientific (World Scientific Series on Nonlinear Science Series A)., Vol. 13. 35- 90 (1996) 16. Kovacic I., Brennan M. J. The Duffing Equation: Nonlinear Oscillators and their Behaviour, John Wiley & Sons, 386 p.