Open access peer-reviewed chapter

Mathematical Models of Oscillators with Memory

By Roman Ivanovich Parovik

Submitted: July 9th 2018Reviewed: October 4th 2018Published: November 12th 2018

DOI: 10.5772/intechopen.81858

Downloaded: 827


The chapter proposes a mathematical model for a wide class of hereditary oscillators, which is a Cauchy problem in the local formulation. As an initial model equation, an integrodifferential equation of Voltaire type was introduced, which was reduced by means of a special choice of difference kernels to a differential equation with nonlocal derivatives of fractional-order variables. An explicit finite-difference scheme is proposed, and questions of its stability and convergence are investigated. A computer study of the proposed numerical algorithm on various test examples of the hereditary oscillators Airy, Duffing, and others was carried out. Oscillograms and phase trajectories are plotted and constructed.


  • mathematical model
  • cauchy problem
  • heredity
  • derivative of fractional order
  • finite-difference scheme
  • stability
  • convergence
  • oscillograms
  • phase trajectory

1. Introduction

In the paper of the Italian mathematician Vito Volterra [1], the notion of heredity (memory), a property of a dynamical system characterized by nonlocality in time, is introduced, which consists in the dependence of its current state on a finite number of previous states. In another paper [2], Volterra investigated the hereditary oscillator—a vibration system with memory, which was written in the form of an integrodifferential equation with a difference kernel, a function of memory. Further, for such an oscillator, Volterra derived the law of total energy, in which an additional term appeared, responsible for the dissipation of energy in the vibrational system. This fact was confirmed in subsequent works.

In papers [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], fractal oscillators were considered, which represent the class of hereditary oscillators with a power-law function of memory. The peculiarity of such oscillators is that their mathematical description can be reduced to differential equations with nonlocal derivatives of fractional constant orders, which are investigated within the framework of the theory of fractional calculus [22, 23, 24].

In papers [6, 7, 8, 9, 11, 12, 13, 14, 18, 19, 20, 21], models of fractal linear oscillators were investigated in the sense of the Gerasimov-Caputo derivative and in papers [4, 10, 16]—in the sense of the Riemann-Liouville derivative. Analytical solutions of model equations in terms of a special function of Mittag-Leffler-type and generalized Wright-type function, oscillograms, and phase trajectories are constructed. It is shown that in the regime of free oscillations, the presence of memory effects in the system leads to attenuation of oscillations as a result of energy dissipation, and with allowance for external periodic action, it is possible to stabilize the amplitude of the oscillations, with the phase trajectories reaching the limit cycle and also the resonance effect.

In papers of the author [25, 26], the fractal parametric resonance (the fractal Mathieu oscillator) was investigated, and the Strutt-Ince diagrams of parametric resonance existence areas were constructed. It is shown that these regions strongly depend on the orders of the value fractional derivatives entering into the initial equation.

In a monograph by the Slovak mathematician Ivo Petras [10], the fractal nonlinear oscillator models whose differential equations contained fractional derivatives in the sense of Riemann-Liouville were considered and analyzed using numerical methods and considered the stability of the rest point of oscillatory systems. However, the stability and convergence of numerical methods have not been considered.

A further continuation of the investigation of hereditary oscillators is associated with the introduction of the derivatives of fractional variable orders in the model equations. This is due to the fact that the orders of fractional derivatives are related to the properties of the medium in which this or that process takes place and changes with time under the influence of external influence. Therefore, papers [27, 28, 29, 30] proposed that the models of fractal nonlinear oscillators (Duffing, Van der Pol, Van der Pol-Duffing, FitzHugh-Nagumo) were proposed and investigated using explicit finite-difference schemes, whose equations contain both the derivatives of the constants and variable fractional orders of the Gerasimov-Caputo and Riemann-Liouville types. With the help of computer experiments, the convergence of finite-difference schemes was shown, and estimates of the computational accuracy of the method were obtained; oscillograms and phase trajectories were constructed. However, the questions of stability and convergence were not formulated in the form of corresponding theorems.

In [31, 32], a new class of fractal oscillators was proposed and investigated; active fractal oscillators (AFOs)—nonlinear oscillators with external influences, which include the fractional Riemann-Liouville integral—were investigated. Such oscillators are constructed on the basis of the scheme of the radioelectronic аutogenerator with a fractional feedback circuit. Authors use the method of equivalent linearization to investigate AFOs and come to the conclusion that the self-oscillator is isochronous.

From the analysis of the above publications on the study of hereditary oscillator, we can conclude that the main tool for their study is numerical methods, for example, finite-difference schemes. In most cases, the authors leave without considering the questions of stability and convergence of finite-difference schemes and, even if they touch, then without formulating the corresponding theorems and proofs. Therefore, the goal of the present paper is to construct a finite-difference scheme for a wide class of hereditary (fractal) linear and nonlinear oscillators, prove its stability and convergence, formulate results in the form of corresponding theorems, and study finite-difference schemes on specific test examples.


2. Formulation of the problem

Consider the following model integrodifferential equation for the function xtC30T, where T>0:


where x¨t=d2x/dt2,ẋt=dx/dt, λ>0—given constant.

Eq. (1) describes a wide class of hereditary, depending on the form of the right-hand side (function fxtt) of linear or nonlinear oscillators.

Definition 1.Functions K1tηand K2tη—difference kernels in Eq. (1)—will be called memory functions, since they define the notion of heredity (memory), which was introduced in the work of the Italian mathematician Vito Volterra [2].

Definition 2.A nonlinear function fxtton the right-hand side of Eq. (1) satisfies a Lipschitz condition with respect to a variable xt:


L—Lipschitz constant.

Eq. (1) describes a broad class of hereditary nonlinear oscillators, depending on the form of the function fxtton its right-hand side, and the parameter λhas the meaning of the coefficient of friction.

Note that in Definition 1, the memory functions K1tηand K2tηcan be chosen arbitrarily, depending on the conditions of the particular problem. We will choose these functions power law, since power laws are often found in various fields of knowledge [33]. We choose the memory functions K1tηand K2tηin the form


where γt,βtC0T, Γt—Euler gamma function.

We give the following definitions.

Definition 3.Derivatives of fractional variables of orders βtand γtGerasimov-Caputo type: we call the following operators of fractional differentiation:


We note that in the case when the functions βtand γtin the relations (4) and (5) are constants, we arrive at the definitions of the fractional derivative in the sense of Gerasimov-Caputo [34, 35], and in the case when these constants β=2and γ=1, the operators of fractional differentiation (4) become classical derivatives of the second and first orders.

Taking into account Definition 3, the model Eq. (1) can be rewritten in a more compact form:


For Eq. (5), the initial conditions in the local formulation are valid:


where α0and α1—given constants. As a result, we arrive at the Cauchy problems (5) and (6), which we will investigate.


3. Explicit finite-difference scheme

We construct an explicit finite-difference scheme for the Cauchy problems (5) and (6). We divide the time interval 0Tinto Nequal parts with a step τ=T/N. We introduce the grid function xtk=xk, where tk=,k=1,,N1.

The derivatives of the variables of fractional orders in Eq. (5) are approximated according to the relations in [36, 37]:


then the formulas will refer to the formula (7)


Here, to shorten the record, βtk=βk,γtk=γk.

Taking into account relation (7), the Cauchy problems (5) and (6) in the difference formulation will have the form


Here, to shorten the record, fk=fxktk. We write the problem (8) explicitly:


We note that the weight coefficients aj,kand bj,khave properties, which we formulate in the form of the following lemmas.

Lemma 1.For anykweights, coefficientsaj,kandbj,k, as well as coefficientsAkandBk, have the following properties:

  1. j=0k1aj,k=k2βk,j=0k1bj,k=k1γk,

  2. 1=a0,k>a1,k>>0,1=b0,k>b1,k>>0,

  3. Ak0,Bk0.

Proof.The first property follows the definition of weight coefficients:



The second property is proven in the following way. We introduce two functions:


where x>0. These functions are decreasing. Really derived from these functions


Therefore, the second property holds. The third property follows also the properties of the gamma function. The lemma is proven.

Let ¯0tβtxηand ¯0tγtxη—approximations of differential operators of Gerasimov-Caputo types 0tβtxηand 0tγtxη. Then, we have the following lemma.

Lemma 2.Approximations ¯0tβtxηand ¯0tγtxηoperators of the Gerasimov-Caputo types 0tβtxηand 0tγtxηsatisfy the following estimates:


where С1and С2—constants that are independent of the parameter τ.

Proof.Using the first property of Lemma 1and Definition 3, we obtain


Similarly, we can obtain the second estimate from Eq. (10).


The lemma is proven.

Investigation.According to Lemma 2, it can be shown that the explicit finite-difference scheme (9) has an error ε=Oτ. This fact will be used in computer experiments in determining the computational accuracy of the numerical method.

Lemma 3.The sums in the finite-difference scheme (9) have the following representations:


Proof.The representation (11) follows the properties of the sum. Indeed, by opening the sums in Eq. (9) and grouping the terms properly, we arrive at the representation (11).

Using Lemma 3, the finite-difference scheme (9) can be rewritten as


Or in matrix form


where the matrix M=mij, i=1,,N1,j=1,,N1:


Theorem 1.An explicit finite-difference scheme (9) converges with the first orderx¯kxk=Oτif the following condition is satisfied:


Proof.Let X¯k=x¯0x¯N2Tbe the exact solution of system (8) and the error vector ek+1=X¯k+1Xk+1,e0=0.. Then, system (8), with allowance for Lemma 2, can be written as follows:




We note that for any k, the inequality holds LkL. Consider the norm for the matrix M:M=maxij=1k1mij, we obtain


According to Lemma 1, we note that the inequality holds 2A1A1+B10. Suppose that the condition is satisfied Ai1Bi1, then the second diagonal element satisfies the inequality 12A1A1+B12in the matrix M, and the remaining diagonal elements are equal to the inequality 0Ai12ai2,i1Bi1bi2,i1Ai1+Bi11; these elements with iN1, in view of properties 2 and 3, tend to be zero.

The remaining elements of the matrix (16) also possess these properties. We also note that the matrix Mis a matrix with a diagonal predominance for small values λ.

Therefore, the sum of the elements of the second row in the matrix Msatisfies the condition 12A1A1+B1+A1a0,1B1b0,1A1+B13. Further, by virtue of properties 2 and 3 of Lemma 1, it is obvious that the sum of the remaining terms also satisfies these conditions. Therefore, the following estimate is valid for the norm: 1M3.

Note that for the values of the parameter λ1the norm M1, however, the condition number μM1is violated and the diagonal transformation is violated; therefore, it is necessary to decrease the step τ.

Further, for any constant С>0independent of τ, and the error rate, the following estimate holds:


We introduce the notation in Eq. (17): sk=3+LAk+Bk,s=. Then, we obtain the following estimate:


Substituting into Eq. (18) r=k1, we obtain


From the second initial condition (6) it follows: e1e0and С0=p=1ksp.

Now according to our assumption Ai1Bi1, which leads us to the relation


Condition (19) begins to work at such values λ, when many of conditionalities arise μМ1, and for sufficiently small values λ, it suffices that the step satisfies the inequality τ1. Therefore, we arrive at the relation (14). The theorem is proven.

We note that in [38] the authors used the classical Lax theorem, which holds for local finite-difference schemes, to prove the convergence of the scheme. For nonlocal finite-difference schemes, the convergence must be proven independently.

We consider the stability of an explicit finite-difference scheme (4). Suppose that Xkand two Ykare different solutions of the matrix Eq. (12) with initial conditions X0and Y0.

Theorem 2.An explicit finite-difference scheme (9) is conditionally stable ifcondition (14) is satisfied and the estimate holdsYkXkCY0X0for anyk, whereС>0does not depend on the stepτ.

Proof.We introduce the notation: ek+1=Yk+1Xk+1.Then, Eq. (12) can be written in the form ek+1=Mek+Fe,k. Here, as it was said in.


According to Theorem 1, we have the following estimate:


Therefore, the following estimate holds


With r=k1, we obtain ek+1С0e1C0e0and С0=p=1ksp.

The last inequality follows the second condition of problem (6). Therefore, if X0there is a perturbation, then it does not lead to a large increase in the error of the numerical solution. However, for large values λ, many of conditionalities μМ1arise, and therefore it is necessary to decrease the step τ; according to Eq. (19), for small values λ, the estimate is valid τ1. Then, the system is stable if condition (14) is satisfied. The theorem is proven.


4. Results of modeling

Consider the work of the explicit finite-difference scheme (9) on specific examples. We show that the scheme (9) has the first order of accuracy. Since in the general case, the exact solution of the Cauchy problems (5) and (6) cannot be written in analytical form, we will use the double conversion method. For this, we introduce two parameters: ξ=maxixix2i—absolute error between the numerical solution xiin step τand the numerical solution x2iin step τ/2. Then, the order of computational accuracy pcan be estimated by the formula


We note that in the case when the fractional parameters in the scheme (9) do not change and have the following values of βk=2and γk=1, we arrive at the classical local explicit finite-difference scheme with the second order of accuracy.

The numerical algorithm (9) was implemented in Maple software.

Example 1.Suppose that the right-hand side in Eq. (1) has the form


Then, Eq. (5) describes a linear hereditary Airy oscillator, which was considered in the author’s papers [21, 39] and has the following form.


We choose the initial condition (6) for simplicity by homogeneous:


We note that the Airy oscillator is used in optics in the simulation of Airy laser beams.

In this case, the explicit finite-difference scheme (9) has a more specific form:


For the explicit finite-difference scheme (20), we choose the following values of the control parameters: T=1. λ=1, δ=5, φ=10, ω=10and βt=1.80.03sinωt,γt=0.80.05cosωt. And during the simulation, we will change the number of nodes Nin the calculation grid.

Note that the values of the selected parameters for Example 1satisfy the conditions of Theorems 1and 2, which is indirectly confirmed by the results of modeling for different values Nof the nodes of the computational grid (Table 1).


Table 1.

Results of numerical simulation.

From Table 1 we can notice that when the number of calculated nodes in the grid doubles in nodes N, the maximum error in absolute value decreases twice, and the order of computational accuracy ptends to unite.

This confirms that the explicit finite-difference scheme (9) and in particular the scheme (20) for Example 1have the first order of accuracy, and since condition (14) is satisfied, then convergence with the same order.

In Figure 1 the oscillogram (Figure 1a) and the phase trajectory (Figure 1b) are shown for Example 1at the parameter value T=10,N=1000.It can be noted that with time the amplitude of the oscillations is established and as a result the phase trajectory reaches the limit cycle. Another situation arises in the case of free oscillations δ=0(Figure 2).

Figure 1.

The oscillogram (a) and the phase trajectory (b) for Example 1 with the parameter valuesT=10,N=1000.

Figure 2.

Oscillogram (a) and phase trajectory (b) for Example 1 with initial conditionsx0=0.1,ẋ0=0.2, andδ=0.

The amplitude of the oscillations decays (Figure 2a), and the phase trajectory twists into a spiral (Figure 2b). The dissipation of energy in this case occurs as a result of the presence of friction with a coefficient λand also the “memory” effect, which gives an additional term in the ratio for the total energy of the oscillatory system (Figure 3).

Figure 3.

The oscillogram (a) and the phase trajectory (b) for Example 1 with initial conditionsx0=0.1,ẋ0=0.2, andδ=0,λ=0.

This fact is confirmed by the results of [2]. Consider the following example of a nonlinear hereditary oscillator.

Example 2.Let that in Eq. (1) the right-hand side has the form.


and we choose the initial conditions (6) to be homogeneous:


In this case, Eq. (5) describes the Duffing fractional oscillator [18]:


The explicit finite-difference scheme (9) for this case has the form


For the explicit finite-difference scheme (21), we take the values of the control parameters as follows: T=1, λ=0.3, δ=2, and φ=ω=1.

Remark.Note that this choice of control parameter values is ensured by the condition (14) for Theorems 1and 2. The results of numerical simulation for Example 2are given in Table 2.


Table 2.

Results of numerical simulation.

Note from Table 2 that for Example 2, with an increase in the number of design nodes N, the maximum error ξin absolute value decreases and the order of computational accuracy ptends to unite. This indicates that the explicit finite-difference scheme (21) has the first order of accuracy.

Let’s perform numerical simulation according to the scheme (21) with the values of the following parameters, T=100, N=2000, and δ=50, and leave the remaining parameters unchanged. Let us construct an oscillogram and a phase trajectory (Figure 3).

The oscillogram (Figure 4a) has a constant amplitude of a more complex shape at its minima and maxima, which is reflected in the phase trajectory (Figure 4b). The phase trajectory enters a complex two-loop limit cycle. The presence of such loops, apparently, is associated with the effects of memory in the oscillatory system [40].

Figure 4.

The oscillogram (a) and the phase trajectory (b) forExample 2.

Figure 5 shows the case of free oscillations for Example 2. It is seen that the presence of friction and memory effects in the oscillatory system intensify energy dissipation, which leads to damping of the oscillations (Figure 5a) and a phase trajectory—a twisting spiral (Figure 5b). Indeed, if there is no friction λ=0in the oscillatory system, we obtain an oscillogram and a phase trajectory as in Figure 6.

Figure 5.

Oscillogram (a) and phase trajectory (b) forExample 3with initial conditionsx0=0.1,ẋ0=0.2, andδ=0.

Figure 6.

Oscillogram (a) and phase trajectory (b) forExample 2with initial conditionsx0=0.1,ẋ0=0.2, andδ=0,λ=0.

Example 3.Suppose that in Eq. (1) the right-hand side has the form


where bis the spring travel speed; cis the surface adhesion energy; ωis the frequency of free oscillations; and an=2n01cosπnτcosh2πτis the coefficients of the expansion of the Fourier series.

Eq. (5) with the right-hand side of Eq. (22) describes the hereditary stick–slip effect [20]. The stick–slip effect is encountered in tribology problems, for example, when the movement of a load on a spring along a surface is investigated. Due to adhesion, the load adheres to the surface, and due to the tension of the spring, it breaks and slides along it, and its oscillations occur [41, 42]. The stick–slip effect can also be incorporated into the mechanical model of an earthquake in the subduction zone of lithospheres plates [43].

In [43] it was said that in order to obtain a reliable solution it suffices to take the first seven coefficients anin the expansion of the function (22). The values of these coefficients are taken from [43] a1=0.436, a2=0.344, a3=0.164, a4=0.058, a5=0.021, a6=0.004, and a7=0.003. Values of control parameters are βt=1.80.03sinπtγt=0.60.04cosπt, N=3000δ=50, τ=0.05, λ=0.3, b=1, ω=1, and x0=0,ẋ0=0.3.

Figure 7 shows the calculated displacement curves, displacement velocities, and phase trajectory. Figure 7a shows the oscillogram for Example 3. It can be seen that during the separation, the load experiences oscillations and the rate of such oscillations in the potential well attenuates rather slowly (Figure 7b). This effect is the eradication of the process. The phase trajectory in Figure 7c shows that the potential wells are stable focuses.

Figure 7.

Calculated curves obtained fromformula (6)from [20] (curve 1) andformula (9)(curve 2): (a) oscillogram, (b) oscillator speed, and (c) phase trajectory.


5. Conclusion

A mathematical model characterizing a wide class of hereditary oscillators is proposed and studied. The model is a differential Cauchy problem with derivatives of fractional-order variables of the Gerasimov-Caputo types (5) and (6). Using the theory of finite-difference schemes, a nonlocal explicit finite-difference scheme (9) was constructed with the first order of accuracy. Questions of its stability and convergence, which are formulated in the form of corresponding theorems, were studied.

The main result of the paper can be formulated as follows: an explicit finite-difference scheme is conditionally stable and converges if criterion (14) is satisfied. With the help of computational examples, it was shown that the scheme (9) has the first order of accuracy. It is confirmed that in the case of free oscillations, the presence of friction and heredity increases dissipation of energy, which leads to attenuation of oscillations.

One of the continuations of the investigation of the Cauchy problems (5) and (6) is a generalization of it:


Another continuation of the research is related to the introduction of other memory functions K1tτ,K2tτinto the model Eq. (1), which leads to different model equations with different derivatives of fractional orders, and also the Cauchy problems (5) and (6) can be written in terms of the local fractional derivative [44–46].

The question of the stability of the rest points of dynamical systems described by the Cauchy problems (5) and (6) is also interesting, by analogy with the papers [47, 48].



The work was carried out according to the state task within the framework of research work Vitus Bering Kamchatka State University on the topic “Application of fractional calculus in the theory of oscillatory processes” No.AAAA-A17-117031050058-9 and with the support of the grant of the President of the Russian Federation MK-1152.2018.1.

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Roman Ivanovich Parovik (November 12th 2018). Mathematical Models of Oscillators with Memory, Oscillators - Recent Developments, Patrice Salzenstein, IntechOpen, DOI: 10.5772/intechopen.81858. Available from:

chapter statistics

827total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Quantum Harmonic Oscillator

By Coşkun Deniz

Related Book

First chapter

Application of IR Thermography for Studying Deformation and Fracture of Paper

By Tatsuo Yamauchi

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us