On Aspects of Continuous Approximation of Diatomic Lattice

Abstract

This paper is devoted to the continualization of a diatomic lattice, taking into account natural intervals of wavenumber changes. Continualization refers to the replacement of the original pseudo-differential equations by a system of PDEs that provides a good approximation of the dispersion relations. In this regard, the Padé approximants based on the conditions for matching the values of the dispersion relations of the discrete and continuous models at several characteristic points are utilized. As a result, a sixth-order unconditionally stable system with modified inertia is obtained. Appropriate boundary conditions are formulated. The obtained continuous approximation accurately describes the amplitude ratios of neighboring masses. It is also shown that the resulting continuous system provides a good approximation for the natural frequencies.

1. Introduction

A diatomic lattice was proposed as a useful model in the theory of heat capacity by Born and von Kármán [1,2]. This model is currently used to describe various metamaterials [3,4,5,6]. Therefore, the dynamics of the diatomic lattice has become a subject of extensive investigations, in which both discrete model and its continuous approximations have been considered [3,5,7,8,9,10,11,12].

It is worth mentioning the relationship between the discrete and continuous models. This issue was first discussed by the ancient Greeks. After Newton’s Calculus, the continuous models had been a priority for a long time, and continuous mathematics accumulated the huge number of tools for their analysis. The situation began to change with the advent of modern computers. Kolmogorov wrote about the relations between discrete and continuous mathematics [13]: “It is very likely that with the development of modern computational techniques it will be understood that in many cases it makes sense to study real phenomena without using the intermediate step of their stylization in the form of infinite and continuous mathematics and proceeding directly to discrete models. Pure mathematics was successfully developed mainly as a science of the infinite. Obviously, this state of affairs is deeply rooted in our consciousness, which works with great ease with an intuitively clear understanding of unbounded sequences, limiting processes, continuous and even smooth manifolds, etc.”.

Continuous models can be treated as an asymptotic approximation to discrete models. Their strength is that, as Kolmogorov noted, they are most accessible to our understanding (or perhaps our understanding is distorted by traditional mathematical education, which tends to bias thinking in the direction of continuous mathematics). The advantage of continuous models is the possibility of using the powerful arsenal of continuous mathematics for their analysis. At the same time, when using the continuous models to describe a discrete structure, one must consider the natural limitations of these approximations to avoid various “paradoxes” and artifacts [14,15,16,17].

Interestingly, continualization can lead to important insights, even in cases where it does not allow a complete understanding of a discrete model. For example, “When Kruskal derived the Korteweg–de Vries (K-dV) equation from the FPU lattice, he bypassed the blockade imposed by the discrete to reach the enlightenment of the continuum. In this moment of the journey of the FPU problem, the sweet fruit of solitonic integrability, one of the most beautiful scientific discoveries of the second half of the last century, almost made us forget that bypassing the blockade did not eliminate it. As later studies have clearly shown, crossing the weak excitation threshold causes discretion to re-emerge in full force, or in all its glory if you will, displacing K-dV and its siblings. The maturity of fifty years of experience now makes us realize that Kruskal’s seminal leap is not a license to ignore the discrete system, but rather a bridge to a synergetic methodology wherein we mediate between this true system, for which our analytical means are quite poor, and an idealized continuum system about which we can say many things.” [18].

The connection between discrete and continuous approximations is established by the methods of discretization and continualization (homogenization) [3,4,8,14,15,19,20,21,22]. Discrete and continuous models are fundamentally different from each other. Thus, when we replace one model by another, we can usually only achieve a good approximation of some properties of the original object. Continualization often focuses on the comparison of dispersion relations [3,14,15,22]. In addition, one can assume the correspondence of the asymptotic behavior of wave processes in discrete and continuous media [19], the equality of the group impedance matrices of the model and the modeled structure [4]. Homogenization of the discrete structure can also mean that a periodic structure is replaced by a periodic structure of the same dimension and geometry whose periodicity cell is the simplest discrete mechanical oscillatory system with the smallest possible number of discrete parameters [3,4].

In the paper [20], the possibilities of deterministic discretization and chaotic continualization were shown using the example of the logistic equation. This confirms the ambiguity of the discretization and continualization processes.

Continualization is commonly understood as the approximation of non-local (pseudo-differential or integral) operators by local (differential) ones. The simplest possibility is a polynomial approximation [16] (the term “standard continuation” is also used for this approximation [22]). It can be based on the Taylor series expansion of non-local operators or the analysis of dispersion relations [16,17]. A continuous approximation based on the first roots of the dispersion equation is also possible [16,17]. These procedures are only suitable for describing the lower part of the spectrum or long waves and can lead to ill-posed problems. In this regard, methods of regularization and extension for the range of applicability of continuum models have become widespread. Padé approximants (non-standard continualization [22]) [5,8,14,15,18,19,20,21] are often used for this purpose. They allow the construction of improved long-wave approximations that approximate well the asymptotic behavior of discrete wave processes or describe frequencies in the first Brillouin zone. However, it must be taken into account that the use of such approximations has significant limitations. As mentioned in [16], p. 65, “other approximation models that are not associated with the polynomial approximation are also possible. An approximate description of the dispersion curve over a wide range of wavelengths using an appropriate function of the wavenumber k (for example when interpolating experimental data) can be of interest in a number of cases. This model is good for describing non-decaying waves, but in general it is incorrect to continue it in the complex region and use it for boundary problems”.

Our paper is devoted to the linear diatomic lattice. This classical problem still arouses interest despite its long history [23,24,25,26,27]. In particular, [23,24,25] analyzed the question of the range of wavenumbers that allow a correct description of the wave processes in the diatomic lattice. The peculiarities of diatomic lattice modeling mentioned in [24,25] are considered in our work.

There is extensive literature dealing with diatomic lattices. We only refer to studies related to continualization using Padé approximants [8,21] and the construction of higher-order continuous models [6].

Our paper is concerned with the continuization of diatomic lattices based on the multipoint Padé approximants (Padé approximants of the second kind [28,29]). The proximity of the dispersion curves of discrete and continuous models is used as a criterion for approximation. It is shown that each natural frequency is also well approximated in this case. An important feature of the continuous approximation constructed in this work is that it is considered in the most natural periodic interval of wavenumbers.

2. Discrete System

Consider a finite lattice consisting of 2N + 2 alternating masses $m_1$ and $m_2$ connected by springs of stiffness c (see Figure 1). We construct the final chain so that the number of particles with masses $m_1$ and $m_2$ is equal. Since the total number of particles is assumed to be large, this restriction is not important.

The oscillation of this system is described by a system of mixed differential-difference equations:

$$\begin{array}{l}{m}_1{\ddot{u}}_n+c\left(2u_n-u_{n+1}-u_{n-1}\right)=0,\hspace{1em}n=1,3,5,\dots ,2N+1,\\ m_2{\ddot{u}}_n+c\left(2u_n-u_{n+1}-u_{n-1}\right)=0,\hspace{1em}n=0,2,4,\dots ,2N,\end{array}$$

(1)

where $(\dot{u})\equiv \frac{du}{dt}$.

The lattice is fixed at the ends

$$u_0=u_{2N+1}=0.$$

(2)

Substituting the solution in the form

$$\begin{array}{l}{u}_n=A\sin \beta {\mathrm{e}}^{i\omega t},\hspace{1em}n=1,3,5,\dots ,2N+1,\\ u_n=B\sin \beta {\mathrm{e}}^{i\omega t},\hspace{1em}n=0,2,4,\dots ,2N\end{array}$$

(3)

into Equation (1) and eliminating the amplitudes $A$ and $B$, we find the relations:

$$\begin{array}{l}\left(-{\omega }^2+{\omega }_1^2\right)A-\frac{1}{2}{\omega }_1^2\cos \frac{n\pi }{2N+1}B=0,\\ \left(-{\omega }^2+{\omega }_2^2\right)B-\frac{1}{2}{\omega }_2^2\cos \frac{n\pi }{2N+1}A=0,\end{array}$$

(4)

where $\beta =\frac{n\pi }{2N+1}$, ${\omega }_i=\sqrt{2c/m_i}$, $i=1,2$.

It is worth noting that, for the infinite lattice, the variable $\beta$ should be considered as continuously changing in the interval $0\le \beta \le \pi$.

The dispersion equation can be written as follows

$$\begin{array}{l}{\omega }^4-\left({\omega }_1^2+{\omega }_2^2\right){\omega }^2+{\omega }_1^2{\omega }_2^2{\sin }^2\hspace{0.17em}\left[\pi n/\left(2N+1\right)\right]=0,\\ n=0,1,2,\dots ,2N-1,2N,2N+1.\end{array}$$

(5)

The question of interest is as follows: in which interval of wavenumber $\beta =\frac{n\pi }{2N+1}$ changes should we consider solutions for the dispersion Equation (5)? As indicated in the articles [23,24,25], we can limit ourselves to the interval $0\le \beta \le \pi /2$ when determining the oscillation frequencies. For a complete study of the system, the interval $0\le \beta \le \pi$ should be taken into account.

Consider various limiting cases. In the following, we assume $2N\gg 1$, which is justified since we consider the possibilities of a continuous approximation to the discrete system. For $1/N\ll 1$, Equation (5) can be approximately replaced by the following:

$${\omega }^4-\left({\omega }_1^2+{\omega }_2^2\right)\hspace{0.17em}{\omega }^2+\frac{{\pi }^2{n}^2}{{\left(2N+1\right)}^2}{\omega }_1^2{\omega }_2^2=0.$$

(6)

For n = N, n = [N/2] (here symbol [] means the integer part of the number) and n = 2N + 1, we obtain from (5)

$${\omega }^2\approx {\omega }_1^2;\hspace{0.33em}{\omega }^2\approx {\omega }_2^2,$$

(7)

$${\omega }^2\approx 0.5\hspace{0.17em}\left({\omega }_1^2+{\omega }_2^2\pm \sqrt{{\omega }_1^4+{\omega }_2^4}\right),$$

(8)

$${\omega }^2=0 \mathrm{and} {\omega }^2={\omega }_1^2+{\omega }_2^2,$$

(9)

respectively.

3. Continuous Approximations

We suppose that the masses are uniformly distributed along the interval of length L of the X-axis (Figure 1). Then, the distance between the masses is $H=L/\left(2N+1\right)$. Assuming that the number of masses is large, namely $2N\gg 1$, we can introduce a natural small parameter $h=H/L\ll 1$. Let us introduce dimensionless spatial coordinate $x=X/L.$ We also introduce functions of spatial and time arguments $U_j\left(x,\hspace{0.17em}t\right),\hspace{0.33em}j=1,2$ as follows

$$\begin{gathered} U_1\left(hk,t\right)=u_k\left(t\right),k=1,3,5,\dots , \\ {U}_2\left(hk,t\right)=u_k\left(t\right),k=0,2,4,\dots \end{gathered}$$

(10)

“The basic idea is to establish a one-to-one correspondence between the functions of discrete arguments and a class of analytical functions, as well as between the operations on them. Let the function $u\left(nh\right)$ be given by its values, which are generally complex, at the points $nh$ (nodes) of the X-axis. Let us consider the problem of interpolation of $u\left(nh\right)$ by a smooth function $u\left(x\right)$. Obviously, $u\left(x\right)$ is defined within any $\Psi \left(x\right)$, which vanishes in all nodes, i.e., $u\left(x\right)+\Psi \left(x\right)$ will also be an interpolating function. It is obvious to try to select the smoothest interpolating functions from the entirety of interpolating functions filtering out rapidly oscillating components” [16]. (See also the explanation of this question in [30], p. 59). In [16] the concept of the quasicontinuum is introduced and the Whittaker–Shannon–Kotel’nikov function (sinc-function) of the form

$$sinc(x)=\frac{\sin \hspace{0.17em}\left(\pi x/h\right)}{\pi x}$$

(11)

is utilized.

The functions $U_j\left(x,t\right),\hspace{0.17em}j=1,2$ are assumed to be infinitely differentiable with respect to the spatial coordinate. Hence, they can be expanded into the Taylor series

$$U_j\left(x_0\pm h,t\right)=U_j\left(x_0,t\right)+{\displaystyle \sum _{i=1}^{\infty }{\left(-1\right)}^i{h}^i}{\left.\frac{d^{\hspace{0.17em}i\hspace{0.17em}}{U}_j\left(x,t\right)}{d{x}^i}\right|}_{\hspace{0.17em}x=x_0}.$$

(12)

Rewrite system (1) in the form

$$\begin{array}{l}{m}_1{\ddot{U}}_1\left(kh,t\right)+2c\hspace{0.17em}\left[U_1\left(kh,t\right)-U_2\left(kh,t\right)\right]\\ \hspace{1em}-c\hspace{0.17em}\left[U_2\left(\left(k+1\right)h,t\right)-2U_2\left(kh,t\right)+U_2\left(\left(k-1\right)h,t\right)\right]=0,\hspace{1em}k=1,3,5,\dots \\ m_2{\ddot{U}}_2\left(kh,t\right)+2c\hspace{0.17em}\left[U_2\left(kh,t\right)-U_1\left(kh,t\right)\right]\\ \hspace{1em}-c\hspace{0.17em}\left[U_1\left(\left(k+1\right)h,t\right)-2U_1\left(kh,t\right)+U_1\left(\left(k-1\right)h,t\right)\right]=0,\hspace{1em}k=0,2,4,\dots \end{array}$$

(13)

Using the pseudo-differential operators and sinc-function (11), we can represent system (13) as follows

$$\begin{array}{l}\left\{m_1{\ddot{U}}_1\left(x,t\right)+2c\hspace{0.17em}\left[U_1\left(x,t\right)-U_2\left(x,t\right)\right]\right\}\Phi \left(x\right)+4c{\sin }^2\left(-\frac{ih}{2}\frac{\partial }{\partial x}\right)U_2\left(x,t\right)=0,\\ \left\{m_2{\ddot{U}}_2\left(x,t\right)+2c\hspace{0.17em}\left[U_2\left(x,t\right)-U_1\left(x,t\right)\right]\right\}\Phi \left(x\right)+4c{\sin }^2\left(-\frac{ih}{2}\frac{\partial }{\partial x}\right)U_1\left(x,t\right)=0.\end{array}$$

(14)

Here, $4{\sin }^2\left(-\frac{ih}{2}\frac{\partial }{\partial x}\right)$ is the pseudo-differential operator [25], $\Phi \left(x\right)={\displaystyle \sum _{n=0}^{2N+1}\sin c\left(x-nh\right).}$ For continualization, we use the expansion

$$\begin{array}{l}{\sin }^2\left(-\frac{ih}{2}\frac{\partial }{\partial x}\right)=-\frac{1}{2}{\displaystyle \sum _{k=1}\frac{h^{2k}}{\left(2k\right)!}\frac{{\partial }^{2k}}{\partial x^{2k}}}\\ \hspace{1em}=-\frac{h^2}{4}\frac{{\partial }^2}{\partial x^2}\left(1+\frac{h^2}{12}\frac{{\partial }^2}{\partial x^2}+\frac{h^4}{360}\frac{{\partial }^4}{\partial x^4}+\frac{h^6}{20160}\frac{{\partial }^6}{\partial x^6}+\dots \right)\end{array}$$

(15)

and the equality

$$\Phi \left(x\right)\to 1 \mathrm{at} h\to 0.$$

(16)

In the zero approximation, where only terms of $h^0$ order are considered, we obtain a system of two coupled oscillators

$$\begin{array}{l}{\ddot{U}}_1+{\omega }_1^2\left(U_1-U_2\right)=0,\\ {\ddot{U}}_2+{\omega }_2^2\left(U_2-U_1\right)=0.\end{array}$$

(17)

Keeping the terms up to $h^4$ in expansion (12), we obtain

$$\begin{array}{l}{\ddot{U}}_1+{\omega }_1^2\left(U_1-U_2\right)-\frac{1}{2}{\omega }_1^2{h}^2\left[U_{2xx}+(1/12)h^2{U}_{2xxxx}\right]=0,\\ {\ddot{U}}_2+{\omega }_2^{2}2c\left(U_2-U_1\right)-\frac{1}{2}{\omega }_2^2{h}^2\left[U_{1xx}+(1/12)h^2{U}_{1xxxx}\right]=0.\end{array}$$

(18)

The comparison of the calculation results using Equation (18) with the solution of the dispersion Equation (5) demonstrates the low accuracy of this continuum approximation. In addition, as shown in [27,31], for Equation (18), the effect of the short wavelength instability take place. This artifact can lead to instability or the divergence of numerical algorithms when used to solve Equation (18).

Increasing the accuracy of continuous approximation and the short wavelength instability suppression can be achieved using Padé approximants [27,31]

$$1+\frac{h^2}{12}\frac{{\partial }^2}{\partial x^2}\sim 1/\left(1-\frac{h^2}{12}\frac{{\partial }^2}{\partial x^2}\right).$$

(19)

Hence, we obtain a “model with modified inertia” (the term was proposed in [32])

$$\begin{array}{l}\left(1-\frac{h^2}{12}\frac{{\partial }^2}{\partial x^2}\right)\left[{\ddot{U}}_1+{\omega }_1^2\left(U_1-U_2\right)\right]-\frac{1}{2}{\omega }_1^2{h}^2{U}_{2xx}=0,\\ \left(1-\frac{h^2}{12}\frac{{\partial }^2}{\partial x^2}\right)\left[{\ddot{U}}_2+{\omega }_2^2\left(U_2-U_1\right)\right]-\frac{1}{2}{\omega }_2^2{h}^2{U}_{1xx}=0.\end{array}$$

(20)

The boundary conditions for system (20) are formulated as follows

$$U_j=0 \mathrm{at} x=0,1.$$

(21)

Boundary value problem (20), (21) describes well-posed continuous approximations of original discrete systems.

The natural vibrations of the continuous system can be described as

$$U_j=C_j\sin \hspace{0.17em}\left(n\pi x\right)\exp \hspace{0.17em}\left(i\omega \hspace{0.17em}t\right).$$

(22)

System (20) approximates the oscillation frequencies of a discrete system quite well [8]. It is sufficient to consider the system on the interval $0\le \beta \le \pi /2$. However, as mentioned in [23,24,25], the original system should be considered on the interval $0\le \beta \le \pi$. In this case, system (20) does not provide a good approximation to the dispersion curve of the discrete system.

Let us construct a higher-order approximation system using the conditions for matching the solutions to the dispersion equations of the discrete and continuous systems at n = N, n = [N/2] and n = 2N + 1 (see (7)–(9)). This rational approximation is sometimes called the Padé approximation of the second kind [28,29]. As a result, we obtain a continuous system with modified inertia

$$\begin{array}{l}\left(1+d_1{h}^2\frac{{\partial }^2}{\partial x^2}+d_2{h}^4\frac{{\partial }^4}{\partial x^4}+d_3{h}^6\frac{{\partial }^6}{\partial x^6}\right)\hspace{0.17em}\left[\frac{{\partial }^2}{\partial t^2}{U}_1+{\omega }_1^2\left(U_1-U_2\right)\right]-\frac{1}{2}{\omega }_1^2{h}^2\frac{{\partial }^2}{\partial x^2}{U}_2=0,\\ \left(1+d_1{h}^2\frac{{\partial }^2}{\partial x^2}+d_2{h}^4\frac{{\partial }^4}{\partial x^4}+d_3{h}^6\frac{{\partial }^6}{\partial x^6}\right)\hspace{0.17em}\left[\frac{{\partial }^2}{\partial t^2}{U}_2+{\omega }_2^2\left(U_2-U_1\right)\right]-\frac{1}{2}{\omega }_2^2{h}^2\frac{{\partial }^2}{\partial x^2}{U}_1=0,\end{array}$$

(23)

where

$$\begin{array}{l}{d}_1=\frac{73}{140}+\frac{61}{9{\pi }^2}-\frac{32\sqrt{2}}{35};\hspace{1em}{d}_2=\frac{1}{63\hspace{0.17em}{\pi }^4}\hspace{0.17em}\left[812+\left(327-288\sqrt{2}\right){\pi }^2\right];\\ d_3=\frac{16}{315\hspace{0.17em}{\pi }^6}\left[140+\left(87-72\sqrt{2}\right){\pi }^2\right].\end{array}$$

(24)

Unfortunately, system (22) requires additional boundary conditions. This is a typical difficulty that arises when higher-order models, derived originally for infinite media, are applied to bounded domains [5,14,15]. To overcome this fundamental difficulty, various approaches are used.

In [33], this problem was studied for 1D linear lattice, when the exact solution can be easily obtained. It is proposed to introduce boundary conditions so as to achieve the most suitable frequencies and mode approximations. The limitations of this method are clear.

An often used technique is the formation of some edge zone. Thus, in [34,35], for second-grade elasticity, it is proposed to introduce the so-called ortho-fiber—a thin, notional fiber of a material, starting at the surface and extending inward along the direction of normal. The length of the ortho-fiber is small compared to the macroscopic characteristic length of the problem and tends to zero.

For continuous systems (beams, plates), the correct formulation of boundary conditions requires the consideration of the methods of fastening the plates used in practice (for example, embedding in a groove). For the problem of bending a plate, when deriving and justifying the conditions of simply support, it is assumed that there is a narrow zone of pinching of the plate edge, with the length of it tending towards zero as some power the thickness of the plate [36].

Similar problems occur in grid methods (see overview in [5]), where they can be solved by introducing some fictitious points outside the considered domain.

All these techniques are based on the simple consideration that, in reality, the conditions of ends fixation have a fairly complex structure, which should be taken into account when studying continuous or discrete systems.

Using this ideology, we introduce fictitious masses at points $k=-1,-2,-3,\dots$; $k=2N+2,2N+3,2N+4,\dots$ (Figure 2).

Particularly in the case of a periodic spatial extension (“simple support”), one obtains (see Figure 2):

$$u_i=-u_{-i},\hspace{0.33em}i=1,2,3,\dots ;u_{2N-j}=-u_{2N+2+j},\hspace{0.33em}j=0,1,2,3,\dots$$

(25)

In the continualization problems, this technique was applied in [14,15]. Using it for the considered case, we arrive at the following boundary conditions

$$U_1=U_{1xx}=U_{1xxxx}=0,\hspace{0.33em}{U}_2=U_{2xx}=U_{2xxxx}=0 \mathrm{at} x=0,1.$$

(26)

From system (23), we obtain the dispersion equation for the continuous system

$$\begin{array}{l}{\omega }^4-\left({\omega }_1^2+{\omega }_2^2\right){\omega }^2+{\omega }_1^2{\omega }_2^2\frac{{\left(\frac{n\pi }{2N+1}\right)}^2}{1-d_1{\left(\frac{n\pi }{2N+1}\right)}^2+d_2{\left(\frac{n\pi }{2N+1}\right)}^4-d_3{\left(\frac{n\pi }{2N+1}\right)}^6}\hspace{0.17em}\times \\ \hspace{1em}\times \hspace{0.17em}\left[1-\frac{{\left(\frac{n\pi }{2N+1}\right)}^2}{4\hspace{0.17em}\left[1-d_1{\left(\frac{n\pi }{2N+1}\right)}^2+d_2{\left(\frac{n\pi }{2N+1}\right)}^4-d_3{\left(\frac{n\pi }{2N+1}\right)}^6\right]}\right]=0.\end{array}$$

(27)

The analysis of Equation (27) shows that the system (23), (24) is unconditionally stable.

Note that the boundary value problem of the continuous approximation (23), (26) remains the same for the case of an infinite lattice.

4. Numerical Results

The results of comparing the dispersion curves for the discrete and continuous systems at ${\omega }_2^2/{\omega }_1^2=1.5$ are shown in Figure 3.

The solid curves show branches with normal dispersion, and the dotted curves show those with anomalous dispersion. The calculations were carried out at the points $n=0,1,2,\dots ,2N,2N+1$. The value $2\beta /\pi$ is plotted along the abscissa.

As can be seen from Figure 3, the continuous system (23), (24) provides a good approximation to the dispersion curve of the discrete system over the entire range of wavenumber changes. Recall that we assume $2N\gg 1$, and the distance between the masses is small. Hence, the dispersion curves in Figure 3 are shown as continuous lines. Note that the curves corresponding to the discrete and continuous models practically merge on the scale of the graph.

The variables $U_1, U_2$ have no obvious physical meaning (variables describing low-frequency oscillations of the centers of mass of the cells and high-frequency oscillations associated with internal degrees of freedom have physical meaning [9]).

The ratio of the amplitudes of mass vibrations $C_1$ and $C_2$ is given by the following relations:

$$C_1/C_2=\frac{\left(-A_2+A_1\right){\omega }_1^2}{A_1\left(-{\omega }^2+{\omega }_1^2\right)} \mathrm{at} {\omega }^2\le {\omega }_1^2$$

(28)

and

$$C_2/C_1=\frac{\left(-A_2+A_1\right){\omega }_2^2}{A_1\left(-{\omega }^2+{\omega }_2^2\right)} \mathrm{at} {\omega }_2^2\le {\omega }^2,$$

(29)

where $A_1=1-d_1{\left(\frac{n\pi }{2N+1}\right)}^2+d_2{\left(\frac{n\pi }{2N+1}\right)}^4-d_3{\left(\frac{n\pi }{2N+1}\right)}^6;\hspace{0.33em}{A}_2=\frac{1}{2}{\left(\frac{n\pi }{2N+1}\right)}^2$.

The calculation results at ${\omega }_2^2/{\omega }_1^2=1.5$ are presented in Figure 4. The obtained curves qualitatively coincide with the corresponding curves for the discrete system [24] and make it possible to estimate the ratios of the amplitudes of the neighboring masses.

In the limiting case $m_1=m_2$, the dispersion curve for the Lagrange lattice is obtained (solid curve in Figure 5). The dotted curve corresponds to the displacements with zero amplitude [24,25].

Dispersion curves for the discrete (1) (curves 1 and 2) and continuous (23), (24) (curves 3 and 4) systems are presented at ${\omega }_2^2/{\omega }_1^2=1$. Note that the curves corresponding to the discrete and continuous models practically merge on the scale of the graph.

We also compare the eigenfrequency dependence with respect to the mode number for the finite diatomic lattice, both with the discrete and the continuous approximation proposed in the paper. Such a comparison is uninformative because the frequency values are of different orders. Comparing the values of frequencies ${\omega }_d$ and ${\omega }_c$ in some averaged sense looks more reasonable. For that, we use the following expression for the average approximation error [37]:

$$A=\frac{100\%}{{\overline{\omega }}_d}\sqrt{\frac{{\displaystyle \sum _{i=0}^{2N+1}\hspace{0.17em}{\left({\omega }_{id}-{\omega }_{ic}\right)}^2}}{2N+2}},$$

(30)

where ${\overline{\omega }}_d=\frac{{\displaystyle \sum _{i=0}^{2N+1}\hspace{0.17em}{\omega }_{id}}}{2N+2}$.

The calculations were carried out at ${\omega }_2^2/{\omega }_1^2=1.5$. The maximum error was obtained for the acoustic branch (Figure 3). The calculation results for various values of N are given in

Table 1. The average error of the continuous approximations for frequencies.

N	5	10	15	20	25	30	35	40	45	50	∞
A, %	0.4465	0.5837	0.6272	0.6471	0.6579	0.6645	0.6688	0.6718	0.6740	0.6755	0.6824

.

When calculating the limit, summation is replaced by integration for a sufficiently large N as follows

$$\begin{array}{l}{\overline{\omega }}_d\approx \frac{1}{2}{\displaystyle \underset{0}{\overset{2}{\int }}{\omega }_d}\left(n_1\right)\hspace{0.17em}d{n}_1,\\ \frac{{\displaystyle \sum _{i=0}^{2N+1}\hspace{0.17em}{\left({\omega }_{id}-{\omega }_{ic}\right)}^2}}{2N+2}\approx \frac{1}{2}{\displaystyle \underset{0}{\overset{2}{\int }}{\left({\omega }_{id}\left(n_1\right)-{\omega }_{ic}\left(n_1\right)\right)}^2}\hspace{0.17em}d{n}_1,\hspace{1em}{n}_1=\frac{n\pi }{N+0.5}.\end{array}$$

(31)

Hence, we obtain the following estimate:

$$A\approx \frac{\sqrt{2\hspace{0.17em}{\displaystyle \underset{0}{\overset{2}{\int }}{\left({\omega }_{id}\left(n_1\right)-{\omega }_{ic}\left(n_1\right)\right)}^2}\hspace{0.17em}d{n}_1}}{{\displaystyle \underset{0}{\overset{2}{\int }}{\omega }_d}\left(n_1\right)\hspace{0.17em}d{n}_1}\cdot 100\%$$

(32)

It is obvious that the approximation of frequencies using the continuous approximation is quite satisfactory.

5. Conclusions

The continualization of a diatomic lattice was considered in the natural interval of wavenumber changes. It provided a consistent limiting transition to the Lagrange lattice.

Note that the original system of 2N + 2 difference-differential Equation (1), each of which contains the second-order derivative with respect to time, can be reduced to a system of N + 1 difference-differential equations with the fourth-order derivative with respect to time. After conducting its continualization, we obtain one continuous equation with the fourth-order derivative with respect to time. This result agrees with the fourth-order equation with respect to the time obtained from the continuous approximation (23).

The continualization consisted in the replacement of the original pseudo-differential equations by a system of PDEs that provided a good approximation for the dispersion relations. The use of the Padé approximants guaranteed the matching of the dispersion relation values of the discrete and continuous models at several characteristic points. It was not possible to limit ourselves to the second-order system of PDEs, which led to the need for the correct formulation of boundary conditions. An adequate result was obtained using a sixth-order system of PDEs with respect to spatial variables with modified inertia, which is unconditionally stable.

The application of the continuous approximations allowed correctly describing not only the frequency spectrum in the first Brillouin zone, but also the amplitude ratios of neighboring masses.

The technique proposed in this paper can be used to construct continuous models for flexural waves in the diatomic lattices [12]. The obtained continuous model can be also applied to the nonlinear problems for the diatomic lattice described in [9,11,26].