From Mathematics to Oscillators

Chapter 1

Abstract

trajectory

3

1. Introduction

Mathematical Models of

Oscillators with Memory

The chapter proposes a mathematical model for a wide class of hereditary oscil-

Keywords: mathematical model, cauchy problem, heredity, derivative of fractional

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 vibra-

In papers [3–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 inves-

In papers [6–9, 11–14, 18–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 Wrighttype 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

tigated within the framework of the theory of fractional calculus [22–24].

order, finite-difference scheme, stability, convergence, oscillograms, phase

lators, 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 finitedifference 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.

tional system. This fact was confirmed in subsequent works.

Roman Ivanovich Parovik

#### Chapter 1

## Mathematical Models of Oscillators with Memory

Roman Ivanovich Parovik

#### Abstract

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 finitedifference 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.

Keywords: 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–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–24].

In papers [6–9, 11–14, 18–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 Wrighttype 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–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 selfoscillator 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 x tð Þ <sup>∈</sup> <sup>C</sup><sup>3</sup> ð0; TÞ, where T > 0:

Mathematical Models of Oscillators with Memory DOI: http://dx.doi.org/10.5772/intechopen.81858

$$\int\_{0}^{t} K\_{1}(t - \eta) \ddot{\boldsymbol{x}}^{\cdot}(\eta) d\eta + \lambda \int\_{0}^{t} K\_{2}(t - \eta) \dot{\boldsymbol{x}}^{\cdot}(\eta) d\eta = f(\boldsymbol{x}(t), t), \tag{1}$$

where x t €ðÞ¼ <sup>d</sup> <sup>=</sup>dt x t <sup>=</sup>dt, <sup>λ</sup> > 0—given constant. <sup>2</sup> x <sup>2</sup> , \_ðÞ¼ dx

Eq. (1) describes a wide class of hereditary, depending on the form of the righthand side (function fðx tð Þ; tÞ) of linear or nonlinear oscillators.

Definition 1. Functions K1ðt � ηÞ and K2ðt � ηÞ—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 fð Þ x tð Þ; t on the right-hand side of Eq. (1) satisfies a Lipschitz condition with respect to a variable x tð Þ:

$$|f(\mathbf{x}\_1(t), t) - f(\mathbf{x}\_2(t), t)| \le L|\mathbf{x}\_1(t) - \mathbf{x}\_2(t)|,\tag{2}$$

L—Lipschitz constant.

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

Note that in Definition 1, the memory functions K1ðt � ηÞ and K2ðt � ηÞ 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 K1ðt � ηÞ and K2ðt � ηÞ in the form

$$K\_1(t - \eta) = \frac{(t - \eta)^{1 - \beta(t)}}{\Gamma(2 - \beta(t))},\\ K\_2(t - \eta) = \frac{(t - \eta)^{-\gamma(t)}}{\Gamma(2 - \gamma(t))},\\ \mathbf{1} < \boldsymbol{\beta}(t) < \mathbf{2}, \; \mathbf{0} < \gamma(t) < \mathbf{1},\tag{3}$$

where γð Þt , βð Þt ∈C½0; T�, Γð Þt —Euler gamma function.

We give the following definitions.

Definition 3. Derivatives of fractional variables of orders βð Þt and γð Þt Gerasimov-Caputo type: we call the following operators of fractional differentiation:

$$\partial\_{0t}^{\beta(t)}\mathbf{x}(\eta) = \frac{\mathbf{1}}{\Gamma(2-\beta(t))} \Big|\_{0}^{t} \frac{\ddot{\mathbf{x}}(\eta)d\eta}{(t-\eta)^{\beta(t)-1}},\\\partial\_{0t}^{\gamma(t)}\mathbf{x}(\eta) = \frac{\mathbf{1}}{\Gamma(1-\gamma(t))} \Big|\_{0}^{t} \frac{\dot{\mathbf{x}}(\eta)d\eta}{(t-\eta)^{\gamma(t)}}.\tag{4}$$

We note that in the case when the functions βð Þt and γð Þt in 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 β ¼ 2 and γ ¼ 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:

$$
\partial\_{0t}^{\beta(t)} \varkappa(\eta) + \lambda \partial\_{0t}^{\chi(t)} \varkappa(\eta) = f(\varkappa(t), t). \tag{5}
$$

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

$$\mathfrak{x}(\mathbf{0}) = a\_{\mathbf{0}}, \dot{\mathfrak{x}}(\mathbf{0}) = a\_{\mathbf{0}} \tag{6}$$

where α<sup>0</sup> and α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 ½0; T] into N equal parts with a step τ ¼ T=N. We introduce the grid function x tð <sup>k</sup>Þ ¼ xk, where tk ¼ kτ, k ¼ 1, …, N - 1.

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

$$\partial\_{0t}^{\beta(t)}\mathcal{X}(\eta) \approx \mathcal{A}\_k \sum\_{j=0}^{k-1} a\_{j,k} (\varkappa\_{k-j+1} - 2\varkappa\_{k-j} + \varkappa\_{k-j-1}),\tag{7}$$

$$\partial\_{0t}^{\gamma(t)}\mathcal{X}(\eta) \approx \mathcal{B}\_k \sum\_{j=0}^{k-1} b\_{j,k} (\varkappa\_{k-j+1} - \varkappa\_{k-j-1}),$$

then the formulas will refer to the formula (7)

$$a\_{j,k} = \left(j+1\right)^{2-\beta\_k} - j^{2-\beta\_k}, \\ b\_{j,k} = \left(j+1\right)^{1-\gamma\_k} - j^{1-\gamma\_k},$$

$$A\_k = \frac{\pi^{-\beta\_k}}{\Gamma(3-\beta\_k)}, \\ B\_k = \frac{\lambda \pi^{-\gamma\_k}}{2\Gamma(2-\gamma\_k)},$$

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

$$A\_k \sum\_{j=0}^{k-1} a\_{j,k} \left(\mathbf{x}\_{k-j+1} - \mathbf{2}\mathbf{x}\_{k-j} + \mathbf{x}\_{k-j-1}\right) + B\_k \sum\_{j=0}^{k-1} b\_{j,k} \left(\mathbf{x}\_{k-j+1} - \mathbf{x}\_{k-j-1}\right) = f\_k,\tag{8}$$
 
$$\mathbf{x}\_0 = a\_0,\\ \mathbf{x}\_1 = a\_1 + \tau a\_0,$$

Here, to shorten the record, f <sup>k</sup> ¼ f xð <sup>k</sup>; tkÞ. We write the problem (8) explicitly:

$$\begin{split} \mathbf{x}\_{k+1} &= \frac{1}{A\_k + B\_k} \quad 2A\_k \mathbf{x}\_k - (A\_k - B\_k)\mathbf{x}\_{k-1} - A\_k \sum\_{j=1}^{k-1} a\_{j,k} \left( \mathbf{x}\_{k-j+1} - 2\mathbf{x}\_{k-j} + \mathbf{x}\_{k-j} \right) \Bigg( \frac{1}{A\_k} \\ &- \frac{B\_k}{A\_k + B\_k} \sum\_{j=1}^{k-1} b\_{j,k} \left( \mathbf{x}\_{k-j+1} - \mathbf{x}\_{k-j-1} \right) + f\_k. \end{split} \tag{9}$$

We note that the weight coefficients aj,k and bj,k have properties, which we formulate in the form of the following lemmas.

Lemma 1. For any k weights, coefficients aj,k and bj,k, as well as coefficients Ak and Bk, have the following properties:

<sup>1</sup> <sup>β</sup><sup>k</sup> <sup>1</sup> <sup>γ</sup><sup>k</sup> 1. ∑<sup>k</sup> j¼ - <sup>0</sup> aj,k <sup>¼</sup> <sup>k</sup><sup>2</sup>- , <sup>∑</sup><sup>k</sup> j¼ - <sup>0</sup> bj,k <sup>¼</sup> <sup>k</sup><sup>1</sup>- , 2. 1 ¼ a0,k > a1,k > … > 0, 1 ¼ b0,k > b1,k > … > 0, 3. Ak ≥ 0, Bk ≥ 0:

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

$$\begin{array}{c} \sum\_{j=0}^{k-1} a\_{j,k} = \sum\_{j=0}^{k-1} \left[ (j+1)^{2-\beta\_k} - j^{2-\beta\_k} \right] \begin{pmatrix} \mathbf{1} - \mathbf{0} + \mathbf{2}^{2-\beta\_k} - \mathbf{1} + \mathbf{3}^{2-\beta\_k} - \mathbf{2}^{2-\beta\_k} + \dots \\ + (k-1)^{2-\beta\_k} + k^{2-\beta\_k} - (k+1)^{2-\beta\_k} = k^{2-\beta\_k} \end{array}$$

$$\begin{split} \sum\_{j=0}^{k-1} b\_{j,k} &= \sum\_{j=0}^{k-1} \left[ \left( j + \mathbf{1} \right)^{1-\gamma\_k} - j^{1-\gamma\_k} \right] = \mathbf{1} - \mathbf{0} + \mathbf{2}^{1-\gamma\_k} - \mathbf{1} + \mathbf{3}^{1-\gamma\_k} - \mathbf{2}^{1-\gamma\_k} + \dots \\ & \dots + (k-1) \overset{1-\gamma\_k}{\cdot} + k^{1-\gamma\_k} - (k-1) \overset{1-\gamma\_k}{\cdot} = k^{1-\gamma\_k} \cdot \underbrace{\mathbf{1}}\_{\text{even. }\dots\text{even. }\dots\text{odd.}} \end{split}$$

$$\rho(\mathfrak{x}) = (\mathfrak{x} + \mathfrak{1})^{2 - \beta\_k} - \mathfrak{x}^{2 - \beta\_k} \text{ and } \eta(\mathfrak{x}) = (\mathfrak{x} + \mathfrak{1})^{1 - \gamma\_k} - \mathfrak{x}^{1 - \gamma\_k},$$

$$\rho'(\mathbf{x}) = (2 - \beta\_k) \left[ (\mathbf{x} + \mathbf{1})^{1 - \beta\_k} - \mathbf{x}^{1 - \beta\_k} \right] < 0,\\ \eta'(\mathbf{x}) = (\mathbf{1} - \gamma\_k) \left[ (\mathbf{x} + \mathbf{1})^{1 - \gamma\_k} - \mathbf{x}^{1 - \gamma\_k} \right] < \mathbf{0}.$$

$$\left|\partial\_{0t}^{\beta(t)}\mathcal{x}(\eta) - \overline{\partial}\_{0t}^{\beta(t)}\mathcal{x}(\eta)\right| \leq \mathcal{C}\_{1}\tau,\\ \left|\partial\_{0t}^{\gamma(t)}\mathcal{x}(\eta) - \overline{\partial}\_{0t}^{\gamma(t)}\mathcal{x}(\eta)\right| \leq \mathcal{C}\_{2}\tau,\tag{10}$$

\_ \_ <sup>τ</sup>1-γ<sup>k</sup> <sup>k</sup>1-γ<sup>k</sup> <sup>τ</sup>1-γ<sup>k</sup> <sup>τ</sup>1-<sup>γ</sup> <sup>k</sup>-<sup>1</sup> <sup>k</sup> <sup>k</sup>-<sup>1</sup> <sup>γ</sup>ð Þ<sup>t</sup> <sup>∑</sup> bj,k x t \_<sup>ð</sup> - <sup>j</sup>τÞ þ <sup>O</sup>ð Þ<sup>τ</sup> <sup>∑</sup> bj,k <sup>∂</sup>0<sup>t</sup> <sup>x</sup>ðηÞ ¼ <sup>½</sup> ] ¼ x t \_<sup>ð</sup> - <sup>j</sup>τÞ þ <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup><sup>k</sup> <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup><sup>k</sup> <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup><sup>k</sup> <sup>Þ</sup> <sup>Þ</sup> <sup>Þ</sup> <sup>j</sup>¼<sup>0</sup> <sup>j</sup>¼<sup>0</sup> τ<sup>1</sup>-γ<sup>k</sup> t <sup>1</sup>-γ<sup>k</sup> <sup>τ</sup><sup>1</sup>-<sup>γ</sup> <sup>k</sup>-<sup>1</sup> <sup>k</sup> <sup>k</sup>-<sup>1</sup> Þ þ <sup>O</sup> <sup>τ</sup><sup>2</sup> <sup>O</sup>ðτÞ ¼ <sup>∑</sup> bj,k x t \_<sup>ð</sup> - <sup>j</sup><sup>τ</sup> <sup>¼</sup> <sup>∑</sup> bj, <sup>k</sup> x t \_<sup>ð</sup> - <sup>j</sup>τÞ þ <sup>O</sup>ð Þ<sup>τ</sup> : <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup><sup>k</sup> <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup>k<sup>Þ</sup> <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup>k<sup>Þ</sup> <sup>Þ</sup> <sup>j</sup>¼<sup>0</sup> <sup>j</sup>¼<sup>0</sup> ðjþ ð <sup>1</sup>Þ<sup>τ</sup> <sup>k</sup>-<sup>1</sup> <sup>∑</sup> 1 <sup>∑</sup> <sup>γ</sup>ð Þ<sup>t</sup> <sup>1</sup> <sup>1</sup> <sup>ξ</sup>-<sup>γ</sup> <sup>x</sup>ðηÞ ¼ kx t \_ð - <sup>ξ</sup>Þd<sup>ξ</sup> <sup>¼</sup> bj,k <sup>∂</sup> x t \_ - <sup>η</sup><sup>j</sup> , <sup>0</sup><sup>t</sup> <sup>Γ</sup>ð<sup>1</sup> - <sup>γ</sup>k<sup>Þ</sup> <sup>Γ</sup>ð<sup>1</sup> - <sup>γ</sup>k<sup>Þ</sup> kj¼0 j¼0 jτ η<sup>j</sup> ∈½ jτ;ð j þ 1Þτ]: <sup>τ</sup>1-γ<sup>k</sup> <sup>k</sup>-<sup>1</sup> h i <sup>γ</sup>ð Þ<sup>t</sup> <sup>γ</sup>ð Þ<sup>t</sup> <sup>∂</sup> <sup>0</sup><sup>t</sup> <sup>x</sup>ðηÞ - <sup>∂</sup>0<sup>t</sup> <sup>x</sup>ð Þ<sup>η</sup> <sup>¼</sup> <sup>∑</sup> bj,k x tð - <sup>j</sup>τÞ - x t - <sup>η</sup><sup>j</sup> <sup>þ</sup> <sup>O</sup>ð Þ<sup>τ</sup> <sup>¼</sup> <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup>k<sup>Þ</sup> <sup>j</sup>¼<sup>0</sup> <sup>k</sup>-<sup>1</sup> <sup>τ</sup><sup>1</sup>-γ<sup>k</sup> <sup>k</sup><sup>1</sup>-γ<sup>k</sup> <sup>∑</sup> <sup>τ</sup><sup>1</sup>-γ<sup>k</sup> <sup>¼</sup> bj,k ·OðτÞ þ <sup>O</sup>ð Þ<sup>τ</sup> <sup>¼</sup> <sup>O</sup>ðτÞ þ <sup>O</sup>ð Þ<sup>τ</sup> <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup>k<sup>Þ</sup> <sup>j</sup>¼<sup>0</sup> <sup>Γ</sup>ð<sup>2</sup> - <sup>γ</sup><sup>k</sup> <sup>Þ</sup> ¼ OðτÞ þ OðτÞ ¼ Oð Þτ :

The lemma is proven.

Investigation. According to Lemma 2, it can be shown that the explicit finitedifference 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. Thesums in thefinite-differencescheme (9) havethefollowing representations:

$$\begin{split} & \sum\_{j=1}^{k-1} a\_{j,k} (\mathbf{x}\_{k-j+1} - 2\mathbf{x}\_{k-j} + \mathbf{x}\_{k-j-1}) = a\_{1,k} \mathbf{x}\_k + a\_{k-1,k} \mathbf{x}\_0 \\ & \quad + (a\_{k-2,k} - 2a\_{k-1,k}) \mathbf{x}\_1 + (a\_{2,k} - 2a\_{1,k}) \mathbf{x}\_{k-1} + \sum\_{j=2}^{k-2} (a\_{j+1,k} - 2a\_{j,k} + a\_{j-1,k}) \mathbf{x}\_{k-j}, \\ & \quad - \sum\_{j=1}^{k-1} b\_{j,k} (\mathbf{x}\_{k-j+1} - \mathbf{x}\_{k-j-1}) = b\_{1,k} \mathbf{x}\_k - b\_{k-1,k} \mathbf{x}\_0 + b\_{2,k} \mathbf{x}\_{k-1} - b\_{k-2,k} \mathbf{x}\_1 \\ & \quad + \sum\_{j=2}^{k-2} (b\_{j+1,k} - b\_{j-1,k}) \mathbf{x}\_{k-j}. \end{split} \tag{11}$$

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

$$\begin{aligned} \mathbf{x}\_1 &= a\_0 + \tau a\_1, \\ \mathbf{x}\_2 &= \frac{2A\_1}{A\_1 + B\_1} \mathbf{x}\_1 - \frac{A\_1 - B\_1}{A\_1 + B\_1} \mathbf{x}\_0 + \frac{f\_1}{A\_1 + B\_1}, k = 1, \\ \mathbf{x}\_{k+1} &= \frac{1}{A\_k + B\_k} \\ \begin{bmatrix} (A\_k(2 - a\_{1,k}) - B\_k b\_{1,k}) \mathbf{x}\_k - (A\_k(a\_{2,k} - 2a\_{1,k} + a\_{0,k}) + B\_k(b\_{2,k} - b\_{0,k})) \mathbf{x}\_{k-1} + f\_k \end{bmatrix} - \mathbf{x}\_1 \\ - \frac{(A\_k(a\_{k-2,k} - 2a\_{k-1,k}) - B\_k b\_{k-2,k}) \mathbf{x}\_1 - (A\_k a\_{k-1,k} - B\_k b\_{k-1,k}) \mathbf{x}\_0}{A\_k + B\_k} \\ - \frac{1}{A\_k + B\_k} \sum\_{j=2}^{k-2} \left[ A\_k \left( a\_{j+1,k} - 2a\_{j,k} + a\_{j-1,k} \right) + B\_k \left( b\_{j+1,k} - b\_{j-1,k} \right) \right] \mathbf{x}\_{k-j}, k = 2, \dots, n-1, \end{aligned}$$

Or in matrix form

$$X\_{k+1} = M X\_k + F\_k,\tag{12}$$

$$\tau \le \tau\_0 = \min \quad \mathbf{1}, \left( \frac{\mathfrak{A}\Gamma(2-\chi\_{i-1})}{\mathfrak{A}\Gamma(3-\beta\_{i-1})} \right)^{\frac{1}{\beta\_{i-1}-\chi\_i}} \binom{\chi}{i} \mathbf{j} = \mathbf{2}, \dots, N-1. \tag{14}$$

$$e\_{k+1} = Me\_k + F\_{\varepsilon,k} + O(\tau),\tag{15}$$

$$\begin{aligned} F\_{\epsilon,k} &= \frac{1}{A\_k + B\_k} \left( |f(\mathbf{x}\_1, t\_k) - f(\overline{\mathbf{x}}\_1, t\_k)|, \dots, |f(\mathbf{x}\_k, t\_k) - f(\overline{\mathbf{x}}\_k, t\_k)| \right)^T \\ &\le \frac{1}{A\_k + B\_k} (L\_1 e\_1, \dots, L\_k e\_k) = \Delta F\_k e\_k, \Delta F\_k = \frac{1}{A\_k + B\_k} \text{diag} \left( L\_1, \dots, L\_k \right)^T. \end{aligned}$$

According to Lemma <sup>1</sup>, we note that the inequality holds <sup>2</sup>A<sup>1</sup> <sup>≥</sup> 0. Suppose that <sup>A</sup>1þB<sup>1</sup> the condition is satisfied Ai-<sup>1</sup> ≥ Bi-1, then the second diagonal element satisfies the inequality <sup>1</sup> <sup>≤</sup> <sup>2</sup>A<sup>1</sup> <sup>≤</sup> <sup>2</sup> in the matrix <sup>M</sup>, and the remaining diagonal elements are <sup>A</sup>1þB<sup>1</sup> Ai-<sup>1</sup>ð2-ai-<sup>2</sup>,i-<sup>1</sup>Þ-Bi-<sup>1</sup>bi-<sup>2</sup>,i-<sup>1</sup> equal to the inequality <sup>0</sup> <sup>≤</sup> <sup>≤</sup> 1; these elements with Ai-1þBi-<sup>1</sup> i ! N - 1, 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 M is a matrix with a diagonal predominance for small values λ.

Therefore, the sum of the elements of the second row in the matrix M satisfies <sup>A</sup>1a0, <sup>1</sup>-B1b0, <sup>1</sup> the condition <sup>1</sup> <sup>≤</sup> <sup>2</sup>A<sup>1</sup> <sup>þ</sup> <sup>≤</sup> 3. Further, by virtue of properties <sup>2</sup> and <sup>3</sup> of <sup>A</sup>1þB<sup>1</sup> <sup>A</sup>1þB<sup>1</sup> 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: 1 ≤ k k M <sup>∞</sup> ≤ 3. Note that for the values of the parameter λ ≫ 1 the norm k k M <sup>∞</sup> ! 1, however, the condition number μð Þ M ≫ 1 is violated and the diagonal transformation is violated; therefore, it is necessary to decrease the step τ.

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

$$\left\|\left|e\_{k+1}\right|\right\|\_{\infty} \le \left\|\Delta F\_k + M\right\|\_{\infty} \left\|e\_k\right\|\_{\infty} + C\tau \le \left(3 + \frac{L}{A\_k + B\_k}\right) \left\|e\_k\right\|\_{\infty} + C\tau. \tag{17}$$

 We introduce the notation in Eq. (17): sk <sup>¼</sup> <sup>3</sup> <sup>þ</sup> <sup>L</sup> , s <sup>¼</sup> <sup>C</sup>τ. Then, we obtain AkþBk the following estimate:

$$\|e\_{k+1}\|\_{\infty} \le s\_k \|e\_k\|\_{\infty} + s \le s\_k \left(s\_{k-1} \|e\_{k-1}\|\_{\infty} + s\right) + s = s\_k s\_{k-1} \|e\_{k-1}\|\_{\infty} + s(s\_k + \mathbf{1})$$

$$\le s\_k s\_{k-1} (s\_{k-2} \|e\_{k-2}\|\_{\infty} + s) + s(s\_k + \mathbf{1}) = s\_k s\_{k-1} s\_{k-2} \|e\_{k-2}\|\_{\infty} + s(s\_k s\_{k-1} + s\_k + \mathbf{1})$$

$$\le s\_k s\_{k-1} s\_{k-2} (s\_{k-3} \|e\_{k-3}\|\_{\infty} + s) + s(s\_k s\_{k-1} + s\_k + \mathbf{1}) = s\_k s\_{k-1} s\_{k-2} s\_{k-3} \|e\_{k-3}\|\_{\infty} \tag{18}$$

$$\le s\_k s\_{k-1} s\_{k-2} (s\_{k-3} \|e\_{k-3}\|\_{\infty} + s) + s(s\_k s\_{k-1} + s\_k + \mathbf{1}) = s\_k s\_{k-1} s\_{k-2} s\_{k-3} \|e\_{k-3}\|\_{\infty}$$

$$+ s(s\_k s\_{k-1} s\_{k-2} + s\_k s\_{k-1} + s\_k + \mathbf{1}) \le \dots \le s\_k s\_{k-1} \cdot \dots \cdot s\_{k-r} \|e\_{k-r}\| $$

$$+ s(s\_k s\_{k-1} \cdot \dots \cdot s\_{k-r+1} + \dots + s\_k + \mathbf{1}).$$

Substituting into Eq. (18) r ¼ k - 1, we obtain

$$||\mathfrak{e}\_{k+1}||\_{\infty} \le \mathfrak{s}\_k \mathfrak{s}\_{k-1} \cdot \dots \cdot \mathfrak{s}\_1 ||\mathfrak{e}\_1|| + \mathfrak{s} (\mathfrak{s}\_k \mathfrak{s}\_{k-1} \cdot \dots \cdot \mathfrak{s}\_2 + \dots + \mathfrak{s}\_k + \mathfrak{1}) \le \mathbb{C}\_0 ||\mathfrak{e}\_0|| + O(\mathfrak{r}).$$

Q k From the second initial condition (6) it follows: k k <sup>e</sup><sup>1</sup> <sup>≤</sup>k k <sup>e</sup><sup>0</sup> and <sup>С</sup><sup>0</sup> <sup>¼</sup> sp. <sup>p</sup>¼<sup>1</sup> Now according to our assumption Ai-<sup>1</sup> ≥ Bi-1, which leads us to the relation

$$\tau \le \left( \frac{2\Gamma(2-\gamma\_{i-1})}{\lambda\Gamma(3-\beta\_{i-1})} \right)^{\frac{1}{\beta\_{i-1}-\gamma\_{i-1}}}, i = 2, \ldots, N-1. \tag{19}$$

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.

Mathematical Models of Oscillators with Memory DOI: http://dx.doi.org/10.5772/intechopen.81858

We consider the stability of an explicit finite-difference scheme (4). Suppose that Xk and two Yk are different solutions of the matrix Eq. (12) with initial conditions X<sup>0</sup> and Y0.

Theorem 2. An explicit finite-difference scheme (9) is conditionally stable if condition (14) is satisfied and the estimate holds jYk � Xkj≤C Yj <sup>0</sup> � X0j for any k, where С > 0 does not depend on the step τ.

Proof. We introduce the notation: ekþ<sup>1</sup> ¼ Ykþ<sup>1</sup> � Xkþ1: Then, Eq. (12) can be written in the form ekþ<sup>1</sup> ¼ Mek þ Fe,k. Here, as it was said in.

<sup>1</sup> <sup>T</sup> <sup>≤</sup> <sup>1</sup> Fe,k <sup>¼</sup> ðj <sup>f</sup>ðx1; tkÞ � <sup>f</sup>ðx1; tkÞj; …; <sup>j</sup> <sup>f</sup>ðxk; tkÞ � <sup>f</sup>ðxk; tkÞjÞ AkþBk AkþBk ðL1e1; …; LkekÞ ¼ ΔFkek

According to Theorem 1, we have the following estimate:

$$||M + \Delta F\_k|| \le \left( 3 + \frac{L}{A\_k + B\_k} \right) = s\_k.$$

Therefore, the following estimate holds

$$\|e\_{k+1}\|\_{\infty} \le \|M + \Delta F\_k\| \|\|e\_k\|\|\_{\infty} \le \left(\mathfrak{I} + \frac{L}{A\_k + B\_k}\right) \|e\_k\|\|\_{\infty}$$

$$\ell = s\_k \|\|e\_k\|\|\_{\infty} \le s\_k s\_{k-1} \|e\_{k-1}\|\|\_{\infty} \le s\_k s\_{k-1} s\_{k-2} \|e\_{k-2}\|\|\_{\infty} \le \dots \le s\_k s\_{k-1} \cdot \dots \cdot s\_{k-r} \|\|e\_{k-r}\|\|\_{\infty}$$

$$\text{With } r = k - 1, \text{ we obtain } \|\|e\_{k+1}\|\|\_{\infty} \le C\_0 \|\|e\_1\|\| \le C\_0 \|\|e\_0\|\| \text{ and } C\_0 = \prod\_{p=1}^k s\_p.$$

The last inequality follows the second condition of problem (6). Therefore, if X<sup>0</sup> there 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 μð Þ М ≫ 1 arise, 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: ξ ¼ max<sup>i</sup> jxi � x2<sup>i</sup>j—absolute error between the numerical solution xi in step τ and the numerical solution x2<sup>i</sup> in step τ=2. Then, the order of computational accuracy p can be estimated by the formula

$$p = \log\_2(\xi) / \log\_2(\pi/2).$$

We note that in the case when the fractional parameters in the scheme (9) do not change and have the following values of β<sup>k</sup> ¼ 2 and γ<sup>k</sup> ¼ 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

$$f(\mathbf{x}(t), t) = \delta \sin\left(\boldsymbol{\varrho}t\right) + t\boldsymbol{\varkappa}(t).$$

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

$$
\partial\_{0t}^{\beta(t)} \varkappa(\eta) + \lambda \partial\_{0t}^{\prime(t)} \varkappa(\eta) - t\varkappa(t) = \delta \cos\left(\eta t\right).
$$

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

$$\mathfrak{x}(\mathbf{0}) = \dot{\mathfrak{x}}(\mathbf{0}) = \mathbf{0}.$$

We note thatthe Airy oscillatoris used in opticsin the simulation of Airy laser beams. In this case, the explicit finite-difference scheme (9) has a more specific form:

$$\mathbf{x}\_0 = \mathbf{x}\_1 = \mathbf{0},$$

$$\mathbf{x}\_{k+1} = \frac{1}{A\_k + B\_k} ((2A\_k - k\tau)\mathbf{x}\_k - (A\_k - B\_k)\mathbf{x}\_{k-1}) \tag{20}$$

$$-\frac{A\_k}{A\_k + B\_k} \sum\_{j=1}^{k-1} a\_{j,k} \left(\mathbf{x}\_{k-j+1} - 2\mathbf{x}\_{k-j} + \mathbf{x}\_{k-j-1}\right)$$

$$-\frac{B\_k}{A\_k + B\_k} \sum\_{j=1}^{k-1} b\_{j,k} \left(\mathbf{x}\_{k-j+1} - \mathbf{x}\_{k-j-1}\right) + \delta \sin\left(qk\tau\right).$$

For the explicit finite-difference scheme (20), we choose the following values of the control parameters: T ¼ 1. λ ¼ 1, δ ¼ 5, φ ¼ 10, ω ¼ 10 and βðÞ¼ t 1:8 � 0:03 sin ð Þ ωt , γðÞ¼ t 0:8 � 0:05 cosð Þ ωt . And during the simulation, we will change the number of nodes N in the calculation grid.

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

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 p tends to unite.

This confirms that the explicit finite-difference scheme (9) and in particular the scheme (20) for Example 1 have 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 1 at 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).

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).

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.

$$\begin{array}{c c c c} \hline \text{N} & \text{g} & \text{p} & \text{p} \\ \hline \text{640} & 0.0003331017 & & 1.119146497 \\ \hline \text{1280} & 0.0001745618 & & 1.02636795 \\ \hline \text{2560} & 0.0000906971 & & 1.08981915 \\ \hline \end{array}$$

$$f(\mathfrak{x}(t),t) = \delta \sin\left(\mathfrak{q}t\right) - a\mathfrak{x}(t) + b\mathfrak{x}^3(t),$$

Table 1. Results of numerical simulation.

Mathematical Models of Oscillators with Memory DOI: http://dx.doi.org/10.5772/intechopen.81858

Figure 1.

The oscillogram (a) and the phase trajectory (b) for Example 1 with the parameter values T ¼ 10, N ¼ 1000:

Figure 2.

Oscillogram (a) and phase trajectory (b) for Example 1 with initial conditions xð0Þ ¼ 0:1, x\_ð0Þ ¼ 0:2, and δ ¼ 0.

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

$$\mathfrak{x}(\mathbf{0}) = \dot{\mathfrak{x}}(\mathbf{0}) = \mathbf{0}.$$

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

$$
\partial\_{\boldsymbol{\Omega}t}^{\beta(t)} \mathfrak{x}(\boldsymbol{\eta}) + \lambda \partial\_{\boldsymbol{\Omega}t}^{r(t)} \mathfrak{x}(\boldsymbol{\eta}) + b\boldsymbol{\varkappa}^{\beta}(\boldsymbol{t}) - a\boldsymbol{\varkappa}(\boldsymbol{t}) = \delta \sin \left( \boldsymbol{\varrho}t \right).
$$

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

#### Figure 3.

The oscillogram (a) and the phase trajectory (b) for Example 1 with initial conditions xð0Þ ¼ 0:1, x\_ð0Þ ¼ 0:2, and δ ¼ 0, λ ¼ 0.


#### Table 2.

Results of numerical simulation.

$$\mathbf{x}\_0 = \mathbf{x}\_1 = \mathbf{0},$$

$$\mathbf{x}\_{k+1} = \frac{1}{A\_k + B\_k} \left( (2A\_k + 1)\mathbf{x}\_k - \mathbf{x}\_k^3 - (A\_k - B\_k)\mathbf{x}\_{k-1} \right) \tag{21}$$

$$-\frac{A\_k}{A\_k + B\_k} \sum\_{j=1}^{k-1} a\_{j,k} \left( \mathbf{x}\_{k-j+1} - 2\mathbf{x}\_{k-j} + \mathbf{x}\_{k-j-1} \right)$$

$$-\frac{B\_k}{A\_k + B\_k} \sum\_{j=1}^{k-1} b\_{j,k} \left( \mathbf{x}\_{k-j+1} - \mathbf{x}\_{k-j-1} \right) + \delta \sin \left( \eta k \pi \right).$$

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 1 and 2. The results of numerical simulation for Example 2 are given in Table 2.

Mathematical Models of Oscillators with Memory DOI: http://dx.doi.org/10.5772/intechopen.81858

Figure 4. The oscillogram (a) and the phase trajectory (b) for Example 2.

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 p tends 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 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 λ ¼ 0 in the oscillatory system, we obtain an oscillogram and a phase trajectory as in Figure 6.

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

$$f(\mathbf{x}(t), t) = bt + c \sum\_{n=1}^{7} a\_n \sin \left( n\mathbf{x}(t) \right) - a \boldsymbol{\vartheta}^{(t)} \boldsymbol{\varkappa}(t), \tag{22}$$

where b is the spring travel speed; c is the surface adhesion energy; ω is the Ð <sup>1</sup> cosðπnτÞd<sup>τ</sup> frequency of free oscillations; and an <sup>¼</sup> <sup>2</sup><sup>n</sup> is the coefficients of the cosh <sup>2</sup>ðπτ<sup>Þ</sup> <sup>0</sup> 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

Oscillogram (a) and phase trajectory (b) for Example 3 with initial conditions xð0Þ ¼ 0:1, x\_ð0Þ ¼ 0:2, and δ ¼ 0.

Figure 6. Oscillogram (a) and phase trajectory (b) for Example 2 with initial conditions xð0Þ ¼ 0:1, x\_ð0Þ ¼ 0:2, and δ ¼ 0, λ ¼ 0.

Mathematical Models of Oscillators with Memory DOI: http://dx.doi.org/10.5772/intechopen.81858

#### Figure 7.

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

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 an in the expansion of the function (22). The values of these coefficients are taken from [43] a<sup>1</sup> ¼ 0:436, a<sup>2</sup> ¼ 0:344, a<sup>3</sup> ¼ 0:164, a<sup>4</sup> ¼ 0:058, a<sup>5</sup> ¼ 0:021, a<sup>6</sup> ¼ 0:004, and a<sup>7</sup> ¼ 0:003. Values of control parameters are βðÞ¼ t 1:8 � 0:03 sin ðπtÞ γðÞ¼ t 0:6 � 0:04 cosðπtÞ, N ¼ 3000 δ ¼ 50, τ ¼ 0:05, λ ¼ 0:3, b ¼ 1, ω ¼ 1, and xð0Þ ¼ 0, x\_ð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.

#### 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

#### Oscillators - Recent Developments

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 finitedifference 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:

$$
\partial\_{0t}^{\beta(\mathbf{x}(t),t)} \varkappa(\eta) + \lambda(\varkappa(t),t) \partial\_{0t}^{\gamma(\mathbf{x}(t),t)} \varkappa(\eta) = f(\varkappa(t),t), \\
\varkappa(\mathbf{0}) = a\_{0}, \dot{\varkappa}(\mathbf{0}) = a\_{1}.
$$

Another continuation of the research is related to the introduction of other memory functions K1ðt � τÞ, K2ðt � τÞ 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 ofthe stability ofthe rest points of dynamicalsystems described by the Cauchy problems (5) and (6) is also interesting, by analogy with the papers [47, 48].

#### Acknowledgements

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.

#### Author details

Roman Ivanovich Parovik1,2\*

1 Institute of Cosmophysical Research and Radio Wave Propagation FEB RAS, v. Paratunka, Kamchatskiy kray, Russia

2 Vitus Bering Kamchatka State University, Petropavlovsk-Kamchatskiy, Russia

\*Address all correspondence to: romanparovik@gmail.com

© 2018 The Author(s). Licensee IntechOpen. This chapteris distributed underthe terms oftheCreative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Mathematical Models of Oscillators with Memory DOI: http://dx.doi.org/10.5772/intechopen.81858

### References

[1] Volterra V. Sur les équations intégrodifférentielles et leurs applications. Acta Math. 1912;35(1):295-356

[2] Volterra V. Theory of Functionals and of Integral and Integro-differential Equations. New York: Dover Publication Inc; 2005. 226 p

[3] Mainardi F. Fractional relaxationoscillation and fractional diffusionwave. Chaos, Solitons & Fractals. 1996; 7(9):146-1477. DOI: 10.1016/0960-0779 (95)00125-5

[4] Meilanov R, Yanpolov M. Features of the phase trajectory of a fractal oscillator. Technical Physics Letters. 2002;28(1):30-32. DOI: 10.1134/ 1.1448634

[5] Achar B, Hanneken J, Clarke T. Response characteristics of a fractional oscillator. Physica A: Statistical Mechanics and its Applications. 2002; 309(3–4):275-288. DOI: 10.1016/ S0378-4371(02)00609-X

[6] Nakhusheva V. Differential Equations of Mathematical Models of Non-Local Processes. Moscow: Nauka; 2006. p. 173

[7] Al-Rabtah A, Ertürk V, Momani S. Solutions of a fractional oscillator by using differential transform method. Computers & Mathematics with Applications. 2010;59(3):1356-1362. DOI: 10.1016/j.camwa.2009.06.036

[8] Rand R, Sah S, Suchrsky M. Fractional Mathieu equation. Communications in Nonlinear Science and Numerical Simulation. 2010;15: 3254-3262. DOI: 10.1016/j. cnsns.2009.12.009

[9] Afanas'ev V, Daniel M. Polish stabilization of the inertial effects of the fractal oscillator. Technical Physics Letters. 2010;36(7):1-6

[10] Petras I. Fractional-Order Nonlinear Systems. Modeling, Analysis and Simulation. Beijing and Springer-Verlag. Berlin Heidelberg: Springer; 2011. 218 p

[11] Zurigat M. Solving fractional oscillators using laplace homotopy analysis method. Annals of the University of Craiova, Mathematics and Computer Science Series. 2011;38(4): 1-11

[12] Gomez-Aguilar JF, Rosales-Garc'ıab JJ, Bernal-Alvaradoa JJ, Cordova-Fraga T, Guzman-Cabrera R. Fractional mechanical oscillators. Revista Mexicana de F'ısica. 2012;58:348-352

[13] Duan J-S. The periodic solution of fractional oscillation equation with periodic input. Adv. Math. Phys. 2013; 2013:869484. DOI: 10.1155/2013/869484

[14] Parovik R. Fractal parametric oscillator as a model of a nonlinear oscillation system in natural mediums. International Journal of Communications, Network and System Sciences. 2013;6(3):134-138

[15] Xu Y, Agrawal O. Models and numerical solutions of generalized oscillator equations. Journal of Vibration and Acoustics. 2014;136:051005. DOI: 10.1115/1.4027241

[16] Syta A, Litak G, Lenci S, Scheffler M. Chaotic vibrations of the Duffing system with fractional damping. Chaos: An Interdisciplinary Journal of Nonlinear Science. 2014;24(1):013107. DOI: 10.1063/1.4861942

[17] Blaszczyk T. A numerical solution of a fractional oscillator equation in a nonresisting medium with natural boundary conditions. Romanian Reports in Physics. 2015;67(2):350-358

[18] Parovik R. Mathematical modeling of nonlocal oscillatory Duffing system with fractal friction. Bulletin KRASEC. Physical and Mathematical Sciences. 2015;10(1):16-21. DOI: 10.18454/ 2313-0156-2015-10-1-16-21

[19] Parovik R. Mathematical modeling of the hereditary oscillator. Computer Research and Modeling. 2015;7(5): 1001-1021

[20] Parovik R. On a credit oscillatory system with the inclusion of stick-slip. E3S Web of Conferences. 2016;11: 00018. DOI: 10.1051/e3sconf/ 20161100018

[21] Parovik R. Mathematical modelling of hereditarity airy oscillator with friction. Bulletin of South Ural State University. Series Mathematical Modelling, Programming & Computer Software. 2017;10(1):138-148. DOI: 10.14529/mmp170109

[22] Nakhushev A. Fractional Calculus and its Applications. Moscow: Fizmatlit; 2003. 272 p

[23] Kilbas A, Srivastava H, Trujillo J. Theory and Applications of Fractional Differential Equations. Amsterdam: Elsevier; 2006. 523 p

[24] Uchaikin V. Fractional Derivatives for Physicists and Engineers. Vol. I. Background and Theory. Beijing/Berlin: Higher Education Press/Springer-Verlag; 2013. 373 p

[25] Parovik R. Phase trajectories of the fractal parametric oscillator. Journal of Engineering Research and Applications. 2013;3(5):1520-1523. DOI: 10.4236/ ijcns.2013.63016

[26] Parovik R. Charts Strutt-Ince for generalized Mathieu equation. Bulletin KRASEC. Physical and Mathematical

Sciences. 2012;1(4):30-34. DOI: 10.18454/2079-6641-2012-4-1-29-30

[27] Drobysheva I. Mathematical modeling of nonlinear hereditary oscillators on the example of Duffing oscillator with fractional derivatives in the sense of Riemann-Liouville. Bulletin KRASEC. Physical and Mathematical Sciences. 2016;13(2): 39-45. DOI: 10.18454/2313-0156- 2016-13-2-39-45

[28] Kim V. Duffing oscillator with external harmonic action and variable fractional Riemann-Liouville derivative characterizing viscous friction. Bulletin KRASEC. Physical and Mathematical Sciences. 2016;13(2):46-49. DOI: 10.18454/2313-0156-2016-13-2-46-49

[29] Lipko O. Mathematical model of propagation of nerve impulses with regard hereditarity. Bulletin KRASEC. Physical and Mathematical Sciences. 2017;17(1):33-43. DOI: 10.18454/ 2313-0156-2017-16-1-52-60

[30] Novikova E. Van der Pol-Duffing oscillator with the effect of hereditary. Bulletin KRASEC. Physical and Mathematical Sciences. 2017;17(2): 65-75. DOI: 10.18454/2079-6641- 2017-18-2-65-75

[31] Zaitsev V, Ar K, Yarovoy G. Selfoscillations dynamics of active fractional oscillator. Theoretical Physics. 2013;14: 11-18

[32] Zaitsev V, Karlov A, Yarovoy G. Self-oscillations dynamics of fractional Thomson oscillator. Physics of Wave Processes and Radio Systems. 2012; 15(1):64-68

[33] Schroeder M. Fractals, Chaos, Power Laws: Minutes From an Infinite Paradise. New York: V.H. Freeman; 1991

[34] Gerasimov A. Generalization of linear deformation laws and their application to internal friction

Mathematical Models of Oscillators with Memory DOI: http://dx.doi.org/10.5772/intechopen.81858

problems. AS USSR. Applied Mathematics and Mechanics. 1948;12: 529-539

[35] Caputo M. Linear models of dissipation whose Q is almost frequency independent-II. Geophysical Journal International. 1967;13(5):529-539

[36] Parovik R. Finite-difference schemes for fractal oscillator with a variable fractional order. Bulletin KRASEC. Physical and Mathematical Sciences. 2015;11(2):85-92. DOI: 10.18454/2313-0156-2015-11-2-85-92

[37] Parovik R. Explicit finite-difference scheme for the numerical solution of the model equation of nonlinear hereditary oscillator with variable-order fractional derivatives. Archives of Control Sciences. 2016;26(3):429-435. DOI: 10.1515/acsc-2016-0023

[38] Xu Y, Erturk V. A finite difference technique for solving variable-order fractional integro-differential equations. Iranian Mathematical Society. 2014; 40(3):699-712

[39] Parovik R. Mathematical Modeling of Hereditary Linear Oscillators. Petropavlovsk-Kamchatsky: Vitus Bering Kamchatka State University; 2015. 178 p

[40] Parovik R. Fractional calculus in the theory of oscillating systems. Modern Science Technologies. 2017;1:61-68

[41] Daub EG, Carlson JM. Stick-slip instabilities and shear strain localization in amorphous materials. Physical Review E. 2009;80(6):066113

[42] Rekhviashvili S. Sh. Dimensional phenomena in condensed matter physics and nanotechnology. Nalchik, Russia: KBNts RAS; 2014

[43] Scholz Ch H. The mechanics of earthquakes and faulting. Cambridge: Cambridge university press; 2002

[44] Xiao-Jun Yang, Tenreiro Machado JA. A new fractional operator of variable order: Application in the description of anomalous diffusion. Physica A: Statistical Mechanics and its Applications. 2017;481:276-283

[45] Gao Feng, Xiao-Jun Yang. A new family of the local fractional PDEs. Fundamenta Informaticae. 2017; 151(1-4):63-75

[46] Boyadzhiev D, Hristo Kiskinov H, Veselinova M, Zahariev A. Stability analysis of linear distributed order fractional systems with distributed delays. Fractional Calculus and Applied Analysis. 2017;20(4):914-935

[47] Gallegos JA, Manuel A. Robustness and convergence of fractional systems and their applications to adaptive schemes. Fractional Calculus and Applied Analysis. 2017;20(4):895-913

[48] Xiao-Jun Yang, Tenreiro Machado JA, Cattani Carlo, Gao Feng. On a fractal LC-electric circuit modeled by local fractional calculus. Communications in Nonlinear Science and Numerical Simulation. 2017;47:200-206

Chapter 2

Abstract

Coşkun Deniz

the capability of these methods.

1. Introduction

form:

23

Quantum Harmonic Oscillator

Quantum harmonic oscillator (QHO) involves square law potential (x2

Keywords: Schrodinger equation, quantum mechanics, JWKB, MAF

<sup>H</sup>^ j i <sup>φ</sup> <sup>¼</sup> <sup>E</sup>j i <sup>φ</sup> ) �ℏ<sup>2</sup>

U rð Þ¼ <sup>1</sup> 2 mw<sup>2</sup> r <sup>2</sup> <sup>¼</sup> kr<sup>2</sup> 2m

Time-independent Schrodinger equation (TISE) is an eigenvalue problem in the

<sup>∇</sup><sup>2</sup> <sup>þ</sup> U rð Þ ˜ °

where the terms are in the usual meanings, namely, ∇2, the Laplacian operator; H^ , Hamiltonian operator (kinetic energy plus potential energy operators); m, mass; ℏ, Planck's constant divided by 2π; φ, wave function (eigenfunction); E, total energy (eigenvalue); and U(r), function of potential energy [1–7]. Quantum harmonic oscillator (QHO) is described by the TISE in (1) for the square law potential:

φ<sup>n</sup> ¼ Enφ (1)

≥ 0 (2)

2m

Schrodinger equation and is a fundamental problem in quantum mechanics. It can be solved by various conventional methods such as (i) analytical methods where Hermite polynomials are involved, (ii) algebraic methods where ladder operators are involved, and (iii) approximation methods where perturbation, variational, semiclassical, etc. techniques are involved. Here we present the general outcomes of the two conventional semiclassical approximation methods: the JWKB method (named after Jeffreys, Wentzel, Kramers, and Brillouin) and the MAF method (abbreviated for "modified Airy functions") to solve the QHO in a very good precision. Although JWKB is an approximation method, it interestingly gives the exact solution for the QHO except for the classical turning points (CTPs) where it diverges as typical to the JWKB. As the MAF method, it enables very approximate wave functions to be written in terms of Airy functions without any discontinuity in the entire domain, though, it needs careful treatment since Airy functions exhibit too much oscillatory behavior. Here, we make use of the parity conditions of the QHO to find the exact JWKB and approximate MAF solutions of the QHO within

) in the

#### Chapter 2

## Quantum Harmonic Oscillator

Coşkun Deniz

### Abstract

Quantum harmonic oscillator (QHO) involves square law potential (x<sup>2</sup> ) in the Schrodinger equation and is a fundamental problem in quantum mechanics. It can be solved by various conventional methods such as (i) analytical methods where Hermite polynomials are involved, (ii) algebraic methods where ladder operators are involved, and (iii) approximation methods where perturbation, variational, semiclassical, etc. techniques are involved. Here we present the general outcomes of the two conventional semiclassical approximation methods: the JWKB method (named after Jeffreys, Wentzel, Kramers, and Brillouin) and the MAF method (abbreviated for "modified Airy functions") to solve the QHO in a very good precision. Although JWKB is an approximation method, it interestingly gives the exact solution for the QHO except for the classical turning points (CTPs) where it diverges as typical to the JWKB. As the MAF method, it enables very approximate wave functions to be written in terms of Airy functions without any discontinuity in the entire domain, though, it needs careful treatment since Airy functions exhibit too much oscillatory behavior. Here, we make use of the parity conditions of the QHO to find the exact JWKB and approximate MAF solutions of the QHO within the capability of these methods.

Keywords: Schrodinger equation, quantum mechanics, JWKB, MAF

### 1. Introduction

Time-independent Schrodinger equation (TISE) is an eigenvalue problem in the form:

$$
\hat{H}|\rho\rangle = E|\rho\rangle \Rightarrow \left[\frac{-\hbar^2}{2m}\nabla^2 + U(r)\right]\rho\_n = E\_n\rho\tag{1}
$$

where the terms are in the usual meanings, namely, ∇2, the Laplacian operator; H^ , Hamiltonian operator (kinetic energy plus potential energy operators); m, mass; ℏ, Planck's constant divided by 2π; φ, wave function (eigenfunction); E, total energy (eigenvalue); and U(r), function of potential energy [1–7]. Quantum harmonic oscillator (QHO) is described by the TISE in (1) for the square law potential:

$$U(r) = \frac{1}{2}m\nu r^2 r^2 = \frac{kr^2}{2m} \ge 0\tag{2}$$

pffiffiffiffiffiffiffiffiffi where <sup>w</sup> <sup>¼</sup> <sup>k</sup>=<sup>m</sup> is the natural angular momentum (associated with the angular frequency f ¼ w=ð Þ 2π ). Since U rð Þ≥ 0, our eigenvalue problem (or bound-state problem) requires En ≥0 to give the following:

$$\nabla^2 \rho\_n(E\_n \ge 0, r) + f(E\_n \ge 0, r)\rho\_n(E\_n \ge 0, r) = 0;$$

$$f(E\_n \ge 0, r) = \hbar^2(E\_n \ge 0, r) = \frac{2m}{\hbar^2}[E\_n - U(r)] = \frac{2m}{\hbar^2} \left[E\_n - \frac{1}{2}m\omega^2 r^2\right] =: \begin{cases} 0 \text{ or } (+) \\ (-) \end{cases}$$

The QHO is a very good approximation in solving systems of diatomic molecules vibrating under the spring constant [1, 2, 5] and finds various modern physics applications such as in [8–10] as stated in a famous quotation: "the career of a young theoretical physicist consists of treating the harmonic oscillator in ever-increasing levels of abstraction by Sidney Coleman" [10, 11]. Here, U(r) is central potential which can be given in Cartesian coordinates (x, y, z) where solutions involve Hermite polynomials as in [1–7, 12] or in spherical coordinates (r, θ, ϕ) where solutions involve spherical harmonics as in [1, 2, 13]. For simplicity, it is widely studied in one dimension (say, in x only), and higher dimensional systems are called isotropic harmonic oscillators in 2D or in 3D. The QHO can be solved by various conventional methods such as the following: (i) by analytical methods where some analytic functions involving Hermite polynomials are involved [1–5]; (ii) algebraic methods where ladder operators are involved, that is, [1, 2]; and (iii) by approximation methods such as perturbation methods, JWKB method, variational methods, etc., that is, [1–6, 14]. Brownian study of QHO as an open dynamic quantum system in terms of quantum Langevin equation was studied in [15–17]. We study here one dimensional and non-frictional, that is, undamped case, and present its solution by the two following conventional semiclassical approximation methods: (i) the JWKB method (named after the authors, Jeffreys, Wentzel, Kramers, and Brillouin, who contributed to the theory) [1–7, 14] and (ii) the MAF method (abbreviated from modified Airy function) [3, 18–23].

JWKB method is known to give exact eigenenergies for the QHO, but eigenfunctions fail at and around the classical turning points (CTPs) where f ¼ 0 (or, equivalently, En ¼ U rð Þ) in (3) as typical to the JWKB method [1–7, 14]. These discontinuities prevent us from using continuity at the boundaries by equating the JWKB solutions of two neighboring regions directly at the CTPs to find the eigenenergy-dependent coefficients in the general JWKB eigenfunctions (wave functions). It also prohibits the use of normalizability of the eigenfunctions between �∞ and ∞. To surmount the problem, parity conditions of the problem regarding the symmetry of the QHO in the dimensionless form are used, and advanced computational software such as Mathematica can be used to achieve these calculations [3, 4, 14, 24]. Moreover, asymptotic matching is required in the JWKB solutions to maintain the normalizability except for the CTPs as discussed above [4, 7, 14]. As to the MAF method, it does not exhibit discontinuities at the CTPs, though highly oscillating behavior of the Airy functions requires careful handling in finding their zeros and the parity treatment used in the JWKB solution seems straightforward to be also applicable to the MAF solution of the QHO [3, 19–22]. Although it was originally suggested in 1931 by Langer in [25], finding zeros of highly oscillatory Airy functions became practical as the advances in computational software and the MAF method became widespread by the 1990s [3, 18–23]. In this work, we present the general outcomes of the conventional JWKB and MAF methods as two semiclassical conventional methods and solve the QHO by using the parity condition of the problem in the dimensionless form pedagogically. We will discuss the treatment of parity matching and asymptotic matching in solving the QHO by these semiclassical methods.

#### 2. Exact solution of the QHO in 1D by the analytic method

The QHO in (3) is a bound-state problem which can be written in 1D for the 2 2 <sup>¼</sup> kx<sup>2</sup> potential function in 1D (U xð Þ ¼ <sup>2</sup> <sup>1</sup> mw x ≥ 0) as follows: <sup>2</sup><sup>m</sup>

$$\rho\_n^{'}(\mathbf{x}) + f(E\_n \ge \mathbf{0}, \mathbf{x})\rho\_n(\mathbf{x}) = \mathbf{0}$$

$$f(E\_n \ge \mathbf{0}, \mathbf{x}) = \hbar^2 (E\_n \ge \mathbf{0}, \mathbf{x}) = \frac{2m}{\hbar^2} [E\_n - U(\mathbf{x})] = \frac{2m}{\hbar^2} \left[ E\_n - \frac{1}{2}mw^2\mathbf{x}^2 \right] =: \begin{cases} \mathbf{0} \, or \, (+) \\ (-) \end{cases}$$

or, simply,

$$
\rho\_n^\*(\mathbf{x}) + \frac{2m}{\hbar^2} \left( E\_n - \frac{1}{2} m w^2 \mathbf{x}^2 \right) \rho\_n(\mathbf{x}) = \mathbf{0} \tag{5}
$$

whose solution by various conventional approaches (such as analytical, algebraic, approximation, etc.) is given in any fundamental textbooks, that is, [1–3, 5] and whose results can be summarized as follows [14]:

i. Change of variable in (4) and (5):

$$
\gamma(\mathbf{x}) = \beta \mathbf{x} = \sqrt{\frac{mw}{\hbar}} \mathbf{x} \tag{6}
$$

ii. TISE for the QHO in 1D in dimensionless form:

$$\rho\_n^\*(\boldsymbol{y}) + k^2[\lambda\_n(\boldsymbol{E}\_n \ge \mathbf{0}), \boldsymbol{y}]\rho\_n(\boldsymbol{y}) = \mathbf{0}$$

$$\Leftrightarrow \left\{ f[\lambda\_n(\boldsymbol{E}\_n) \ge \mathbf{0}, \boldsymbol{y}] =: k^2[\lambda\_n(\boldsymbol{E}\_n) \ge \mathbf{0}, \boldsymbol{y}] = \lambda\_n^2 - \boldsymbol{y}^2, \lambda\_n^2 = \frac{2E\_n}{\hbar w} =: \begin{cases} \mathbf{0} \, or \, (+) \\ (-) \end{cases} \right\} \tag{7}$$

Note that here <sup>f</sup> <sup>≕</sup> <sup>k</sup><sup>2</sup> is <sup>a</sup> function of <sup>λ</sup>nð Þ En &<sup>y</sup> and <sup>λ</sup><sup>n</sup> <sup>≥</sup> <sup>0</sup> since En <sup>≥</sup>0. Moreover, fðλn; yÞ is an even function as shown in Figure 1 (λ is chosen as continuous including the discrete energy values assuming that the eigenenergies have not been found yet).

iii. Exact eigenenergies:

$$\Lambda\_{EX} = \colon \Lambda = \lambda\_n^2 = 2n + 1 \Rightarrow \left| E\_n = \left( n + \frac{1}{2} \right) \hbar w, n = 0, 1, 2, \dots \right| \tag{8}$$

iv. Exact eigenfunctions (wave functions) in y:

$$\varphi\_n(\boldsymbol{y}) = \boldsymbol{\nu}(\boldsymbol{\lambda}\_n, \boldsymbol{y}) = \sqrt{\frac{1}{\sqrt{\pi} 2^n n!}} H\_n(\boldsymbol{y})^{-\frac{\boldsymbol{\gamma}^2}{2}}, n = \mathbf{0}, 1, 2, \dots \tag{9}$$

v. By using (6), we have the wave functions in x:

Figure 1. Graphs of fðλ ≥0; yÞ and gðλ≥0; yÞ ¼ 0.

$$\rho\_n(\infty) = \circ \rho(\beta, \lambda\_n, \infty) = \sqrt{\frac{\beta}{\sqrt{\pi} 2^n n!}} H\_n(\beta \infty)^{-\frac{(\beta \alpha)^2}{2}}, n = 0, 1, 2, \dots \tag{10}$$

We used two different symbols (φ and ψ) to label the functions in two different independent variables (in x and in y, respectively). Exact wave functions in (9) (via exact eigen energies in (8)) are given for even and odd n values in Figures 3 and 4 along with the JWKB solutions for comparison. Hn in (9) and (10) is Hermite polynomials with indice n (named after the French mathematician Charles


Table 1. Some properties of Hermite polynomials. Hermite). Some of the properties of Hermite polynomials are tabulated in Table 1, and calculation of conversion factor β which exhibits a quantization, namely,

$$\beta = \sqrt{\frac{m\omega}{\hbar}}\Big|\_{w=\sqrt{\frac{k}{m}}} = \left(\frac{mk}{\hbar^2}\right)^{1/4} = \frac{k}{\hbar w} = \sqrt{\frac{(n+1/2)k}{E\_n}} =: \beta\_n. \tag{11}$$

is given along with the related Mathematica codes in [14].

### 3. A review of the JWKB solution of the QHO

2D plot of Figure 1 is schematically given in Figure 2 for the QHO under study (in the dimensionless form) from which we have the following outcomes [14]:

#### 3.1 JWKB eigenenergies of the QHO

JWKB eigenenergies can be found by applying the Bohr-Sommerfeld quantization formula given by [1–7, 14]:

$$\int\_{\mathcal{Y}\_1}^{\mathcal{Y}\_2} k\left[\tilde{\lambda}\_f(E\_n), \mathbf{y}\right] d\mathbf{y} = \left(n + \frac{1}{2}\right)\pi \tag{12}$$

as follows:

$$\begin{split} \int\_{\gamma 1 = -\lambda\_n}^{\gamma\_1 = -\lambda\_n} k\left[\bar{\lambda}\_f(E\_n), \mathbf{y}\right] d\mathbf{y} &= 2 \int\_0^{\gamma\_1 = \lambda\_n} k\left[\bar{\lambda}\_f(E\_n), \mathbf{y}\right] d\mathbf{y} = \left(n + \frac{1}{2}\right)\pi \\ \Rightarrow \int\_{-\lambda\_n}^{\lambda\_n} \sqrt{\bar{\lambda}\_f^2(E\_n) - \mathbf{y}^2} d\mathbf{y} &= 2 \int\_0^{\lambda\_n} k\left[\bar{\lambda}\_f(E\_n), \mathbf{y}\right] d\mathbf{y} = \left(n + \frac{1}{2}\right)\pi \\ \mathbf{y} &= \sin \theta \Rightarrow 2 \int\_{-\pi/2}^{\pi/2} \bar{\lambda}\_f^2(E\_n) \cos^2 \theta d\theta = \left(n + \frac{1}{2}\right)\pi \Leftrightarrow \bar{\lambda}\_f^2(E\_n) \frac{\pi}{2} = \left(n + \frac{1}{2}\right)\pi \\ \Rightarrow \boxed{\bar{\lambda}\_f(E\_n) = \sqrt{2n + 1} \left(\text{or} \right) \bar{E}\_n &:= E\_{f|WKB,n} = \left(n + \frac{1}{2}\right)\hbar \mathbf{n} \end{split} \tag{13}$$

Figure 2. Schematic 2D sketch of <sup>f</sup>½λnð Þ <sup>≥</sup>0; <sup>y</sup>� <sup>≕</sup> <sup>k</sup><sup>2</sup> En <sup>½</sup>λnðEn <sup>≥</sup>0Þ; <sup>y</sup>� for <sup>a</sup> given <sup>λ</sup>n.

which is already the exact solution given in (8) [1–4, 6, 7, 14]. Results are given along with the MAF solutions in Table 2 for comparison. Note that we use the following notation for the symmetrical (or even parity (EP)) and antisymmetrical (or odd parity (OP)) solutions:

$$\begin{aligned} \tilde{E}\_{fn} = \begin{cases} \tilde{E}\_{f,n\_t} = (n\_t + \mathbf{1}/2)\pi, n\_t = \mathbf{0}, 2, 4, \dots \\ \tilde{E}\_{f,n\_d} = (n\_d + \mathbf{1}/2)\pi, n\_d = \mathbf{1}, 3, 5, \dots \end{cases} \end{aligned} \tag{14}$$

where the subscripts "J, ns" represent J, JWKB, and ns, symmetrical indices (ns = even), and similarly, "J, na" represents J, JWKB, and na, antisymmetrical indices (ns ¼odd).

#### 3.2 JWKB solution of eigenfunctions (wave functions) of the QHO

Conventional first-order JWKB solution of the QHO given in the normal form in (4) or (7) is as follows:

$$\rho\_f(\lambda\_n, y) = \frac{c\_{f1} \exp\left[-i \int^{\eta\_{\nu} \cdot y} k(\lambda\_n, y) dy\right]}{\sqrt{k(\lambda\_n, y)}} + \frac{c\_{f2} \exp\left[i \int^{\eta\_{\nu} \cdot y} k(\lambda\_n, y) dy\right]}{\sqrt{k(\lambda\_n, y)}}\tag{15}$$

where yt is either of the classical turning points (CTPs: either "y1, on the left" or "y2, on the right" depending on the region under question) and cJ1&cJ<sup>2</sup> are arbitrary JWKB constants. Once solution in any region is found (say, φJII), solution in the adjacent region (say, φJIII) can directly be found by the conventional JWKB connection formulas given in [3, 4, 14] without calculating it via (15). The integrals here are the definite integrals whose upper and lower values should be chosen as the related turning point (either y<sup>1</sup> or y2) and the variable y should be in the correct ascending integration order. Normally, constant coefficients in the general solutions are determined from normalization by applying the boundary conditions of the


Table 2. JWKB and MAF eigenenergies.

#### Quantum Harmonic Oscillator DOI: http://dx.doi.org/10.5772/intechopen.85147

eigenvalue problem. However, it is useless since the boundary conditions correspond to the CTPs at which (and also in a narrow region around) the conventional first-order JWKB solutions typically diverge [3, 4, 7, 14]. This might be thought as a violation of continuity requirement of the acceptable wave function properties concerning continuity, but higher order JWKB approximation can fix it. Now, due to the discontinuities at the boundaries between the adjacent regions (such as between I and II or between II and III), the unidirectional JWKB connection formulas (surely, for the first-order JWKB) given in the literature [3, 4, 7, 14] cannot be used to find the constant coefficients in the general solution in (18) and (19). Note that these connection formulas can be used to determine the structure of the JWKB solutions in all regions (I, II, and III), but they cannot be used to find the constant coefficients (which will be a function of eigenenergy) as explained.

However, we are fortunately not helpless: since <sup>f</sup>½<sup>λ</sup> ð Þ≥0; <sup>y</sup>� <sup>≕</sup> <sup>k</sup><sup>2</sup> ½λ ð ≥0Þ; y� in (7) is an even function (see Figure 2), we should have even and odd-parity solutions. If we start by considering the exact solutions in (9) and (10) and considering them to be approximate to the JWKB solution (shown with tilde and subscript J), we have the following outcomes [14]: <sup>n</sup> En <sup>n</sup> En

� � � � � � � pffiffiffi ψ λ<sup>ð</sup> <sup>n</sup>; <sup>y</sup>Þ ¼ <sup>½</sup>φ β<sup>ð</sup> ; <sup>λ</sup>n; <sup>x</sup> ! <sup>y</sup>=βÞ� <sup>¼</sup><sup>1</sup>OR : φ β<sup>ð</sup> ; <sup>λ</sup>n; <sup>x</sup>Þ ¼ βψ λ<sup>ð</sup> <sup>n</sup>; <sup>y</sup> ! <sup>β</sup>x<sup>Þ</sup> <sup>β</sup> (16) 8 >>>>>>>>< 8 < iÞφnðβ; λ; 0޼� p, p . 0 E:P: : φnðβ; λ; �xÞ ¼ φnðβ; λ; xÞ, n ¼ 0, 2, 4, … ) : <sup>ð</sup>β;λ;x<sup>Þ</sup> <sup>n</sup> ii<sup>Þ</sup> <sup>∂</sup><sup>φ</sup> <sup>¼</sup> <sup>0</sup> <sup>∂</sup><sup>x</sup> <sup>x</sup>¼<sup>0</sup> 8 <sup>&</sup>lt; >>>>>>>>: iÞ λnðβ; λ; 0Þ ¼ 0 O:P: : φnðβ; λ; �x޼�φnðβ; λ; xÞ, n ¼ 1, 3, 5, … ) : <sup>ð</sup>β;λ;x<sup>Þ</sup> <sup>n</sup> ii<sup>Þ</sup> <sup>∂</sup><sup>φ</sup> ¼ � q, <sup>q</sup> . <sup>0</sup> <sup>∂</sup><sup>x</sup> <sup>x</sup>¼<sup>0</sup> 3 7 7 5 2 6 6 4 8 >>>>>>>>>< >>>>>>>>>: �<sup>2</sup> �<sup>p</sup> <sup>β</sup><sup>~</sup> <sup>≃</sup>β<sup>J</sup>:E:P: <sup>¼</sup> , n <sup>¼</sup> <sup>0</sup>, <sup>2</sup>, <sup>4</sup>, … <sup>φ</sup> <sup>ð</sup>1; <sup>λ</sup>; <sup>y</sup><sup>Þ</sup> <sup>n</sup> <sup>y</sup>¼<sup>0</sup> <sup>2</sup> <sup>β</sup> <sup>≃</sup>β<sup>~</sup> <sup>J</sup> <sup>¼</sup> �<sup>q</sup> <sup>β</sup><sup>~</sup> <sup>O</sup>:P: <sup>≃</sup>β<sup>J</sup>:O:P: <sup>¼</sup> , n <sup>¼</sup> <sup>1</sup>, <sup>3</sup>, <sup>5</sup>, … <sup>∂</sup><sup>φ</sup> <sup>ð</sup>1; <sup>λ</sup>; <sup>y</sup><sup>Þ</sup> <sup>n</sup> <sup>∂</sup><sup>y</sup> <sup>y</sup>¼<sup>0</sup> ) (17)

where p&q are positive real constants regarding the even-parity (EP) and oddparity (OP) initial values of the physical system. Remember that we use φ for the xsystem and ψ for the y-system as shown in (16). In finding the constant coefficients, we can take �p ¼ �q ¼ 1, and alternating sign can be modified as a parity matching as follows [14]:

$$
\bar{\boldsymbol{\Psi}}\_{(J,E,\mathbb{P})}^{\left(\text{par.m.}\right)}(\lambda\_{n},\boldsymbol{y}) = (-1)^{\binom{4}{2}} \times \begin{cases}
\bar{\boldsymbol{\Psi}}\_{f,I}^{\left(\text{asy.m.}\right)}(\lambda\_{n},\boldsymbol{y}) = \bar{\boldsymbol{\Psi}}\_{f,\text{III}}^{\left(\text{asy.m.}\right)}(\lambda\_{n},-\boldsymbol{y}) \text{ for } -\infty \le \boldsymbol{y} \le -\lambda\_{n} \\
\bar{\boldsymbol{\Psi}}\_{f,\text{II}}(\lambda\_{n},\boldsymbol{y}) = \begin{cases}
\bar{\boldsymbol{\Psi}}\_{f,\text{II}}^{\left(\lambda\_{n},-\boldsymbol{y}\right)}(\lambda\_{n},-\boldsymbol{y}) \text{ for } -\lambda\_{n} < \boldsymbol{y} \le \mathbf{0} \\
\bar{\boldsymbol{\Psi}}\_{f,\text{II}}^{\left(\lambda\_{n},\boldsymbol{y}\right)}(\lambda\_{n},\boldsymbol{y}) \text{ for } 0 < \boldsymbol{y} \le \lambda\_{n} \\
\bar{\boldsymbol{\Psi}}\_{f,\text{III}}^{\left(\lambda\_{n},\boldsymbol{y}\right)}(\lambda\_{n},\boldsymbol{y}) \text{ for } \lambda\_{n} \le \boldsymbol{y} < \infty \\
\end{cases}
\end{cases} \tag{18}$$

$$
\bar{\boldsymbol{\Psi}}\_{(J,O,\mathcal{P})}^{\left(\text{par.m.}\right)}(\lambda\_{n},\boldsymbol{y}) = (-1)^{\left(\frac{\boldsymbol{\lambda}-\boldsymbol{1}\right)}{\boldsymbol{\lambda}}} \times \begin{cases}
\bar{\boldsymbol{\Psi}}\_{f,I}^{\left(\text{asy.m.}\right)}(\lambda\_{n},\boldsymbol{y}) = -\bar{\boldsymbol{\Psi}}\_{f,\text{III}}^{\left(\text{asy.m.}\right)}(\lambda\_{n},-\boldsymbol{y}) \text{ for } -\infty < \boldsymbol{y} \le -\boldsymbol{\lambda}\_{n} \\
\bar{\boldsymbol{\Psi}}\_{f,\text{II}}(\lambda\_{n},\boldsymbol{y}) = \begin{cases}
\bar{\boldsymbol{\Psi}}\_{f,\text{II}}(\lambda\_{n},\boldsymbol{y}) \text{ for } 0 < \boldsymbol{y} \le \boldsymbol{\lambda}\_{n} \\
\bar{\boldsymbol{\Psi}}\_{f,\text{III}}^{\left(\text{asy.m.}\right)}(\lambda\_{n},\boldsymbol{y}) \text{ for } \lambda\_{n} \le \boldsymbol{y} < \mathbf{0}
\end{cases}
\end{cases}
\tag{19}
$$

� � � � where the superscripts (par.m.) and (asy.m.) represent parity matched and asymptotically matched JWKB solutions, respectively. Eqs. (18) and (19) tells that we will take �p ¼ �q ! 1 to find the solutions in 0≤y≤ y<sup>2</sup> ¼ λ<sup>n</sup> ∪ y<sup>2</sup> ¼ λ<sup>n</sup> ≤y , ∞ firstly by using (21) for the asymptotic matching and then extending it to the second quadrant according to the parity under question. Note that asymptotically matched general ðJWKBÞ<sup>1</sup> solution can be obtained as follows (see [3, 4, 7, 14] for details):

$$\tilde{\boldsymbol{\rho}}\_{f}^{(\text{ay.m.})}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y}) = \begin{cases} \tilde{\boldsymbol{\rho}}\_{f\text{I}}^{(\text{ay.m.})}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y}) = \text{either } \tilde{k}\_{1}\tilde{\boldsymbol{\rho}}\_{1}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y} < \boldsymbol{y}\_{t1}) \text{ or } \tilde{k}\_{2}\tilde{\boldsymbol{\rho}}\_{2}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y} < \boldsymbol{y}\_{t1})\\ \tilde{\boldsymbol{\rho}}\_{f\text{II}}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y}) = \tilde{\boldsymbol{\rho}}\_{f}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y}\_{t1} < \boldsymbol{y} < \boldsymbol{y}\_{t2})\\ \tilde{\boldsymbol{\rho}}\_{f\text{III}}^{(\text{ay.m.})}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y}) = \text{either } \tilde{k}\_{1}\tilde{\boldsymbol{\rho}}\_{f1}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y}\_{t2} < \boldsymbol{y}) \text{ or } \tilde{k}\_{2}\tilde{\boldsymbol{\rho}}\_{f2}(\boldsymbol{\lambda}\_{\text{v}},\boldsymbol{y}\_{t2} < \boldsymbol{y}) \end{cases} \tag{20}$$

so that they exhibit the following asymptotic behaviors:

$$\tilde{\rho}\_{\mathcal{I}}^{(\text{ay},\text{m})}(\lambda\_{n},\mathcal{y}) = \begin{cases} \lim\_{\mathcal{Y}\to\to\infty} \left[ \tilde{\rho}\_{\mathcal{I}} \left( \lambda\_{n},\mathcal{Y} \le \mathcal{Y}\_{t\_{1}} \right) = \tilde{\rho}\_{\mathcal{I}\mathcal{I}} (\lambda\_{n},\mathcal{y}) \right] = \mathbf{0} \\\lim\_{\mathcal{Y}\to\infty} \left[ \tilde{\rho}\_{\mathcal{I}} \left( \lambda\_{n},\mathcal{Y}\_{t2} \le \mathcal{Y} \right) = \tilde{\rho}\_{\mathcal{I}\text{III}} (\lambda\_{n},\mathcal{y}) \right] = \mathbf{0} \end{cases} \tag{21}$$

#### 3.2.1 Even-parity (EP) wave functions

When initial values at x ¼ y ¼ 0 for the EP case in (17), namely (by using (16)),

$$\left\{\ddot{\boldsymbol{\nu}}\_{\boldsymbol{f},\boldsymbol{II}}(\boldsymbol{\lambda}\_{\boldsymbol{n}},\boldsymbol{y})\Big|\_{\boldsymbol{\mathcal{Y}}=\boldsymbol{0}}=\mathbf{1}/\sqrt{\boldsymbol{\beta}};\Big{\boldsymbol{\partial}}\_{\boldsymbol{\mathcal{Y}}}\ddot{\boldsymbol{\mu}}\_{\boldsymbol{f},\boldsymbol{II}}(\boldsymbol{\lambda}\_{\boldsymbol{n}},\boldsymbol{y})\Big|\_{\boldsymbol{\mathcal{Y}}=\boldsymbol{0}}=\mathbf{0}\right\},\tag{22}$$

is applied to the JWKB solution in (15), we find the following:

$$(i)\ \partial\_{\boldsymbol{\eta}}\tilde{\boldsymbol{\eta}}\_{\boldsymbol{I},\boldsymbol{I}}(\lambda\_{n},\boldsymbol{y})\big|\_{\boldsymbol{y}=\boldsymbol{0}}=\mathbf{0}\Rightarrow\tilde{\boldsymbol{\eta}}\_{\boldsymbol{II}}(\tilde{\lambda}\_{n}=\lambda\_{n},\boldsymbol{y})=\frac{\boldsymbol{A}(\lambda\_{n})}{\sqrt{\mathbb{K}(\lambda\_{n},\boldsymbol{y})}}\cos\left[\int\_{0}^{\boldsymbol{I}}\mathbb{k}(\lambda\_{n},\boldsymbol{y})d\boldsymbol{y}\right],\tag{23}$$

where the second complementary solution (in the sine form) has been canceled and calculation of the integral in the cosine term can be calculated by the similar change of variable as in (13) whose result will give ηð y; 0Þ (see Eq. (18) below and apply ηðy; 0Þ ¼ η λð <sup>n</sup> ! y; y ! 0Þ).

$$\left. \left( i \right) \left. \ddot{\boldsymbol{\nu}}\_{l, \Pi} (\lambda\_n, \boldsymbol{\chi}) \right|\_{\boldsymbol{\chi} = 0} = \mathbf{1} / \sqrt{\boldsymbol{\beta}} \Rightarrow \mathbf{A} (\lambda\_n) = \sqrt{\frac{\lambda\_n}{\boldsymbol{\beta}}} \tag{24}$$

Quantum Harmonic Oscillator DOI: http://dx.doi.org/10.5772/intechopen.85147

and by using (16), we have.

$$\left(\tilde{\boldsymbol{\mu}}\_{\boldsymbol{f},\Pi}(\boldsymbol{\beta},\lambda\_{n},\mathbf{x}) = \sqrt{\boldsymbol{\beta}} \tilde{\boldsymbol{\mu}}\_{\boldsymbol{f},\Pi}(\lambda\_{n},\boldsymbol{y} \to \boldsymbol{\beta}\mathbf{x}) = \sqrt{\frac{\lambda\_{n}}{k(\lambda\_{n},\boldsymbol{\beta}\mathbf{x})}} \cos\left[\int\_{0}^{\boldsymbol{\beta}\mathbf{x}} k(\lambda\_{n},\mathbf{y}) d\mathbf{y}\right], \text{for } \mathbf{0} \le \mathbf{x} \le \lambda\_{n}/\boldsymbol{\beta}. \tag{25}$$

Now, by applying the JWKB connection formula with a small phase term α, we get.

$$\bar{\boldsymbol{\nu}}\_{\text{III}}(\bar{\boldsymbol{\lambda}}\_{\text{n}} = \boldsymbol{\lambda}\_{\text{n}}, \boldsymbol{y}) = \frac{\boldsymbol{A}(\boldsymbol{\lambda}\_{\text{n}})}{\sqrt{\kappa(\boldsymbol{\lambda}\_{\text{n}}, \boldsymbol{y})}} \begin{cases} \cos\left[a(\boldsymbol{\lambda}\_{\text{n}})\right] \exp\left[\zeta(\boldsymbol{\lambda}\_{\text{n}}, \boldsymbol{y})\right] \\ + \frac{1}{2} \sin\left[a(\boldsymbol{\lambda}\_{\text{n}})\right] \exp\left[-\zeta(\boldsymbol{\lambda}\_{\text{n}}, \boldsymbol{y})\right] \end{cases} \text{, for } \boldsymbol{\lambda}\_{\text{n}} < \boldsymbol{y} < \boldsymbol{\infty}, \tag{26}$$

and the asymptotically matched (modified) wave function in region III via (20) and (21) of [3, 4, 7, 14] gives:

$$\tilde{\boldsymbol{\psi}}\_{\text{III}}^{(\text{ay.}\text{m},\text{})} \left( \tilde{\boldsymbol{\lambda}}\_{n} = \boldsymbol{\lambda}\_{n}, \boldsymbol{y} \right) = \frac{\boldsymbol{A} \left( \boldsymbol{\lambda}\_{n} \right)}{2 \sqrt{\kappa \left( \boldsymbol{\lambda}\_{n}, \boldsymbol{y} \right)}} \sin \left[ a \left( \boldsymbol{\lambda}\_{n} \right) \right] \exp \left[ -\zeta \left( \boldsymbol{\lambda}\_{n}, \boldsymbol{y} \right) \right], \text{for } \boldsymbol{\lambda}\_{n} < \boldsymbol{y} < \infty \tag{27}$$

Abbreviations we use for the EP JWKB solutions here (and also for the OP solutions in the next subsection) are as follows [14]:

$$\begin{cases} & a(\lambda\_n) = \int\_0^{\lambda\_n} k(\lambda\_n, y) dy + \frac{\pi}{4} = \eta(\lambda\_n, 0) + \frac{\pi}{4} \\ & \quad \cdot \left(\eta(\lambda\_n, y) = \int\_y^{\lambda\_n} k(\lambda\_n, y) dy \to -\frac{\lambda\_n^2 \pi}{4} - \frac{\lambda\_n^2}{2} \sin^{-1}\left(\frac{y}{\lambda\_n}\right) - \frac{y}{2} \sqrt{\lambda\_n^2 - y^2} \\ & \quad \cdot \left(\lambda\_n, y\right) = \int\_y^{\lambda} \kappa(\lambda\_n, y) dy \to -\frac{y}{2} \sqrt{y^2 - \lambda\_n^2} - \frac{1}{2} \lambda\_n^2 \ln \left| \frac{y + \sqrt{y^2 - \lambda\_n^2}}{\lambda\_n} \right| \end{cases} \tag{28}$$

� <sup>ψ</sup><sup>ð</sup>asy:m:<sup>Þ</sup> Since we have already calculated <sup>ψ</sup><sup>~</sup> <sup>ð</sup>β; <sup>λ</sup>n; <sup>x</sup><sup>Þ</sup> and <sup>~</sup> J,II III <sup>~</sup><sup>λ</sup> <sup>¼</sup> <sup>λ</sup>n; <sup>y</sup><sup>Þ</sup> in the first <sup>n</sup> quadrant (0≤ y≤λn), JWKB solutions in the other regions can be easily written as in (18). JWKB wave functions regarding the EP case are given in Figure 3 along with the exact solutions for comparison.

#### 3.2.2 Odd-parity (OP) wave functions

Similarly, by using the boundary conditions for the OP case in (17), namely (by using (16)),

$$\left. \left\{ \tilde{\boldsymbol{\nu}}\_{f,\Pi}(\boldsymbol{\lambda}\_{\boldsymbol{n}}, \boldsymbol{\chi}) \Big|\_{\boldsymbol{\mathcal{Y}}=\boldsymbol{0}} = \mathbf{0} ; \left. \partial\_{\boldsymbol{\mathcal{Y}}} \tilde{\boldsymbol{\nu}}\_{f,\Pi}(\boldsymbol{\lambda}\_{\boldsymbol{n}}, \boldsymbol{\mathcal{y}}) \right|\_{\boldsymbol{\mathcal{Y}}=\boldsymbol{0}} = \mathbf{1} / \sqrt{\boldsymbol{\beta}^{3}} \right\},\tag{29}$$

and starting with region II.

Figure 3. Exact and JWKB solutions of EP wave functions (for p ¼ 1).

Quantum Harmonic Oscillator DOI: http://dx.doi.org/10.5772/intechopen.85147

$$\begin{aligned} \left. \begin{pmatrix} \text{i} \\ \end{pmatrix} \right. \quad \left. \ddot{\boldsymbol{\nu}}\_{I,H} (\lambda\_n, \boldsymbol{\chi}) \right|\_{\boldsymbol{\chi}=0} = 0 \Rightarrow \left. \ddot{\boldsymbol{\nu}}\_{H} (\tilde{\lambda}\_n = \lambda\_n, \boldsymbol{\chi}) \right. = \frac{B(\lambda\_n)}{\sqrt{k(\lambda\_n, \boldsymbol{\chi})}} \sin \left[ \int\_0^\gamma k(\lambda\_n, \boldsymbol{\chi}) d\boldsymbol{\chi} \right], \text{for } 0 < \boldsymbol{\chi} < \lambda\_n \end{aligned} \tag{30}$$

$$\left. \begin{pmatrix} \mathrm{ii} \end{pmatrix} \right. \quad \left. \partial\_{\mathbf{y}} \tilde{\boldsymbol{\mu}}\_{f,\Pi}(\lambda\_n, \mathbf{y}) \right|\_{\boldsymbol{\eta}=\mathbf{0}} = \mathbf{1} / \sqrt{\boldsymbol{\beta}^3} \Rightarrow \mathbf{B}(\lambda\_n) = \mathbf{1} / \sqrt{\lambda\_n \boldsymbol{\beta}^3} \tag{31}$$

connecting to region III in the first quadrant (0 ≤y≤ λn) via the JWKB connection formula.

$$\Psi\_{\boldsymbol{\lambda},\boldsymbol{\Pi}\boldsymbol{\Pi}}(\boldsymbol{\tilde{\lambda}}\_{\boldsymbol{n}}=\boldsymbol{\lambda}\_{\boldsymbol{n}},\boldsymbol{\chi}) = \frac{B(\boldsymbol{\lambda}\_{\boldsymbol{n}})}{\sqrt{\kappa(\boldsymbol{\lambda}\_{\boldsymbol{n}},\boldsymbol{\chi})}} \begin{Bmatrix} \sin\left[a(\boldsymbol{\lambda}\_{\boldsymbol{n}})\right] \exp\left[\boldsymbol{\zeta}(\boldsymbol{\lambda}\_{\boldsymbol{n}},\boldsymbol{\chi})\right] \\ -\frac{1}{2}\cos\left[a(\boldsymbol{\lambda}\_{\boldsymbol{n}})\right] \exp\left[-\boldsymbol{\zeta}(\boldsymbol{\lambda}\_{\boldsymbol{n}},\boldsymbol{\chi})\right] \end{Bmatrix}, \text{for } \boldsymbol{\lambda}\_{\boldsymbol{n}} \leqslant \boldsymbol{\chi} \leqslant \boldsymbol{\infty} \tag{32}$$

whose asymptotic matching gives.

$$\bar{\Psi}\_{l,\text{III}}^{(\text{m.})}(\lambda\_n, y) = -\frac{B(\lambda\_n)}{2\sqrt{\kappa(\lambda\_n, y)}} \cos\left[a(\lambda\_n)\right] \exp\left[-\zeta(\lambda\_n, y)\right], \text{for } \lambda\_n \le y \le \infty \tag{33}$$

ð Þ <sup>ψ</sup> <sup>m</sup>: Again, since we have already obtained <sup>ψ</sup><sup>~</sup> <sup>ð</sup>λn; <sup>y</sup><sup>Þ</sup> and <sup>~</sup> <sup>ð</sup>λn; <sup>y</sup><sup>Þ</sup> in <sup>0</sup><sup>≤</sup> <sup>y</sup>≤λn, J,II J,III JWKB solutions in the second quadrant can be written in terms of them as shown in (19). JWKB wave functions regarding the OP case are given in Figure 3 along with the exact solutions for comparison.

#### 4. The MAF method

If we follow the QHO in dimensionless form given in (7), we have the following properties in MAF theories [3, 18–23]:

#### 4.1 General structure of the MAF approximation to the bound-state wave functions

Formal MAF method suggests a solution to the TISE in (7) in terms of Airy functions as follows:

$$\begin{aligned} \Psi\_{\text{MAF}}(\lambda\_n, y) &= \check{\nu}\_M(\lambda\_n, y) = \begin{cases} F(\lambda\_n, y) Ai[\xi(\lambda\_n, y)] \\ or \\ G(\lambda\_n, y) Bi[\xi(\lambda\_n, y)] \end{cases} \\ \Rightarrow \check{\nu}\_M(\lambda\_n, y) &= a\_1 F(\lambda\_n, y) Ai[\xi(\lambda\_n, y)] + a\_2 G(\lambda\_n, y) Bi[\xi(\lambda\_n, y)] \end{aligned} \tag{34}$$

where Ai and Bi represent the Airy functions (namely, Aið Þ x and Bið Þ x are the ″ linearly independent solutions of the Airy differential equation y xð Þ� xyð Þ¼ x 0 in x), a1&a<sup>2</sup> are the arbitrary constants which will be found from boundary values, and F&G are the functions to be determined. Note that the first variable λ<sup>n</sup> is the eigenenergies (constant values quantized by index n) which will also be determined soon. So, for now, we can consider all these functions as one dimensional in only y for simplification (say, ψ~ ðλn; yÞ ≕ ψ~ ð Þy , ξ λð <sup>n</sup>; yÞ ≕ ξð Þy , Fðλn; yÞ ≕ F yð Þ, etc.). If we <sup>J</sup> <sup>J</sup> choose one of the linearly independent solutions, say F yð Þ:Ai½ξð Þy �, to substitute in the TISE in (4), then it gives:

$$\frac{F'(\mathbf{y})}{F(\mathbf{y})} + \frac{2F'(\mathbf{y})Ai'[\xi(\mathbf{y})]\xi'(\mathbf{y}) + F(\mathbf{y})Ai'[\xi(\mathbf{y})]\xi'(\mathbf{y})}{F(\mathbf{y})Ai[\xi(\mathbf{y})]} + \left\{\xi(\mathbf{y})[\xi'(\mathbf{y})]^2 + f(\mathbf{y})\right\} = \mathbf{0} \quad \text{(35)}$$

Now, with the choice of the last term in (35) as zero, we find the following:

$$\left[\xi(\mathbf{y})[\xi'(\mathbf{y})]^2 + f(\mathbf{y}) = \mathbf{0} \Rightarrow \xi(\mathbf{y}) = \left[\int\_{\mathcal{Y}\_t}^{\mathcal{Y}} \frac{\mathbf{3}}{2} \sqrt{-f(\mathbf{y})} d\mathbf{y}\right]^{2/3} \tag{36}$$

Here, the property of the Airy functions, Ai″ ð Þ¼ ξ ξAið Þξ , was used [3, 18]. The integral interval in (36) is also chosen tactically in a fashion that it invokes a relationship with the turning point yt (representing the correct order yt<sup>1</sup> or yt<sup>2</sup> to give a non-imaginary result), and it can be written in a more explicit and conventional form (by also using in our two-variable form here) as follows:

$$\xi(\boldsymbol{\lambda}\_{n},\boldsymbol{\chi}) = \begin{cases} \xi\_{I} : \left[ \frac{3}{2} \int\_{\mathcal{Y}}^{\eta\_{11}} \kappa(\boldsymbol{\lambda}\_{n},\boldsymbol{\chi}) d\boldsymbol{\mathcal{Y}} \right]^{2/3}, \boldsymbol{\text{for } \boldsymbol{\chi} \le \boldsymbol{\chi}\_{t1} \\\\ \xi\_{II} : -\left[ \frac{3}{2} \int\_{\mathcal{Y}\_{t1}}^{\eta} k(\boldsymbol{\lambda}\_{n},\boldsymbol{\chi}) d\boldsymbol{\mathcal{Y}} \right]^{2/3} = -\left[ \frac{3}{2} \int\_{\mathcal{Y}}^{\eta\_{12}} k(\boldsymbol{\lambda}\_{n},\boldsymbol{\chi}) d\boldsymbol{\chi} \right]^{2/3}, \boldsymbol{\text{for } \boldsymbol{\chi}\_{t1} \le \boldsymbol{\chi} \le \boldsymbol{\chi}\_{t2} \\\\ \xi\_{III} : \left[ \frac{3}{2} \int\_{\mathcal{Y}\_{t2}}^{\eta} \kappa(\boldsymbol{\lambda}\_{n},\boldsymbol{\chi}) d\boldsymbol{\chi} \right]^{2/3}, \boldsymbol{\text{for } \boldsymbol{\chi}\_{t2} \le \boldsymbol{\chi} \end{cases} \tag{37}$$

where <sup>2</sup> <sup>f</sup>ðλn; <sup>y</sup>≥0Þ ¼ <sup>k</sup><sup>2</sup> ðλn; y޼�κ ðλn; yÞ and yt1&yt<sup>2</sup> are the CTPs at the interface of the regions I � II&II � III, respectively. The remaining terms in (44) and (45) are also made zero as follows:

Starting from the second term, we have.

$$\frac{2F'Ai'(\xi)\xi'}{F} + Ai'(\xi)\xi' = 0 \Rightarrow F(y) = \frac{b\_1}{\sqrt{\xi'(y)}}\tag{38}$$

where b<sup>1</sup> is some constant, and finally, making the first term in (35) zero (which is the only assumption in the MAF method), we have the following:

$$P(\boldsymbol{y}) = \frac{\boldsymbol{F}^{\boldsymbol{\prime}}(\boldsymbol{y})}{\boldsymbol{F}(\boldsymbol{y})} \approx \mathbf{0} \tag{39}$$

Or more correctly in two-variable form in our eigenvalue system.

$$P(\lambda\_n, \mathbf{y}) = \frac{\partial^2\_{\mathbf{y}} F(\lambda\_n, \mathbf{y})}{F(\lambda\_n, \mathbf{y})} \tag{40}$$

can be thought as a measure of the accuracy of the MAF solution, namely, Pðλn; yÞ ! 0, as MAF solution gets more accurate [18].

The same results would also be obtained if we had chosen the other linearly independent solution, G yð Þ:Bið Þy , in (34). Consequently, using the results found in (37) and (38), the general solution suggested in (34) can be written explicitly in the standard form of the MAF formula as follows:

$$\psi\_{MAF}(\mathbf{y}) = \frac{\tilde{c}\_1}{\sqrt{\xi'(\mathbf{y})}} Ai[\xi(\mathbf{y})] + \frac{\tilde{c}\_2}{\sqrt{\xi'(\mathbf{y})}} Bi[\xi(\mathbf{y})] \tag{41}$$

or more correctly in two variables here in our study.

$$\Psi\_{\rm MAP}(\lambda\_n, y) = \frac{\tilde{c}\_1}{\sqrt{\partial\_{\mathcal{Y}} \xi(\lambda\_n, y)}} \text{Ai}[\xi(\lambda\_n, y)] + \frac{\tilde{c}\_2}{\sqrt{\partial\_{\mathcal{Y}} \xi(\lambda\_n, y)}} \text{Bi}[\xi(\lambda\_n, y)] \tag{42}$$

where c<sup>1</sup> ¼ a1:b<sup>1</sup> and c<sup>2</sup> ¼ a2:b<sup>2</sup> are the arbitrary constants to be determined from the boundary values as mentioned and <sup>∂</sup>yξ λ<sup>ð</sup> ; <sup>y</sup><sup>Þ</sup> represents the first derivative of ξ with respect to y. Using the result in (38), the approximation term Pðλn; yÞ in (40) can be rewritten explicitly as follows:

$$P(\lambda\_n, \boldsymbol{y}) = \frac{\mathfrak{Z}}{4} \left[ \frac{\partial\_{\boldsymbol{y}}^2 \xi(\lambda\_n, \boldsymbol{y})}{\partial\_{\boldsymbol{\mathcal{Y}}} \xi(\lambda\_n, \boldsymbol{y})} \right]^2 - \frac{\partial\_{\boldsymbol{y}}^3 \xi(\lambda\_n, \boldsymbol{y})}{2 \partial\_{\boldsymbol{\mathcal{Y}}} \xi(\lambda\_n, \boldsymbol{y})} \tag{43}$$

#### 4.2 MAF solution of eigenenergies

For a symmetrical f as in Figure 2, we have even-parity (EP) and odd-parity (OP) MAF wave functions just as in JWKB method, but now it leads to two different MAF quantization formulas with two different MAF universal constants regarding EP and OP solutions as given in [3] and as we study in this section. We again use the symbolism in (9) (φ⇔ψ) and start with the first quadrant, by applying that lim<sup>y</sup>!<sup>∞</sup> ¼ 0 requires c<sup>2</sup> ¼ 0 in (42), namely,

$$\Psi\_{\rm MAP,n}(\lambda\_n, \boldsymbol{\chi}) = :\tilde{\boldsymbol{\mu}}\_{\rm Mn}(\lambda\_n, \boldsymbol{\chi}) = \frac{\tilde{\boldsymbol{c}}\_1}{\sqrt{\partial\_{\boldsymbol{\chi}} \xi(\lambda\_n, \boldsymbol{\chi})}} Ai[\xi(\lambda\_n, \boldsymbol{\chi})],\tag{44}$$

where the denominator can be written in the following form [3]:

$$\frac{1}{|\sqrt{\partial\_{\mathcal{Y}}\xi(\lambda\_n,\mathcal{Y})}|} = \frac{|\xi(\lambda\_n,\mathcal{Y})|^{1/4}}{\left|k^2(\lambda\_n,\mathcal{Y})\right|^{1/4}}\tag{45}$$

i. Even-parity (EP) eigenenergies: if we apply the EP formulas of the exact solution in (17), by using (16), to the MAF wave functions, we have the following:

$$\left\{\ddot{\boldsymbol{\nu}}\_{\mathrm{Mn}}(\boldsymbol{\lambda}\_{n},\mathbf{0}) = \frac{\pm p}{\sqrt{\beta}}, \partial \frac{\ddot{\boldsymbol{\nu}}\_{\mathrm{Mn}}(\boldsymbol{\lambda}\_{n},\mathbf{y})}{\partial \boldsymbol{\nu}}\Big|\_{\boldsymbol{\nu}=\mathbf{0}} = \mathbf{0}\right\}\tag{46}$$

<sup>1</sup> ð Þ<sup>i</sup> <sup>ψ</sup><sup>~</sup> <sup>M</sup>ðλn; <sup>0</sup>Þ ¼ <sup>p</sup>ffiffiffi ) find <sup>~</sup>c<sup>1</sup> <sup>≕</sup>~c1<sup>s</sup> <sup>β</sup> sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>0</sup> <sup>0</sup> (47) <sup>ξ</sup> <sup>ξ</sup> Ai½ξðλn; <sup>y</sup>Þ� <sup>0</sup> <sup>0</sup> result :) <sup>~</sup>c1<sup>s</sup> <sup>¼</sup> <sup>p</sup>ffiffiffi ) <sup>ψ</sup><sup>~</sup> <sup>M</sup>ðλn; <sup>y</sup>Þ ¼ <sup>β</sup>Ai½ξ λ<sup>ð</sup> <sup>n</sup>; <sup>y</sup>Þ� <sup>β</sup>∂yAi½ ðξ λn; <sup>y</sup>Þ� Aið Þ <sup>ξ</sup><sup>0</sup>

$$\text{H}(\text{ii}) \quad \partial \frac{\tilde{\boldsymbol{\mu}}\_{M}(\boldsymbol{\lambda}\_{n}, \boldsymbol{\mathcal{y}})}{\partial \boldsymbol{\mathcal{y}}} \Big|\_{\boldsymbol{\mathcal{y}}=\boldsymbol{0}} = \mathbf{0} \Rightarrow \text{Ai}'(\boldsymbol{\xi}\_{0}) - \frac{\text{Ai}(\boldsymbol{\xi}\_{0})\boldsymbol{\xi}\_{0}^{\cdot}}{2\boldsymbol{\xi}\_{0}^{\cdot 2}} = \mathbf{0} \tag{48}$$

� � � 0 <sup>0</sup> for simplification, namely, ξ <sup>0</sup> <sup>¼</sup> <sup>∂</sup>yyð<sup>λ</sup> <sup>Þ</sup> <sup>0</sup> ″ where we � used ξ0, ξ0, and ξ ″ <sup>0</sup> ¼ ξ λð <sup>n</sup>; 0Þ, , and <sup>ξ</sup> <sup>y</sup>¼<sup>0</sup> <sup>ξ</sup> , respectively. We assumed here <sup>0</sup> <sup>¼</sup> <sup>∂</sup>yðλn; <sup>y</sup><sup>Þ</sup> <sup>n</sup>; <sup>y</sup> <sup>y</sup>¼<sup>0</sup> φ λð <sup>n</sup>; 0Þ ¼ 1, and we will use then parity correction for φ λð <sup>n</sup>; 0޼�1 case just as in the JWKB calculations. If we take the derivative of (36), we get.

$$\left[\xi^{\prime 3}(\lambda\_n, y) + 2\xi(\lambda\_n, y)\right] \not{q}\_{\mathcal{V}}^{\prime} \xi(\lambda\_n, y) \left[\not{q}\_{\mathcal{V}}^{\prime} \xi(\lambda\_n, y)\right] \not{q}\_{\mathcal{V}} \not{q}\_{\mathcal{V}}^2(\lambda\_n, y) = 0 \tag{49}$$

where the last term vanishes as y ! 0 to give.

$$
\xi^{\prime 2}\_{\ \ 0} + 2\xi\_0 \xi\_0^{\prime} = 0 \Rightarrow \frac{\xi\_0^{\prime}}{2\xi^{\prime 2}\_{\ \ 0}} = -\frac{1}{4\xi\_0} \tag{50}
$$

whose substitution in (48) gives.

$$
\xi\_0 Ai'(\xi\_0) + \frac{1}{4} Ai(\xi\_0) = 0 \tag{51}
$$

Now, by the substitution of ξ<sup>0</sup> ! �Zsn, we have.

$$-Z\_{sn} \text{Ai}'(-Z\text{sn}) + \frac{\text{Ai}(-Z\_{sn})}{\mathbf{4}} = \mathbf{0} \tag{52}$$

where the subscripts sn stand for s, symmetrical solution (EP), and n, quantization order (nth quantization), and Zsn is the nth solution of the differential equation in (52) regarding the symmetrical solution. Now, by using the results in (13), we find the MAF quantization formula regarding the symmetrical solution:

$$\iint\limits\_{\mathbf{0}} \mathbf{k} \left[ \bar{\lambda}\_{M,m}, \mathbf{y} \right] d\mathbf{y} =: \left( \zeta\_m + \frac{1}{2} \right) \boldsymbol{\pi} \Rightarrow \zeta\_m = \quad \frac{4Z\_m^{3/2}}{3\pi} - \left( \frac{1}{2} \right), n = 1, 2, 3, \dots \tag{53}$$

where ζsn is the universal MAF constants regarding the symmetrical solution whose values are given in Table 2 along with the JWKB solutions (which are already exact) for some n values in comparison. Note that we used ns to represent the symmetrical (EP) MAF indices in Table 2.

ii. Odd-parity (OP) eigenenergies: similarly, if we apply the OP formulas of the exact solution in (17), by using (16), to the MAF wave functions, we have the following:

$$\left\{ \ddot{\boldsymbol{\nu}}\_{M}(\boldsymbol{\lambda}, \mathbf{0}) = \mathbf{0}, \left. \partial\_{\boldsymbol{\mathcal{V}}} \ddot{\boldsymbol{\mu}}\_{M}(\boldsymbol{\lambda}, \mathbf{y}) \right|\_{\boldsymbol{\mathcal{y}} = \mathbf{0}} = \frac{\pm q}{\sqrt{\boldsymbol{\beta}^{3}}} \right\} \left( \mathbf{0} \tag{54}$$

$$\mathbf{i}(\mathbf{i}) \quad \ddot{\boldsymbol{\mu}}\_{M}(\boldsymbol{\lambda}, \mathbf{0}) = \mathbf{0} \Rightarrow \mathbf{A} \mathbf{i}(\boldsymbol{\xi}\_{0} \rightarrow -\mathbf{Z}\_{an}) = \mathbf{A} \mathbf{i}(-\mathbf{Z}\_{an}) = \mathbf{0}, n = \mathbf{1}, 2, 3, \dots \tag{55}$$

Quantum Harmonic Oscillator DOI: http://dx.doi.org/10.5772/intechopen.85147

$$\left. \begin{pmatrix} \mathrm{ii} \end{pmatrix} \right. \quad \partial\_{\mathcal{T}} \tilde{\boldsymbol{\mu}}\_{M} (\lambda\_{\boldsymbol{n}\_{a}}, \boldsymbol{y}) \big|\_{\boldsymbol{y} = 0} = \frac{1}{\sqrt{\boldsymbol{\beta}^{3}}} \underset{\boldsymbol{\beta}}{\underset{\boldsymbol{\beta}}{\operatorname{\boldsymbol{\beta}}}} \operatorname{Ai}'(\xi\_{0}) - \frac{\operatorname{Ai}(\xi\_{0}) \tilde{\boldsymbol{\xi}}\_{0}^{\prime}}{2 \boldsymbol{\xi}^{\prime 2}\_{\boldsymbol{0}}} = \frac{1}{\sqrt{\boldsymbol{\beta}^{3}}} \underset{\boldsymbol{\beta}}{\underset{\boldsymbol{\beta}}{\operatorname{\boldsymbol{\beta}}}} \operatorname{find} \vec{c}\_{1} =: \bar{c}\_{1a} \quad \text{(56)}$$

where, similarly, the subscripts an stand for a, antisymmetrical (OP), and n, quantization order (nth quantization), and Zan is the nth solution of the equation in (55) regarding the asymmetrical solutions. Similarly, by using the results in (13), we find the MAF quantization formula regarding the antisymmetrical solution:

$$\int\_{0}^{\gamma\_{1}} k\left[\hat{\lambda}\_{M,an},\gamma\right] d\circ =: \left(\zeta\_{an} + \frac{1}{2}\right)\pi \Rightarrow \zeta\_{an} = \quad \frac{4Z\_{an}^{3/2}}{3\pi} - \left(\frac{1}{2}\right), n = 1, 2, 3, \dots \tag{57}$$

where ζan is the universal MAF constants regarding the antisymmetrical solution whose values are given in Table 2 along with the JWKB solutions (which are already exact) for some n values in comparison. Note that we used na to represent the antisymmetrical (OP) MAF indices in Table 2.

#### 4.3 MAF solution of eigenfunctions

By using a tentative boundary condition with q = 1 for the EP solutions, we have found the result in (47), and we said that we would extend it by considering the parity matching for q ¼ �q. Consequently,

$$\tilde{\boldsymbol{\psi}}\_{(M,E,\boldsymbol{P}\_{\cdot})}^{(\text{par.m.})}(\boldsymbol{\lambda}\_{n},\boldsymbol{\chi}\geq\mathbf{0}) = (-1)^{\binom{n\_{\rm 2}}{2}} \times \frac{\tilde{\boldsymbol{c}}\_{\text{3}}}{\sqrt{\partial\_{\boldsymbol{\mathcal{V}}}\xi(\boldsymbol{\lambda}\_{n},\boldsymbol{\chi})}} \boldsymbol{A} \boldsymbol{i}[\boldsymbol{\xi}(\boldsymbol{\lambda}\_{n},\boldsymbol{\chi})] \; \boldsymbol{n}\_{\boldsymbol{s}} = \boldsymbol{0}, \boldsymbol{2, 4, ...} (\text{even}) \quad \text{(58)}$$

$$\hat{\boldsymbol{\mu}}\_{(M,EP.)}^{(\text{par.m.})}(\lambda\_{n\_{\text{\tiny\text{-}}}},\boldsymbol{\chi}\geq\mathbf{0}) = (-1)^{\binom{n\_{\text{\tiny\text{-}}}{2}}{2}} \times \sqrt{\frac{\xi\_{0}^{\prime}}{\beta \partial\_{\boldsymbol{\beta}} A i[\boldsymbol{\xi}(\lambda\_{n\_{\text{\tiny\text{-}}}},\boldsymbol{\chi})]}} \begin{At} i[\boldsymbol{\xi}(\lambda\_{n\_{\text{\tiny\text{-}}}},\boldsymbol{\chi})]}{Ai(\xi\_{0})}, n\_{\text{\tiny\text{-}}} = 0, 2, 4, ... (even) \tag{59}$$

Similarly, for the antisymmetric parity wave functions, we have.

$$\tilde{\psi}\_{(M,O,P,)}^{(\text{par.m.})}(\lambda\_{u\_4}, y \ge 0) = (-1)^{\binom{n\_4 - 1}{2}} \times \frac{\tilde{c}\_{1d}}{\sqrt{\partial\_Y \xi(\lambda\_u, y)}} \text{Ai}[\xi(\lambda\_u, y)], n\_i = \mathbf{1}, \mathbf{3}, \mathbf{5}, ...(odd) \tag{60}$$

where constant coefficients ~c1<sup>s</sup> and ~c1<sup>a</sup> represent the related symmetric and antisymmetric coefficients, respectively.

#### 5. MAF solution of the QHO

Again, we use the schematic sketch given in Figure 2 for the QHO under study.

#### 5.1 MAF eigenenergies of the QHO

Since we have tactically used (53)–(57) to resemble the MAF quantization formula to the JWKB quantization formula given in (12), by using the result of

calculation of the same integral in (13), we have the following results regarding the MAF eigenenergies of the QHO:

$$E.P.: \bar{\lambda}\_{M,m}(E\_m) = \sqrt{2\zeta\_m + 1} \text{(for)} \\ \bar{E}\_{M,m} = E\_{MAF,m} = \left(\zeta\_m + \frac{1}{2}\right) \left(m = 1, 2, 3, \dots\right) \tag{61}$$

$$O.P.: \tilde{\lambda}\_{M,an}(E\_{an}) = \sqrt{2\zeta\_{an} + 1} \\
\text{(for)} \\
\tilde{E}\_{\overline{M,an}} = E\_{\text{MAP},an} = \because \left(\zeta\_n + \frac{1}{2}\right)\zeta; n = 1, 2, 3, \dots \Bigg\} \\
\text{(62)}$$

MAF eigenenergies are given in Table 2 along with the JWKB solutions (which are already exact) for some n values in comparison. (Note again that we used ns and na to represent the symmetrical (EP) and antisymmetrical MAF indices in Table 2, respectively).

#### 5.2 MAF eigenfunctions of the QHO

For the regions IIb and III, we have the following definitions:

$$f(\lambda\_n, y) = k^2(\lambda\_{\overline{n}}, y) = \begin{cases} \not\succ \left(\bar{\lambda}\_{\mathsf{M}} \quad \bar{\lambda}\_n, y\right) \not\preceq \bar{\lambda}\_n^2 - y^2; 0 \le y \le \bar{\lambda}\_n\\ \not\succ \left(\bar{\lambda}\_{\mathsf{M}} \quad \bar{\lambda}\_n, y\right) \not\preceq y^2 - \bar{\lambda}\_n^2; \bar{\lambda}\_n \le y < \infty \end{cases} \tag{63}$$

Calculation of ξ in (37) for the first quadrant gives.

 � � � � ~ λ<sup>M</sup> ~ ~ 8 � <sup>R</sup> � � �<sup>2</sup>=<sup>3</sup> <sup>h</sup> <sup>R</sup> � � <sup>i</sup><sup>2</sup>=<sup>3</sup> <sup>ξ</sup>IIb : � <sup>3</sup> <sup>λ</sup>M; <sup>y</sup> <sup>λ</sup>M; <sup>y</sup> >> <sup>y</sup> <sup>y</sup> ¼� <sup>3</sup> <sup>t</sup><sup>2</sup> <sup>2</sup> < k dy k dy , for 0≤ y≤y <sup>2</sup> <sup>0</sup> <sup>y</sup> <sup>t</sup><sup>2</sup> ξðλn; yÞ¼ > h <sup>R</sup> � � <sup>i</sup> <sup>&</sup>gt; 2=3 : <sup>3</sup> <sup>y</sup> dy , for <sup>y</sup> <sup>ξ</sup>III : <sup>t</sup><sup>2</sup> <sup>≤</sup><sup>y</sup> ~ ~ ~ ~ ~λM; y ( " !#)<sup>2</sup>=<sup>3</sup> 3<sup>2</sup>=<sup>3</sup> � � <sup>2</sup> <sup>y</sup> <sup>λ</sup>M; <sup>y</sup> <sup>λ</sup> λM; y λ<sup>M</sup> ≤y , ∞ <sup>κ</sup> <sup>2</sup> yt<sup>2</sup> 8 >> >>>> <sup>ξ</sup>IIb : � �2yk <sup>þ</sup> <sup>π</sup> � 2arctan , for <sup>0</sup>≤<sup>y</sup> <sup>≤</sup> >>< <sup>M</sup> 4 k <sup>~</sup> <sup>~</sup> <sup>~</sup> <sup>~</sup> <sup>¼</sup>( " #)<sup>2</sup>=<sup>3</sup> <sup>λ</sup><sup>M</sup> <sup>λ</sup>M; <sup>y</sup> <sup>λ</sup> � � λM;<sup>y</sup> <sup>&</sup>gt;<sup>2</sup> <sup>3</sup> <sup>M</sup>ln >>> <sup>2</sup>=<sup>3</sup> >> yκ þ <sup>&</sup>gt; <sup>y</sup><sup>κ</sup> >: <sup>ξ</sup> , for III : <sup>24</sup><sup>=</sup><sup>3</sup> (64)

<sup>~</sup><sup>λ</sup> <sup>~</sup> <sup>M</sup> (and also EM below) are MAF eigenenergies which become (61) for the symmetric (EP) and (62) for the antisymmetric (OP) case accordingly. Calculation of the constant coefficient in (58) or (59) for the symmetric boundary values given in (46) with q = 1 gives. where

$$E.P.: \bar{c}\_{1\*} = \frac{\sqrt{\bar{2}} \bar{E}\_{s\_n}^{1/6}}{\sqrt{\beta} (3\pi)^{1/6} \text{Ai}\left[\sqrt{\frac{1}{4} \bar{E}\_{s\_n}^{4/3} (3\pi)^{2/3}}\right]} \tag{65}$$

Figure 4. Exact and JWKB solutions of OP wave functions (for q ¼ 1).

Figure 5. Relative and absolute error of EP MAF solutions.

Quantum Harmonic Oscillator DOI: http://dx.doi.org/10.5772/intechopen.85147

Similarly, calculation of the constant coefficient in (60) for the antisymmetric boundary values given in (54) with q = 1 gives.

$$O.P.: \tilde{c}\_{1a} = \frac{\tilde{E}\_{a\_u}^{7/6} (3\pi)^{1/6}}{\sqrt{2\rho^3} \left\{ \bar{E}\_{a\_u}^{4/3} (3\pi)^{2/3} A \text{i}^{\prime} \left[ -\frac{1}{4} \bar{E}\_{a\_u}^{4/3} (3\pi)^{2/3} \right] - A \text{i} \left[ -\frac{1}{4} \bar{E}\_{a\_u}^{4/3} (3\pi)^{2/3} \right] \right\}} \tag{66}$$

Since the MAF solutions of both EP and OP solutions are very close to the exact solutions given in Figures 3 and 4, their absolute and relative error graphs with respect to the exact solution are given in Figures 5 and 6. We can also see that there are no discontinuities at the CTPs in the MAF solutions when compared with the JWKB solutions given in Figures 3 and 4.

#### 6. Conclusion

Here we studied the fundamental outcomes of the two conventional semiclassical approximation methods, namely, JWKB and MAF methods pedagogically, and obtained the solutions of the QHO by these semiclassical methods by using the parity conditions of the expected solutions by using the dimensionless form of the QHO system. We applied the asymptotic matching and parity matching procedure to obtain the correct form of semiclassical solutions. As expected, JWKB solutions diverge at and around the CTPs, whereas MAF solutions do not. As also expected (since being typical), JWKB eigenenergies are exact, whereas MAF eigenenergies are unfortunately not but very accurate as expected from an approximation method. In the MAF method, function p in (40) or in (43) is assumed zero. Indeed, it is very close to zero to give approximate results, and function P in (40) or in (43) can be used as an approximation criterion for the MAF method [3, 18]. However, improved MAF methods (IMAF) or perturbation corrections concerning the nonzero P function seem straightforward to improve the accuracy of the MAF solutions as in [3, 20, 22]. Normally, for an even potential function in the TISE, EP and OP initial values are as given in (17), but due to the conversion factor β in (11) or (16), for the QHO in the dimensionless form (in ψ), we have (22) and (29). In our notation, we have used the notation, φ⇔ψ, where real physical system is in φ and the dimensionless form is in ψ. Since the standard formulation is given according to the real physical systems, JWKB and MAF formulas in the literature such as in [1–7, 19–23] surely correspond to the initial values β ! 1 in our dimensionless form formulation in ψ. Consequently, we hereby present a full JWKB and MAF solution concerning the quantized conversion factor β in (11).

#### Acknowledgements

Author acknowledges special thanks to the IntechOpen for the financial support in the publishment of this chapter.

Quantum Harmonic Oscillator DOI: http://dx.doi.org/10.5772/intechopen.85147

### Author details

Coşkun Deniz Department of Electrical and Electronics Engineering, Faculty of Engineering, Adnan Menderes University, Turkey

\*Address all correspondence to: cdeniz@adu.edu.tr

© 2019 The Author(s). Licensee IntechOpen. This chapteris distributed underthe terms oftheCreative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

### References

[1] Liboff R. Introductory Quantum Mechanics, 4/E. Cornell University: Addison Wesley; 2001

[2] Griffiths DJ. 20 Introduction to Quantum Mechanics, 2nd ed. Upper Saddle River, NJ: Pearson; 2005. ISBN: 0-13-191175-9

[3] Ghatak AK, Gallawa RL, Goyal IC. Modified Airy Functions and WKB Solutions to the Wave Equation. Washington: NIST; 1991

[4] Deniz C. Semiclassical anomalies of the quantum mechanical systems and their modifications for the asymptotic matching. Annals of Physics. 2011; 326(8):1816-1838

[5] Landau LD, Lifshitz EM. Quantum mechanics, non-relativistic theory. In: Course of Theoretical Physics. 2nd ed. Vol. 3. NY: Pergamon; 1965

[6] Hruska M, Keung WY, Sukhatme U. Accuracy of semiclassical methods for shape-invariant potentials. Physical Review A. 1997;55(5):3345-3350

[7] Bender CM, Orszag SA. Advanced Mathematical Methods for Scientists and Engineers Asymptotic Methods and Perturbation Theory. NY: Springer– Verlag; 1999

[8] Aydin MC, Uncu H, Deniz C. A parabolic model for dimple potentials. Physica Scripta. 2013;88(3):035006. Available from: http://iopscience.iop. org/article/10.1088/0031-8949/88/03/ 035006/meta

[9] Gómez-Vergel D, Villaseñor E. The time-dependent quantum harmonic oscillator revisited: Applications to quantum field theory. Annals of Physics. 2009;324:1360-1385

[10] Tong T. Quantum Field Theory, University of Cambridge Part III.

Mathematical Tripos. 2007. Available from: http://www.damtp.cam.ac.uk/ user/tong/qft.html

[11] Gin-Ge Chen B, Derbes D, Griffiths D, Hill B, Sohn R, Ting YS. Lectures of Sidney Coleman on Quantum Field Theory. 1st ed. NJ: World Scientific; 2018

[12] Shankar R. Principles of Quantum Mechanics. Chapter 10, Exercises 10.2.2–10.2.3. 2nd ed. NY: Springer, Plenum Press; 1994. pp. 259-260. ISBN 0-306-44790-8

[13] Shankar R. Principles of Quantum Mechanics. Chapter 12. 2nd ed. NY: Springer, Plenum Press; 1994. pp. 351-352. ISBN 0-306-44790-8

[14] Deniz C. On the exact and JWKB solution of 1D quantum harmonic oscillator by mathematica. Journal of Physics: Conference Series. 2016; 707(012033):1-11. Available from: h ttp://iopscience.iop.org/article/10.1088/ 1742-6596/707/1/012033/pdf

[15] Lindblad G. Brownian motion of a quantum harmonic oscillator. Reports on Mathematical Physics. 1976;10(3): 393-406

[16] Ford GW, Lewis JT, O'Connel RF. Quantum Langevin equation. Physical Review A. 1988;37(11):4419-4428

[17] Stenholm S. Quantum theory of linear friction. Brazillian Journal of Physics. 1997;27(2):214-237

[18] Deniz C, Gerceklilu M. An analysis of the exactly solvable linear differential equations by the modified airy function (MAF) method. Indian Journal of Physics. 2011;85(2):339-357. DOI: 10.1007/s12648-011-0010-1. Available from: http://link.springer.com/article/ 10.1007/s12648-011-0010-1

Quantum Harmonic Oscillator DOI: http://dx.doi.org/10.5772/intechopen.85147

[19] Goyal IC, Gallawa RL, Ghatak AK. Improved variational analysis of inhomogenous optical waveguides using airy functions. Journal of Lightwave Technology. 1993;11(10):1575-1578

[20] Goyal IC, Rajeev J, Ghatak AJ. Planar optical waveguides with arbitrary index profile: An accurate method of analysis. Journal of Lightwave Technology. 1997;15(11):2179-2182

[21] Ghatak AK, Goyal IC, Jindal R, Varshni YP. MAF solution for bounded potential problems. Canadian Journal of Physics. 1998;76:351-359

[22] Goyal IC. Modified airy function solutions to optical waveguide problems. In: Proceedings of the 2002 4th International Conference on Transparent Optical Networks; IEE. 1. 2002. pp. 155-1605. Also available by print ISBN: 0-7803-7375-8/02

[23] Deniz C. On the MAF solution of the uniformly lengthening pendulum via change of independent variable in the Bessels equation. Results in Physics. 2017;733:333-343. Available from: h ttps://www.sciencedirect.com/science/a rticle/pii/S2211379717314092

[24] Boccara N. Essentials of Mathematica with Applications to Mathematics and Physics. NY: Springer; 2007. pp. 475-480

[25] Langer RE. On the asymptotic solutions of ordinary differential equations, with an application to the bessel functions of large order. The Transactions of the American Mathematical Society. 1931;33:23-64

Chapter 3

Systems

Ozkan Ozturk

Abstract

1. Introduction

47

Oscillation Criteria of

contributes several theoretical results to the literature.

nonlinear system, fixed point theorems

and Peterson in 2001 and 2003, respectively.

Two-Dimensional Time-Scale

Oscillation and nonoscillation theories have recently gotten too much attention and play a very important role in the theory of time-scale systems to have enough information about the long-time behavior of nonlinear systems. Some applications of such systems in discrete and continuous cases arise in control and stability theories for the unmanned aerial and ground vehicles (UAVs and UGVs). We deal with a two-dimensional nonlinear system to investigate the oscillatory behaviors of solutions. This helps us understand the limiting behavior of such solutions and

Keywords: oscillation, nonoscillation, two-dimensional systems, time scale,

This chapter analyses the oscillatory behavior of solutions of two-dimensional (2D) nonlinear time-scale systems of first-order dynamic equations. We also investigate the existence and asymptotic properties of such solutions. The tools that we use are the most well-known fixed point theorems to consider the sign of the component functions of solutions of our system. A time scale, denoted by T, is an arbitrary nonempty closed subset of the real numbers R, which is introduced by a German mathematician, Stefan Hilger, in his PhD thesis in 1988 [1]. His primary purpose was to unify continuous and discrete analysis and extend the results to one comprehensive theory. For example, the results hold for differential equations when T ¼ R, while the results hold for difference equations when T ¼ Z. Therefore, there might happen to be two different proofs and maybe similar in most cases. In other words, our essential desire is to combine continuous and discrete cases in one comprehensive theory and remove the obscurity from both. For more details in the theory of differential and difference equations, we refer the books [2–4] to interested readers. As for the time-scale theory, we assume most of the readers are not familiar with the time-scale calculus, and thus we give a concise introduction to the theory of time scales from the books [5, 6] written by Bohner

Two-dimensional dynamical systems have recently gotten too much attention because of their potential in applications in engineering, biology, and physics (see, e.g., [7–11]). For example, Bartolini and Pvdvnowski [12] consider a nonlinear

#### Chapter 3
