=Paper=
{{Paper
|id=Vol-3309/short20
|storemode=property
|title=Computational Methods of Integration of Deterministic Systems, which are Described by Systems of Ordinary Differential Equations
|pdfUrl=https://ceur-ws.org/Vol-3309/short20.pdf
|volume=Vol-3309
|authors=Vladyslav Khaidurov,Tamara Tsiupii,Tetiana Zhovnovach,Artur Zaporozhets,Olena Kharchenko,Serhii Kharchenko
|dblpUrl=https://dblp.org/rec/conf/ittap/KhaidurovTZZKK22
}}
==Computational Methods of Integration of Deterministic Systems, which are Described by Systems of Ordinary Differential Equations==
https://ceur-ws.org/Vol-3309/short20.pdf
Computational Methods of Integration of Deterministic
Systems, which are Described by Systems of Ordinary
Differential Equations
Vladyslav Khaidurova, Tamara Tsiupiib, Tetiana Zhovnovachc, Artur Zaporozhetsa, Olena
Kharchenkod and Serhii Kharchenkoa
a
Institute of General Energy of NAS of Ukraine, 172, Antonovycha st., Kyiv, 03150, Ukraine
b
National University of Life and Environmental Sciences of Ukraine, 15, Heroyiv Oborony st., Kyiv, 03041,
Ukraine
c
Cherkasy Branch of the European University, Smilyanska st., 83, Cherkasy, 18000, Ukraine
d
National Aviation University, Liubomyra Huzara ave., 1, Kyiv, 03058, Ukraine
Abstract
The result of the work is obtaining difference schemes for dynamic systems, which are
described by systems of ordinary differential equations of the second order based on the
classical methods of Bossak, Newmark and the generalized α-method. For the developed
modifications of the methods, the corresponding software was developed in the MATLAB
system of applied mathematics for conducting numerical experiments and testing methods.
Keywords 1
Deterministic system, Newmark’s method, Bossack’s method, generalized α-method,
computational method.
1. Introduction
The development of computational methods for modeling dynamic systems occupies a significant
place in solving scientific and technical tasks in today's conditions. A number of requirements are put
forward to such methods, in particular, the methods that are developed should be stable and those that
obtain a numerical solution with minimal errors [1; 2]. The rapid development of computer technology
makes it possible to apply these methods in practice. Therefore, the work will demonstrate a tool for
obtaining basic mathematical dependencies based on classical methods for the analysis of nonlinear
systems, for example, systems described by Cauchy problems for second-order ordinary differential
equations.
2. The main part of the article
To analyze dynamic systems, let's take a classical nonlinear system in the form of a nonlinear
oscillator [3]. The nonlinear oscillator should be viewed as a generalized model suitable for describing
oscillatory phenomena in systems of different physical nature.
A nonlinear oscillator is a system whose dynamics is described by a second-order differential
equation of the general form [3; 4]:
ITTAP’2022: 2nd International Workshop on Information Technologies: Theoretical and Applied Problems, November 22–24, 2022,
Ternopil, Ukraine
EMAIL: allif0111@gmail.com (A. 1); ts.tamara19@gmail.com (A. 2); ts.tamara19@gmail.com (A. 3); a.o.zaporozhets@nas.gov.ua (A. 4);
olena80@ukr.net (A. 5); nanoavia@ukr.net (A. 6)
ORCID: 0000-0002-4805-8880 (A. 1); 0000-0003-2206-2897 (A. 2); 0000-0003-1037-4383 (A. 3); 0000-0002-0704-4116 (A. 4);
0000-0002-1311-8548 (A. 5); 0000-0001-9808-7607 (A. 6)
©️ 2021 Copyright for this paper by its authors.
Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
CEUR Workshop Proceedings (CEUR-WS.org)
𝑥 ′′ + 𝛿𝑥 ′ + 𝑓(𝑥) = 0, (1)
where 𝑥 – dynamic variable (or generalized coordinate), 𝑓(𝑥) – some non-linear function, 𝛿 –
dissipation parameter. Depending on the different values of this parameter, different transformations of
the phase portrait of the nonlinear oscillator occur. In order to form a complete view of the nonlinear
oscillator near the catastrophe of the assembly, it is necessary to draw a function 𝑓(𝑥) species
equality (1). In equation (1), the function 𝑓(𝑥) has the following form:
𝑓(𝑥) = 𝑎̄ 𝑥 + 𝑏 + 𝑥 3 . (2)
Equality (2) corresponds to an assembly type disaster. When varying parameter values 𝑎̄ and 𝑏 the
behavior of the system will also change. Therefore, the second type of disaster is called an "assembly"
disaster. In the behavior of the phase portrait curves, it is possible to observe both trajectories without
a jump, with smooth development, and with a jump in development.
The catastrophe theory explains the dependence of unstable forms of the phase portrait on the
number of control parameters. Yes, if we have only one control parameter, we can observe a fold
disaster. If we have two control parameters, as in our case 𝑎̄ і 𝑏, a build crash is observed.
Based on equalities (1) and (2), we can form a complete view of the nonlinear oscillator near the
assembly catastrophe:
𝑥 ′′ + 𝛿𝑥 ′ + 𝑎̄ 𝑥 + 𝑏 + 𝑥 3 = 0. (3)
In equation (3) 𝛿 is dissipation parameter, 𝑎̄ and 𝑏 are control parameters [5; 6].
2.1. Newmark’s method
The Newmark method is a numerical integration method for solving differential equations. It is
widely used in the calculations of many engineering, technical, chemical problems for the numerical
evaluation of dynamic characteristics. The method is named after a former professor of civil engineering
at Illinois State University who developed it in 1959 for use in structural dynamics.
To derive the formulas and analyze them, consider the Cauchy problem for an ordinary differential
equation of the second order. Its mathematical representation has the form:
𝑥 ′′ + 𝑔𝑥 ′ + 𝑎𝑥 + 𝑏 + 𝑥 3 = 0;
𝑥(𝑡0 ) = 𝑥0 , 𝑥′ (𝑡0 ) = 𝑥0′.
1
𝑈𝑛+1 = 𝑈𝑛 + ℎ𝑉𝑛 + ℎ2 ( − 𝛽) 𝐴𝑛 + ℎ2 𝛽𝐴𝑛+1 ,
2 (4)
𝑉𝑛+1 = 𝑉𝑛 + ℎ(1 − 𝛾)𝐴𝑛 + ℎ𝛾𝐴𝑛+1 ,
3
𝐴𝑛+1 + 𝑔𝑉𝑛+1 + 𝑎𝑈𝑛+1 + 𝑏 + 𝑈𝑛+1 = 0.
Since the system (4) is nonlinear, then the classical Seidel method can be used for its numerical
solution. For this system, Seidel’s method will look like this:
(𝑘+1) 1 (𝑘)
𝑈𝑛+1 = 𝑈𝑛 + ℎ𝑉𝑛 + ℎ2 ( − 𝛽) 𝐴𝑛 + ℎ2 𝛽𝐴𝑛+1 ,
2
(𝑘+1) (𝑘)
𝑉𝑛+1 = 𝑉𝑛 + ℎ(1 − 𝛾)𝐴𝑛 + ℎ𝛾𝐴𝑛+1 , (5)
(𝑘+1) (𝑘+1) (𝑘+1) (𝑘+1) 3
𝐴𝑛+1 = − (𝑔𝑉𝑛+1 + 𝑎𝑈𝑛+1 + 𝑏 + (𝑈𝑛+1 ) ),
𝑈0 = 𝑥0 , 𝑉0 = 𝑥0′, 𝐴0 = −(𝑔𝑉0 + 𝑎𝑈0 + 𝑏 + 𝑈03 ),
where 𝑘 is iteration number.
The second version of Newmark's method can be presented in the following form:
1
𝑈𝑛+1 = 𝑈𝑛 + ℎ𝑉𝑛 + ℎ2 ( − 𝛽) 𝐴𝑛 + ℎ2 𝛽𝐴𝑛+1 ,
2
𝑉𝑛+1 = 𝑉𝑛 + ℎ(1 − 𝛾)𝐴𝑛 + ℎ𝛾𝐴𝑛+1 ,
3 ).
𝐴𝑛+1 = −(𝑔𝑉𝑛+1 + 𝑎𝑈𝑛+1 + 𝑏 + 𝑈𝑛+1
Then we will have
1 3
𝑈𝑛+1 = 𝑈𝑛 + ℎ𝑉𝑛 − ℎ2 ( − 𝛽) (𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ2 𝛽(𝑔𝑉𝑛+1 + 𝑎𝑈𝑛+1 + 𝑏 + 𝑈𝑛+1 ),
2
3) 3
𝑉𝑛+1 = 𝑉𝑛 − ℎ(1 − 𝛾)(𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛 − ℎ𝛾(𝑔𝑉𝑛+1 + 𝑎𝑈𝑛+1 + 𝑏 + 𝑈𝑛+1 ),
So,
1 3
𝑈𝑛+1 (1 + ℎ2 𝑎𝛽) = 𝑈𝑛 + ℎ𝑉𝑛 − ℎ2 ( − 𝛽) (𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ2 𝛽(𝑔𝑉𝑛+1 + 𝑏 + 𝑈𝑛+1 ),
2
3 ).
𝑉𝑛+1 (1 + 𝑔ℎ𝛾) = 𝑉𝑛 − ℎ(1 − 𝛾)(𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ𝛾(𝑎𝑈𝑛+1 + 𝑏 + 𝑈𝑛+1
Finally, the formulas for mathematical calculation take the form:
1 3 ),
𝑈𝑛 + ℎ𝑉𝑛 − ℎ2 ( − 𝛽) (𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ2 𝛽(𝑔𝑉𝑛+1 + 𝑏 + 𝑈𝑛+1
𝑈𝑛+1 = 2
1 + ℎ2 𝑎𝛽
3
𝑉𝑛 − ℎ(1 − 𝛾)(𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ𝛾(𝑎𝑈𝑛+1 + 𝑏 + 𝑈𝑛+1 )
𝑉𝑛+1 = .
1 + 𝑔ℎ𝛾
Values 𝑈𝑛+1 and 𝑉𝑛+1 are calculated using computational methods. For example, Seidel's method
for this system of equations will have the form:
1 (𝑘) (𝑘) 3
𝑈𝑛 + ℎ𝑉𝑛 − ℎ2 ( − 𝛽) (𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ2 𝛽 (𝑔𝑉𝑛+1 + 𝑏 + (𝑈𝑛+1 ) ) ,
(𝑘+1) 2
𝑈𝑛+1 =
1 + ℎ2 𝑎𝛽
(6)
(𝑘+1) (𝑘+1) 3
𝑉𝑛 − ℎ(1 − 𝛾)(𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ𝛾 (𝑎𝑈𝑛+1 + 𝑏 + (𝑈𝑛+1 ) )
(𝑘+1)
𝑉𝑛+1 = .
1 + 𝑔ℎ𝛾
The third version of the Newmark method looks like this:
3 )
𝑉𝑛 − ℎ(1 − 𝛾)(𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ𝛾(𝑎𝑈𝑛+1 + 𝑏 + 𝑈𝑛+1
𝑉𝑛+1 = ;
1 + 𝑔ℎ𝛾
1 3 ), (7)
𝑈𝑛 + ℎ𝑉𝑛 − ℎ2 ( − 𝛽) (𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ2 𝛽(𝑔𝑉𝑛+1 + 𝑏 + 𝑈𝑛+1
𝑈𝑛+1 = 2 .
1 + ℎ2 𝑎𝛽
It should be noted that in (7) the second equation is substituted into the first equation. Next is the
value 𝑈𝑛+1 , which is then substituted into the equation:
3 )
𝑉𝑛 − ℎ(1 − 𝛾)(𝑔𝑉𝑛 + 𝑎𝑈𝑛 + 𝑏 + 𝑈𝑛3 ) − ℎ𝛾(𝑎𝑈𝑛+1 + 𝑏 + 𝑈𝑛+1
𝑉𝑛+1 = .
1 + 𝑔ℎ𝛾
In this case, Seidel's method solves only one equation [7; 8].
2.2. Generalized 𝜶-method
Formulas for solving a nonlinear oscillator near the assembly catastrophe using the generalized 𝛼-
method are derived similarly to Newmark's method. The difference is that the values of all variables at
𝑛 + 1 time steps are taken as averages, taking into account the average theorem.
The equality that we need to solve by the method of generalized α has the following form:
𝑥 ′′ + 𝛿𝑥 ′ + 𝑎̄ 𝑥 + 𝑏 + 𝑥 3 = 0. (8)
In equation (8), 𝛿 is dissipation parameter, 𝑎̅ and 𝑏 control parameters.
As before, let's mark 𝑥 ′′ = 𝑎, 𝑥 ′ = 𝑣, and let's complete the equality (8). Now let's write down the
system of differential equations at 𝑛 + 1 instant of time according to the method of generalized 𝛼.
𝑣𝑛+1 = 𝑣𝑛 + (1 − 𝛾)𝛥𝑡 𝑎𝑛 + 𝛾𝛥𝑡 𝑎𝑛+1 ,
1
𝑥𝑛+1 = 𝑥𝑛 + 𝛥𝑡 𝑣𝑛 + ( − 𝛽) 𝛥𝑡 2 𝑎𝑛 + 𝛽𝛥𝑡 2 𝑎𝑛+1 , (9)
2
3
𝑎𝑛+1−𝛼𝑚 + 𝛿𝑣𝑛+1−𝛼𝑓 + 𝑎̄ 𝑥𝑛+1−𝛼𝑓 + 𝑏 + 𝑥𝑛+1−𝛼 = 0,
{ 𝑓
where 𝑎𝑛+1−𝛼𝑚 , 𝑣𝑛+1−𝛼𝑓 , 𝑥𝑛+1−𝛼𝑓 are equal to:
𝑎𝑛+1−𝛼𝑚 = (1 − 𝛼𝑚 )𝑎𝑛+1 + 𝛼𝑚 𝑎𝑛 , 𝑣𝑛+1−𝛼𝑓 = (1 − 𝛼𝑓 )𝑣𝑛+1 + 𝛼𝑓 𝑣𝑛,
(10)
𝑥𝑛+1−𝛼𝑓 = (1 − 𝛼𝑓 )𝑥𝑛+1 + 𝛼𝑓 𝑥𝑛.
Moreover, the values of all constants are equal to:
1 1 2 2𝜌∞ − 1 𝜌∞
𝛾= − 𝛼𝑚 + 𝛼𝑓 , 𝛽 = (1 − 𝛼𝑚 + 𝛼𝑓 ) , 𝛼𝑚 = , 𝛼𝑓 = .
2 4 𝜌∞ + 1 𝜌∞ + 1
With the condition that 0 ≤ 𝜌∞ ≤ 1.
Let's rewrite (9) considering (10):
𝑣𝑛+1 = 𝑣𝑛 + (1 − 𝛾)𝛥𝑡 𝑎𝑛 + 𝛾𝛥𝑡 𝑎𝑛+1
1
𝑥𝑛+1 = 𝑥𝑛 + 𝛥𝑡𝑣𝑛 + ( − 𝛽) 𝛥𝑡 2 𝑎𝑛 + 𝛽𝛥𝑡 2 𝑎𝑛+1
2 (11)
(1 − 𝛼𝑚 )𝑎𝑛+1 + 𝛼𝑚 𝑎𝑛 + 𝛿(1 − 𝛼𝑓 )𝑣𝑛+1 + 𝛿𝛼𝑓 𝑣𝑛 +
3 3
{+𝑎̄ (1 − 𝛼𝑓 )𝑥𝑛+1 + 𝛼𝑓 𝑥𝑛 + 𝑏 + (1 − 𝛼𝑓 )𝑥𝑛+1 + 𝛼𝑓 𝑥𝑛 = 0
From equality (11) we see that 𝑥 is a vector, and therefore velocity and acceleration are also vectors.
We obtained a closed system of 3𝑘 algebraic equations, with 3𝑘 unknowns.
Everything that is on the 𝑛 + 1 step is unknown and it is necessary to look for their value. To do
this, we move all variables 𝑛 + 1 steps to the left, and n steps to the right:
𝛾𝛥𝑡 𝑎𝑛+1 − 𝑣𝑛+1 = −𝑣𝑛 − (1 − 𝛾)𝛥𝑡𝑎𝑛
1
𝛽𝛥𝑡 2 𝑎𝑛+1 − 𝑥𝑛+1 = −𝑥𝑛 − 𝛥𝑡𝑣𝑛 − ( − 𝛽) 𝛥𝑡 2 𝑎𝑛
2 (12)
3
(1 − 𝛼𝑚 )𝑎𝑛+1 + 𝛿(1 − 𝛼𝑓 )𝑣𝑛+1 + 𝑎̄ (1 − 𝛼𝑓 )𝑥𝑛+1 + (1 − 𝛼𝑓 )𝑥𝑛+1 =
3
{ = −𝛼𝑚 𝑎𝑛 − 𝛿𝛼𝑓 𝑣𝑛 − 𝛼𝑓 𝑥𝑛 − 𝑏 − 𝛼𝑓 𝑥𝑛
We see that equality (12) is written in vector form. Therefore, to simplify the notations, we introduce
the vectors:
1
𝐴𝑛+1 = −𝑣𝑛 − (1 − 𝛾)𝛥𝑡 𝑎𝑛 , 𝐵𝑛+1 = −𝑥𝑛 − 𝛥𝑡 v𝑛 − ( − 𝛽) 𝛥𝑡 2 𝑎𝑛 ,
2
𝐶𝑛+1 = −𝛼𝑚 𝑎𝑛 − 𝛿𝛼𝑓 𝑣𝑛 − 𝛼𝑓 𝑥𝑛 − 𝑏 − 𝛼𝑓 𝑥𝑛3.
𝛾𝛥𝑡 𝑎𝑛+1 − 𝑣𝑛+1 = 𝐴𝑛+1
{ 𝛽𝛥𝑡 2 𝑎𝑛+1 − 𝑥𝑛+1 = 𝐵𝑛+1 (13)
3
(1 − 𝛼𝑚 )𝑎𝑛+1 + 𝛿(1 − 𝛼𝑓 )𝑣𝑛+1 + 𝑎̄ (1 − 𝛼𝑓 )𝑥𝑛+1 + (1 − 𝛼𝑓 )𝑥𝑛+1 = 𝐶𝑛+1
Next, we divide the first equality (13) by 𝛾𝛥𝑡, and the second equality – by 𝛽𝛥𝑡 2
1 1
𝑎𝑛+1 − 𝑣𝑛+1 = 𝐴 ,
𝛾𝛥𝑡 𝛾𝛥𝑡 𝑛+1
1 1 (14)
𝑎𝑛+1 − 𝑥𝑛+1 = 𝐵 ,
𝛽𝛥𝑡 2 𝛽𝛥𝑡 2 𝑛+1
3
{(1 − 𝛼𝑚 )𝑎𝑛+1 + 𝛿(1 − 𝛼𝑓 )𝑣𝑛+1 + 𝑎̄ (1 − 𝛼𝑓 )𝑥𝑛+1 + (1 − 𝛼𝑓 )𝑥𝑛+1 = 𝐶𝑛+1 .
Now subtract the first from the second equality (14), and from the third the first multiplied by
(1 − 𝛼𝑚 ).
1 1 1 1
𝑣𝑛+1 − 𝑥𝑛+1 = 𝐵𝑛+1 − 𝐴 ,
𝛾𝛥𝑡 𝛽𝛥𝑡 2 𝛽𝛥𝑡 2 𝛾𝛥𝑡 𝑛+1
(15)
1 − 𝛼𝑚 3 1 − 𝛼𝑚
(𝛿(1 − 𝛼𝑓 ) + ) 𝑣𝑛+1 + 𝑎̄ (1 − 𝛼𝑓 )𝑥𝑛+1 + (1 − 𝛼𝑓 )𝑥𝑛+1 = 𝐶𝑛+1 − 𝐴𝑛+1 .
{ 𝛾𝛥𝑡 𝛾𝛥𝑡
For simplification, we introduce 𝐷 = (𝛿(1 − 𝛼𝑓 ) + (1 − 𝛼𝑚 )⁄(𝛾𝛥𝑡)). Rewrite (15):
𝛾 𝛾
𝑣𝑛+1 − 𝑥𝑛+1 = 𝐵 − 𝐴𝑛+1 ,
𝛽𝛥𝑡 𝛽𝛥𝑡 𝑛+1
1 − 𝛼𝑚 (16)
3
𝐷𝑣𝑛+1 + 𝑎̄ (1 − 𝛼𝑓 )𝑥𝑛+1 + (1 − 𝛼𝑓 )𝑥𝑛+1 = 𝐶𝑛+1 − 𝐴𝑛+1 .
{ 𝛾𝛥𝑡
Next, we subtract the first multiplied by 𝐷 from the second equality (16).
𝛾 3
𝑎̄ (1 − 𝛼𝑓 )𝑥𝑛+1 + 𝐷𝑥 + (1 − 𝛼𝑓 )𝑥𝑛+1 =
𝛽𝛥𝑡 𝑛+1
1 − 𝛼𝑚 𝛾 (17)
= 𝐶𝑛+1 − 𝐴𝑛+1 + 𝐷𝐴𝑛+1 − 𝐷𝐵 .
𝛾𝛥𝑡 𝛽𝛥𝑡 𝑛+1
In order to avoid division by zero, let's multiply equality (17) by 𝛽𝛥𝑡:
3
(𝑎̄ 𝛽𝛥𝑡(1 − 𝛼𝑓 ) + 𝛾𝐷)𝑥𝑛+1 + (1 − 𝛼𝑓 )𝛽𝛥𝑡𝑥𝑛+1 =
1 − 𝛼𝑚 (18)
= 𝛽𝛥𝑡𝐶𝑛+1 − 𝛾𝐷𝐵𝑛+1 + (𝛽𝛥𝑡𝐷 − 𝛽) 𝐴𝑛+1 .
𝛾
From here we find 𝑥. From the first equality (16), we find 𝑣:
𝛾 𝛾 𝛾(𝑥𝑛+1 + 𝐵𝑛+1 )
𝑣𝑛+1 − 𝑥𝑛+1 = 𝐵𝑛+1 − 𝐴𝑛+1 , 𝑣𝑛+1 = − 𝐴𝑛+1 .
𝛽Δ𝑡 𝛽Δ𝑡 𝛽Δ𝑡
From the first equality (14), we find 𝑎:
1 1 𝑣𝑛+1 + 𝐴𝑛+1
𝑎𝑛+1 − 𝑣𝑛+1 = 𝐴𝑛+1 , 𝑎𝑛+1 = .
𝛾Δ𝑡 𝛾Δ𝑡 𝛾Δ𝑡
Let's write down the obtained formulas of the generalized α method for solving a nonlinear oscillator
near the assembly catastrophe:
3
(𝑎̄ 𝛽𝛥𝑡(1 − 𝛼𝑓 ) + 𝛾𝐷)𝑥𝑛+1 + (1 − 𝛼𝑓 )𝛽𝛥𝑡𝑥𝑛+1 −
1 − 𝛼𝑚
−𝛽𝛥𝑡𝐶𝑛+1 + 𝛾𝐷𝐵𝑛+1 − (𝛽𝛥𝑡𝐷 − 𝛽) 𝐴𝑛+1 = 0,
𝛾 (19)
𝛾 1
𝑣𝑛+1 = (𝑥𝑛+1 + 𝐵𝑛+1 ) − 𝐴𝑛+1 , 𝑎𝑛+1 = (𝑣 + 𝐴𝑛+1 ).
𝛽𝛥𝑡 𝛾𝛥𝑡 𝑛+1
The solution (19) will be a curve, any changes of which depend on the parameter values 𝛿, 𝑎̄ and 𝑏.
Formulas (19) demonstrate formulas for solving a nonlinear oscillator near the assembly catastrophe by
the generalized 𝛼-method.
2.3. Bossak’s method
Bossak’s method is a continuation of Newmark’s method. The derivation of the formulas for solving
the nonlinear oscillator near the assembly catastrophe is carried out similarly to the Newmark’s method.
The difference is that the values of the second derivatives at 𝑛 + 1 time steps are taken as averages
(taking into account the average theorem).
The equality that we need to solve by Bossak’s method has the form (3). Everything that is on the
𝑛 + 1 step is unknown and it is necessary to find their value. To do this, we will transfer all the variables
𝑛 + 1 steps to the left, and n steps to the right:
𝛾𝛥𝑡 𝑎𝑛+1 − 𝑣𝑛+1 = −𝑣𝑛 − (1 − 𝛾)𝛥𝑡𝑎𝑛
1
𝛽𝛥𝑡 2 𝑎𝑛+1 − 𝑥𝑛+1 = −𝑥𝑛 − 𝛥𝑡𝑣𝑛 − ( − 𝛽) 𝛥𝑡 2 𝑎𝑛 (20)
2
3
{(1 − 𝛼)𝑎𝑛+1 + 𝛿𝑣𝑛+1 + 𝑎̄ 𝑥𝑛+1 + 𝑏 + 𝑥𝑛+1 = −𝛼𝑎𝑛
We see that the equalities in (20) are written in vector form. Therefore, to simplify the notations, we
introduce the vectors
1
𝐴𝑛+1 = −𝑣𝑛 − (1 − 𝛾)𝛥𝑡 𝑎𝑛 , 𝐵𝑛+1 = −𝑥𝑛 − 𝛥𝑡v𝑛 − ( − 𝛽) 𝛥𝑡 2 𝑎𝑛 , 𝐶𝑛+1 = −𝛼𝑚 𝑎𝑛 .
2
𝛾𝛥𝑡 𝑎𝑛+1 − 𝑣𝑛+1 = 𝐴𝑛+1
2𝑎
{ 𝛽𝛥𝑡 𝑛+1 − 𝑥𝑛+1 = 𝐵𝑛+1 (21)
3
(1 − 𝛼)𝑎𝑛+1 + 𝛿𝑣𝑛+1 + 𝑎̄ 𝑥𝑛+1 + 𝑏 + 𝑥𝑛+1 = 𝐶𝑛+1
Next, we divide the first equality (21) by 𝛾𝛥𝑡, and the second equality – by 𝛽𝛥𝑡 2.
1 1
𝑎𝑛+1 − 𝑣𝑛+1 = 𝐴
𝛾𝛥𝑡 𝛾𝛥𝑡 𝑛+1
1 1 (22)
𝑎𝑛+1 − 𝑥𝑛+1 = 𝐵
𝛽𝛥𝑡 2 𝛽𝛥𝑡 2 𝑛+1
3
{(1 − 𝛼)𝑎𝑛+1 + 𝛿𝑣𝑛+1 + 𝑎̄ 𝑥𝑛+1 + 𝑏 + 𝑥𝑛+1 = 𝐶𝑛+1
Now subtract the first from the second equality (22), and subtract the first multiplied by (1 − 𝛼)
from the third.
1 1 1 1
𝑣𝑛+1 − 𝑥𝑛+1 = 𝐵𝑛+1 − 𝐴 ,
𝛾𝛥𝑡 𝛽𝛥𝑡 2 𝛽𝛥𝑡 2 𝛾𝛥𝑡 𝑛+1
(23)
1−𝛼 3 1−𝛼
(𝛿 + ) 𝑣𝑛+1 + 𝑎̄ 𝑥𝑛+1 + 𝑏 + 𝑥𝑛+1 = 𝐶𝑛+1 − 𝐴 .
{ 𝛾𝛥𝑡 𝛾𝛥𝑡 𝑛+1
For simplification, we introduce 𝐷 = (𝛿 + (1 − 𝛼)⁄(𝛾𝛥𝑡)). Rewrite (23):
𝛾 𝛾
𝑣𝑛+1 − 𝑥𝑛+1 = 𝐵 − 𝐴𝑛+1
𝛽𝛥𝑡 𝛽𝛥𝑡 𝑛+1
1 (24)
3
𝐷𝑣𝑛+1 + 𝑎̄ 𝑥𝑛+1 + 𝑏 + 𝑥𝑛+1 = 𝐶𝑛+1 − 𝐴
{ 𝛾𝛥𝑡 𝑛+1
Next, we subtract from the second equality (24) the first multiplied by 𝐷.
3 𝛾 1−𝛼 𝛾
𝑎̄ 𝑥𝑛+1 + 𝑏 + 𝑥𝑛+1 + 𝐷𝑥𝑛+1 = 𝐶𝑛+1 − 𝐴𝑛+1 − 𝐷𝐵 + 𝐷𝐴𝑛+1 . (25)
𝛽𝛥𝑡 𝛾𝛥𝑡 𝛽𝛥𝑡 𝑛+1
In order to avoid division by zero, let's multiply equality (25) by 𝛽Δ𝑡:
3 )
(1 − 𝛼)𝛽
𝛽Δ𝑡(𝑎̄ 𝑥𝑛+1 + 𝑏 + 𝑥𝑛+1 + 𝛾𝐷𝑥𝑛+1 = 𝛽Δ𝑡𝐶𝑛+1 − 𝐴𝑛+1 − 𝛾𝐷𝐵𝑛+1 + 𝐷𝐴𝑛+1 . (26)
𝛾
From here we find 𝑥. From the first equality (24), we find 𝑣:
𝛾 𝛾 𝛾(𝑥𝑛+1 + 𝐵𝑛+1 )
𝑣𝑛+1 − 𝑥𝑛+1 = 𝐵𝑛+1 − 𝐴𝑛+1 , 𝑣𝑛+1 = − 𝐴𝑛+1 .
𝛽𝛥𝑡 𝛽𝛥𝑡 𝛽𝛥𝑡
We find 𝑎 from the first equality:
1 1 1
𝑎𝑛+1 − 𝑣𝑛+1 = 𝐴𝑛+1 , 𝑎𝑛+1 = (𝑣 + 𝐴𝑛+1 ).
𝛾Δ𝑡 𝛾Δ𝑡 𝛾Δ𝑡 𝑛+1
Let's write down the obtained formulas of Bossack’s method for solving a nonlinear oscillator near
the assembly catastrophe:
3 )
𝛽𝛥𝑡(𝑎̄ 𝑥𝑛+1 + 𝑏 + 𝑥𝑛+1 + 𝛾𝐷𝑥𝑛+1 − 𝛽𝛥𝑡𝐶𝑛+1 +
(1 − 𝛼)𝛽
+ 𝐴𝑛+1 + 𝛾𝐷𝐵𝑛+1 − 𝐷𝐴𝑛+1 = 0,
𝛾 (27)
𝛾 1
𝑣𝑛+1 = (𝑥𝑛+1 + 𝐵𝑛+1 ) − 𝐴𝑛+1 , 𝑎𝑛+1 = (𝑣 + 𝐴𝑛+1 ).
𝛽𝛥𝑡 𝛾𝛥𝑡 𝑛+1
Solution (27) is represented as a curve.
2.4. Programs and Calculation experiments
The general structure of the developed software consists of modules, the names of which are shown
in Figure 1. The software was developed in the computer mathematics environment MATLAB.
Module of Initial
Conditions
Program Software for
Module of Results Modelling Dynamic Module of Glogal
Visualization Systems Constants
Module of Solving
Ninlinear Equetions on one
time-step
Figure 1. Structure of Software of Modelling Dynamic Systems
The developed fragment of the program in MATLAB, with the help of which the solution is found
at each time step, has the following form:
t(i+1) = t(i) + h;
A_n = -v(i) - (1-gamma)*h*A(i);
B_n = -x(i) - h*v(i) - (0.5-betta)*h^2*A(i);
C_n = -al_m*A(i) - g*k2*v(i) - a*k2*x(i) - (k2*x(i))^3 - b;
D = betta*h^2*C_n - B_n*k1 - g*k3*h*(gamma*B_n - betta*h*A_n);
A = k3*betta*h^2;
B = 3*k2*k3^2*betta*h^2*x(i);
C = (3*k2^2*k3*x(i)^2+a*k3)*betta*h^2 + k1 + k3*g*gamma*h;
x(i+1) = Newton (x(i));
v(i+1) = gamma/(betta*h) * (x(i+1) + B_n)-A_n;
A(i+1) = 1/(gamma*h) * (v(i+1) + A_n);
i = i+1;
The obtained dependences of the methods for the considered dynamic systems were tested on
different values of the initial conditions and model constants. Some of the results are shown in Figure
2 and Figure 3. Initial Condition is 𝑥(0) = 1; 𝑣(0) = 2.
Figure 2. Result of numeric solution for equation (3) with parameters: left – 𝑎 = 𝑏 = 𝑔 = 1; right –
𝑎 = 3; 𝑏 = 20; 𝑔 = 4
Figure 3. Result of numeric solution for equation (3) with parameters 𝑎 = 0,5; 𝑏 = 1; 𝑔 = 0,2
Classical methods (Euler’s method and Runge-Kutta method) gave less accurate results.
3. General conclusions
In the work, the main dependencies are obtained on the basis of the implicit methods of Newmark,
Bossak and the generalized 𝛼-method for the example of a nonlinear oscillator near the assembly
catastrophe. Curves of changes are plotted for different values of control parameters. The programs of
implicit methods of Newmark, Bossack and the generalized α-method, Euler method and Runge-Kutta
method with automatic step selection for the solution of a nonlinear oscillator near the assembly
catastrophe are also implemented. Appropriate software for modeling deterministic nonlinear systems,
which are described by Cauchy problems of the second order for ordinary differential equations, has
been developed.
4. References
[1] Franke, M.E., Lynch, P.M., Healy, A.J., “A 50-Year History of the Dynamic Systems and Control
Division” ASME J. Dyn. Sys. Meas., and Control. Pp. 223–233.
[2] Craig A. Kluever. Dynamic Systems: Modeling, Simulation, and Control. ISBN: 978-1-119-
63540-6. Wiley; 1st edition. 2015, 496 p.
[3] Kuznetsov A.P. Fluctuations. Catastrophes. Bifurcations. Chaos. Saratov. 2000. 99 p.
[4] Kuznetsov A.P., Kuznetsov S.P., Ryskin N.M. Nonlinear oscillations. M.: Fizmatlit, 2002. 309 p.
[5] Lambert J. Computational Methods in Ordinary Differential Equations, 1973.
[6] Shampine L. Numerical solution of ordinary differential equations, Chapman & Hall, 1994.
[7] Shampine L., Gordon M. Computer Solution of Ordinary Differential Equations, Freeman, 1975.
[8] Shampine L., Gladwell I., Thompson S. Solving ODEs with MATLAB, Cambridge, 2003.