Improved Continuous Models for Discrete Media

The paper focuses on continuous models derived from a discrete microstructure. Various continualization procedures that take into account the nonlocal interaction between variables of the discrete media are analysed.


Introduction
In the recent years new classes of ultra-dispersive and nanocrystalline materials 1-3 , which require a different modern approach than that of a classical continuous media, have been proposed. For example, the nanocrystalline material is represented by a regular or quasiregular lattice with small size bodies domains, granules, fullerenes, nanotubes, or clusters nanoparticles possessing internal degrees of freedom occupying the lattice sites 4 . The situation mentioned is also a characteristic for various problems of nanomechanics 5, 6 , because the transition of a material to the nanostructural state is accompanied by dimensional effects in its mechanical properties. Models established on a classical continuous media cannot govern high frequency vibrations, material behavior in the vicinity of cracks and on the fronts of destruction waves, and during phase transitions 7 . It is not possible to reach and overcome bifurcation points, that is, thresholds of lattice stability under catastrophic deformations 8 . Wave dispersion in granular materials 9-11 represents an important example of microstructure effects too. Microstructural effects are essential for correct description of softening phenomena 12 in damage mechanics [13][14][15][16] and in the theory of plasticity 17, 18 . As it is mentioned in 19 , "From the behaviour of calcium waves in living cells to the discontinuous propagation of action potentials in the heart or chains of neurons and from chains of chemical reactors or arrays of Josephson junctions to optical waveguides, dislocations and the DNA double strand, the relevant models of physical reality are inherently discrete."

Mathematical Problems in Engineering
The effects mentioned may be analyzed within the frame of discrete models, using molecular dynamics 20 , quasicontinuum analysis 21, 22 , or other numerical approaches. However, a sought result is difficult to obtain using high tech computers in an economical way. For example, modern practical problems are still intractable for molecular dynamicsbased analysis, even if the highest computing facility is at disposal.
Situation can be described by Dirac's words 23 : "The physical laws necessary for the mathematical theory of a large part of physics and whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be solvable. Therefore, it becomes desirable that approximate practical methods of applying quantum mechanics should be developed, which can lead to an explanation of the main features of complex atomic systems without too much computations." Hence, refinement of the existing theory of continuous media for the purpose of more realistic predictions seems to be the only viable alternative. In connection with this, one of the most challenging problems in multiscale analysis is that of finding continuous models for discrete, atomistic models. In statistical physics, these questions were already addressed 100 years ago, but many problems remain open even today. Most prominent is the question how to obtain irreversible thermodynamics as a macroscopic limit from microscopic models that are reversible. In this paper we consider another part of this field that is far from thermodynamic fluctuation. We are interested in reversible, macroscopic limits of atomic models. Debye approach is the simplest model of this type 24 , but it does not take into account spatial dispersion.
Therefore, continuous modeling of micro-and nanoeffects plays a crucial role in mechanics. It seems that the simplest approach to realize this idea relies on a modification of the classical modeling keeping both hypothesis of continuity and main characteristic properties of a discrete structure. Here four strategies exist.
Phenomenological approach: additional terms are added to the energy functional or to the constitutive relation. The structure and character of these terms are postulated 25-31 . Phenomenological approach is very useful and accurate enough in the applied sciences when it is necessary to solve problems of practical significance 32-34 . But progress of natural sciences makes us look for ways of substantiated derivation of constitutive relations from "first principles." In this connection one can recall the 6th Hilbert's problem mathematical treatment of the axioms of physics 35 : "Boltzmann's work on the principles of mechanics suggests the problem of developing mathematically the limiting processes, there merely indicated, which lead from the atomistic view to the laws of motion of continua. Conversely, one might try to derive the laws of the rigid bodies' motion by a limiting process from a system of axioms depending upon the idea of continuously varying conditions of a material filling all space continuously, these conditions being defined by parameters. The question on the equivalence of different systems of axioms is always of great theoretical interest." Statistical approach: starting from an inhomogeneous classical continuum, average values of the state variables are computed to produce enhanced field equations 36 . Unfortunately, great mathematical problems can not give possibility of wide usage of this approach.
Homogenization approach is based on Γand G-limit technique 37-43 and gave many important pure mathematical results, but not new insights for physics.
Continualization procedures are based on various approximations of local discrete operator by the nonlocal one. For this aim Taylor series 44, 45 , one-and two-point Padé approximants 46-62 , composite equations 63, 64 , and other approaches are used. Independently of the strategy employed, the constitutive relation of the homogenized material becomes nonlocal. Mathematically this implies integral relations between stresses and strains, in which the stress in a point is not only related to the strain in the same point but also instantly to strains in the neighbouring points. We will analyse these approaches further. The paper is organized as follows. Section 2 is devoted to some general remarks concerning a chain of elastically linear coupled masses. In Sections 3 and 4 the classical continuous approximation and so-called "splash effect" for 1D case are described. Section 5 is devoted to anticontinuum limit. Improved continuous approximations are studied is Sections 6 and 7 for natural and forced oscillations. Waves in 1D discrete and continuous media are compared in Section 8. 2D problems are analysed in Section 9. Nonlinear phenomenons are studied in Section 10. In Section 11 we analyse some possible generalization of described approaches. Section 12 is devoted to Navier-Stokes equations. Section 13 presents brief concluding remarks. In appendixes to the paper we describe one-and two-point Padé approximations, continuum limit of Toda lattice, and correspondence between functions of discrete arguments and approximating analytical functions.

A Chain of Elastically Coupled Masses
In this section we follow paper 65 . We study a chain of n 1 material points with the same masses m, located in equilibrium states in the points of the axis x with coordinates jh j 0, 1, . . . , n, n 1 and suspended by elastic couplings of stiffness c Figure 1 .
Owing to the Hooke's law the elastic force acting on the jth mass is as follows: where y j t is the displacement of the jth material point from its static equilibrium position. Applying the 2nd Newton's law one gets the following system of ODEs governing chain dynamics: In general, the initial conditions have the following form: As it has been shown in 65  where constants C j are defined via solution of the following eigenvalue problem: Function T t satisfies the following equation: mT tt cλT 0.

2.9
A solution to 2.8 has the form T A exp iωt . Hence, 2.8 and 2.9 yield the following Lagrange formula for frequencies ω k of discrete system: , k 1, 2, . . . , n.

2.10
Since all values λ k are distinct, then all eigenvalues are different. Therefore, each of the eigenvalues is associated with one eigenvector C k C Each of the eigenfrequencies 2.10 is associated with a normal vibration y k j t C k j A k cos ω k t B k sin ω k t , k 1, 2, . . . , n.

2.13
A general solution of the BVP 2.3 -2.5 is obtained as a result of superposition of normal vibrations 2.14 Let us study now the problem of chain masses movement under action of a unit constant force on the point number zero. Motion of such system is governed by 2.3 with the following boundary and initial conditions: σ j t σ jt t 0 for t 0.

2.15
In what follows the initial BVP with nonhomogeneous boundary conditions 2.3 , 2.15 will be reduced to that of homogeneous boundary condition for 2.3 with nonhomogeneous initial conditions, and then the method of superposition of normal vibrations can be applied. The formulas of normal forms obtained so far can be applied in a similar way with exchange of y j t for σ j t . In effect, the following exact solution to the BVP 2.3 , 2.15 is obtained 65 :

Classical Continuous Approximations
For large values of n usually continuous approximation of discrete problem is applied. In our case, described by 2.3 , 2.15 , it takes the form of where l n 1 h.

6
Mathematical Problems in Engineering BVP 3.1 -3.3 can be used, for example, for modeling of stresses in van couplings of the rolling stocks 66 .
Having in hand a solution to continuous BVP 3.1 -3.3 , one obtains a solution of a discrete problem due to the following formulas: σ j t σ jh, t , j 0, 1, . . . , n, n 1.

3.4
Formally, the approximation described so far can be obtained in the following way. Let us denote the difference operator occurring in 2.3 as D, that is, Applying the translation operator exp h∂/∂x , one gets 67 Let us explain 3.6 in more details. The McLaurin formula for infinitely many times differentiable function F x has the following form: Observe that exp h∂/∂x belongs to the so-called pseudodifferential operators. Using relations 3.5 -3.7 , we cast 2.3 into pseudo-differential equation of the following form: On the other hand, splitting the pseudo-differential operator into the McLaurin series is as follows: 3.9 Keeping in right hand of 3.9 only the first term, one obtains a continuous approximation 3.1 . Note that an application of the McLaurin series implies that displacements of the neighborhood masses differ slightly from each other. From a physical point of view, it means that we study vibrations of the chain with a few masses located on the spatial period Figure 2 ; that is, the long wave approximation takes place. Note that the vertical axis in Figure 2 represents the displacement in the x direction, since the problem is 1D.
Mathematical Problems in Engineering 7 x Figure 2: Solution form σ σ x, t in a fixed time instant t = const points-discrete system, curvecontinuous system . Relations 3.10 relatively good approximate low frequencies of discrete system 2.10 , whereas the nth frequency α k of a continuous system differs from the corresponding nth frequency ω k of a discrete system of order of 50%. Accuracy of continuous approximations can be improved, what will be discussed further.

Splashes
One can obtain an exact solution to the BVP 3.1 -3.3 , using the d'Alembert method matched with operational calculus 65 : where H · · · is the Heaviside function. From 4.1 one obtains the following estimation: It was believed that estimation 4.2 with a help of relation 3.4 can be applied also to a discrete system 68 . However, analytical as well as numerical investigations 66, 69-73 indicated a need to distinguish between global and local characteristics of a discrete system. In other words during investigation of lower frequency part a transition into continuous model is allowed. However, in the case of forced oscillations solutions to a discrete system may not be continuously transited into solution of a wave equation for h → 0 67 . Numerical investigations show that for given masses in a discrete chain quantity the P j |σ j t | may exceed the values of 1 in certain time instants reported in 69, Table 5.1 .

Mathematical Problems in Engineering
Observe that splash amplitude does not depend on the parameter m/c. In addition, the amplitude of the chain vibrations increases with an increase of n; whereas its total energy does not depend on n. However, this is not a paradox. Namely, amplitude of vibrations has an order of sum of quantities σ j t ; whereas its potential energy order is represented by a sum of squares of the quantities mentioned 65 .
On the other hand, a vibration amplitude of a mass with a fixed number is bounded for n → ∞, but amplitude of vibrations of a mass with a certain number increasing with increase of n tends to infinity for n → ∞ following ln n 65 . "In the language of mechanics what we just said means that when analyzing the so-called "local properties" of a one-dimensional continuous medium, one cannot treat the medium as the limiting case of a linear chain of point masses, obtained when the number of points increases without limit" 66 .
It should be emphasized that a rigorous proof of the mentioned properties has been obtained for a case, when n 1 is a simple digit or it is a power of 2. However, this assumption is not necessary, as it is mentioned in 65 .
Earlier the same effect of continualization was predicted by Ulam, who wrote 74, pages 89, 90 : "The simplest problems involving an actual infinity of particles in distributions of matter appear already in classical mechanics." A discussion on these will permit us to introduce more general schemes which may possibly be useful in future physical theories.
Strictly speaking, one has to consider a true infinity in the distribution of matter in all problems of the physics of continua. In the classical treatment, as usually given in textbooks of hydrodynamics and field theory, this is, however, not really essential, and in most theories serves merely as a convenient limiting model of finite systems enabling one to use the algorithms of the calculus. The usual introduction of the continuum leaves much to be discussed and examined critically. The derivation of the equations of motion for fluids, for example, runs somewhat as follows. One images a very large number N of particles, say with equal masses constituting a net approximating the continuum, which is to be studied. The forces between these particles are assumed to be given, and one writes Lagrange equations for the motion of N particles. The finite system of ordinary differential equations becomes in the limit N ∞ one or several partial differential equations. The Newtonian laws of conservation of energy and momentum are seemingly correctly formulated for the limiting case of the continuum. There appears at once, however, at least possible objection to the unrestricted validity of this formulation. For the very fact that the limiting equations imply tacitly the continuity and differentiability of the functions describing the motion of the continuum seems to impose various constraints on the possible motions of the approximating finite systems. Indeed, at any stage of the limiting process, it is quite conceivable for two neighbouring particles to be moving in opposite directions with a relative velocity which does not need to tend to zero as N becomes infinite; whereas the continuity imposed on the solution of the limiting continuum excludes such a situation. There are, therefore, constraints on the class of possible motions which are not explicitly recognized. This means that a viscosity or other type of constraints must be introduced initially, singling out "smooth" motions from the totality of all possible ones. In some cases, therefore, the usual differential equations of hydrodynamics may constitute a misleading description of the physical process.
Splash effect was observed numerically in 2D linear and 1D nonlinear case A. M. Filimonov, private communication . By the way, due to this effect many papers "justifying" usual continualization can be treated as naive. For example, in 75 the continuous limit was derived based on the hypothesis that the microscopic displacements are equal to the macroscopic ones. In 75 authors supposed that the displacements of the particles which are connected by elastic springs are small in the following sense: |y i − y j | ≤ cε. This inequality can be justified only for lower part of spectrum for natural oscillations.

Anticontinuum Limit
A classical continuous approximation allows for relatively good description of a low part of the vibration spectrum of a finite chain of masses. In what follows we study now another limiting case, so-called anticontinuum limit, that is, completely uncoupled limit for lattice Figure 3 .
In this case one has σ k −1 k Ω, and equation for Ω reads This is so-called anticontinuum limit.
In the case of vibrations close to the saw-tooth one, the short-wave approximation is applied "envelope continualization" 76-79 Figure 4 . Namely, first we use staggered transformation and then 2.3 and 2.15 are reduced to the following BVP: Then, the following relations are applied: k 0, 1, 2, . . . , n, n 1. Using 5.6 , 5.3 , and considering h 2 as a small parameter one gets we take zeroth and first-order approximations only mΩ tt 4cΩ ch 2 Ω xx 0. 5.7 Appropriate boundary and initial conditions for 5.7 follow

5.8
Observe that practically the whole frequency interval of discrete model is well approximated for two limiting cases, that is, for the case of the chain and for the case of the envelope.

Improved Continuous Approximations
In what follows we are going to construct improved continuous approximations. Modeling of such systems nonlocal theories of elasticity requires integral or gradient formulation. The integral formulation may be reduced to a gradient form by truncating the series expansion of the nonlocality kernel in the reciprocal space 80 . In what follows we apply the gradient formulation approach.
If, in the series 3.9 , we keep three first terms, the following model is obtained: However, a nontrivial problem regarding boundary conditions for 6.1 appears 81, 82 . The conditions mentioned can be defined only assuming the chain dynamics behavior for k −1, −2, −3; k N 2, N 3, N 4. In other words, a boundary point is replaced by a boundary domain 44, 45 . In particular, in the case of periodic spatial extension "simple support" one gets σ σ xx σ xxxx 0 for x 0, l. 6.2

Mathematical Problems in Engineering 11
Assuming σ k t 0 for k −1, −2, −3; k N 2, N 3, N 4, instead of the boundary conditions 6.2 we have "clamping" σ σ x σ xxx 0 for x 0, l. 6.3 Comparison of nth frequency of a continuous system 6.1 , 6.2 with that of a discrete system exhibits essential accuracy improvement applying coefficient 2.1 instead of 2 in an exact solution yields an error of ∼5% . Observe that an estimation of the accuracy of continuous approximation on the basis of comparison of discrete and continuous systems frequencies rather simple but yielding reliable results.
In a general case, keeping in series 3.9 N terms, one gets equations of the so-called intermediate continuous models 70 Boundary conditions for 6.4 have the following form: The corresponding BVPs are correct and also stable during numerical solution for odd N. In this case 6.4 is of hyperbolic type 70 . Application of intermediate continuous models allows catching the mentioned splashes effect.
One also can mention the momentous elasticity theories 14, 25-27, 83-87 . For example, taking into account dependence of energy of deformation from the higher gradients of displacements leads to the concept of momentous stresses 88 . Le Roux was the first who showed the importance of the higher gradients of displacements 28 . In momentous theories of elasticity the method of macrocells is also widely used 89-91 . The dynamics of the continuous media is described by the equations of displacement of the centre of mass of a macro cell and by the equations of moments of various orders. The spectrum of this media tends to the complete spectrum of a linear chain with increasing number of moments. Critical reexamination of this approach can be found in 3 .
The construction of intermediate continuous models is mainly based on the development of a difference operator into Taylor series. However, very often application of Padé approximants PAs is more effective for approximation 52, 92 for description of PA see Appendix A . Collins 53 proposed to construct continuous models using PA. Then this approximation was widely used by Rosenau 55-60 . More exactly, Collins and Rosenau did not use PA straightforward; this was done later 61, 62, 93 . Sometimes these continuous models are called quasicontinuum approximations.
If only two terms are left in the series 3.9 , then the PA can be cast into the following form: For justification of this procedure Fourier or Laplace transforms can be used. The corresponding quasicontinuum model reads The boundary conditions for 6.8 have the form σ 0 for x 0, l. 6.9 The error regarding estimation of n-th frequency in comparison to that of a discrete chain is of ∼16.5%. However, 6.8 is of lower order in comparison to 6.1 . Equation 6.8 is usually called Love equation 94 but as Love mentioned 95 this equation was obtained earlier by Rayleigh 96,97 . Term σ xxtt can be treated as the lateral inertia. Equation 6.8 has hyperbolic type. Kaplunov et al. 98 refer to these type theories as theories with modified inertia. Passage to 6.8 can be treating as regularization procedure. The model governed by 6.8 is unconditionally stable and propagating waves cannot transfer energy quicker than the velocity c. However, the model governed by 6.8 predicts that short waves transfer elastic energy with almost zero speed 99 .
It is worth noting that 5.7 also can be regularized using PA as follows 100-102 : Finally, having in hand both long and short waves asymptotics, one may also apply two-point PA description of two-point PA see Appendix A . In what follows we construct two-point PA using the first term of the series 3.9 as a one of limiting cases. We suppose

6.11
Mathematical Problems in Engineering

13
The improved continuous approximation is governed by the following equation: Now we require the n-th frequency of a continuous and discrete system to coincide ω n 2 c m sin nπ 2 n 1 .

6.13
For large value of n one may approximately assume that α n ≈ 2 c/m 6.14 with the boundary conditions 6.9 .

6.15
Now, using 6.14 one gets a 2 0.25 − π −2 . The highest error in estimation of the eigenfrequencies appears for k 0.5 n 1 and consists of 3%.
Observe that equations similar to the 6.14 are already known. Eringen 103-106 using a correspondence between the dispersion curves of the continuous approximation and Born-von Kármán model 82, 107 obtained a 2 0.1521. This value is very close to that proposed in 108, 109 on the basis of a certain physical interpretation, where also model 6.14 is referred to as the "dynamically consistent model." Interesting, that Mindlin and Herrmann used very similar to two-point PA idea for construction their well-known equation for longitudinal waves in rod 110

Forced Oscillations
In order to study forced vibrations we begin with a classical continuous approximation. A solution to 3.1 is sought in the following form: and the function u x, t is defined by the following BVP: Solution to the BVP 7.2 is obtained via Fourier method where α k defined by 3.10 . For the approximation the motion of the chain, one must keep in an infinite sum only n first harmonics Observe that solutions 2.16 and 7.4 differ regarding frequencies α k and ω k relations 2.8 and 3.10 , resp. . In addition, the coefficients in series 2.16 and 7.4 differ from each other, that is, projections onto normal vibrations of discrete and continuous systems, 1/ n 1 ctg πk/2 n 1 and 2/πk, respectively, are strongly different for k 1. This is because during projection into normal vibrations for the discrete system one must use summation due k from 0 to n, whereas for the continuous system-integration due x from 0 to l. The problem occurring so far can be overcome using the Euler-McLaurin formulas 114 where B i are the Bernoulli numbers, having the following values: B 0 1, B 1 −1/2, B 2 1/6, B 3 0, . . . .
In addition, the following formulas hold 115 : The corresponding integrals read Owing to 7.11 , one can construct a simple expression relatively good approximating sum 7.8 for arbitrary j values from j 1 to j n. For this purpose we change second term of the right-hand side of 7.11 in the following way:

Wave Processes in Discrete and Continuous 1D Media
Let us study wave propagation in the infinite 1D medium Figure 5 . Assuming that when time instant t 0, a force P acts on a mass with number 0 in direction of the axis x. Governing equations can be written in the following way 116 :

8.1
We reduce 8.1 to the dimensionless form

8.2
Here τ t c/m; Y j y j /h; P 1 P/ ch . Due to the linearity of the problem one can suppose P 1 1. Let us apply Fourier transform to 8.2 116 . Namely, multiplying left-and right-hand sides of 8.2 by exp iqj , j 0, ±1, ±2, . . . and adding them one obtains Observe that the term 4sin 2 q/2 can be developed in the vicinity of either q 0 classical continuous approximation 4sin 2 q 2 q 2 · · · 8.4 or q π envelope continualization 4sin 2 q 2 4 · · · . 8.5 The obtained solution is exact one and will be further used for the error estimation of improved continuous models.
Applying the classical continuous approximation instead of system 8.2 , the following wave equation with the Dirac delta-function in right hand-side is obtained: Here z x/h. The following relation between discrete and continuous systems holds: Now, applying the Fourier transform in x and the Laplace transform in τ, one obtains So, a wave propagation in a discrete media strongly differs from this phenomenon in the classical continuous media. In what follows we give solutions on a basis of the theories of elasticity described so far. The intermediate continuous model is as follows: Let us explain the occurrence of term sin πz / πz in the right hand of 8.10 instead of Dirac delta-function for detail description see Appendix C .
Applying both Fourier transform in z and Laplace transform in τ, the following equation governing wave velocity propagation is obtained: The quasicontinuum model is as follows: 1 − 1 12 Model 6.14 for our case has the following form: In this case the associated wave velocity propagation reads

∂Y z, τ ∂τ
Since analytical comparison of the obtained wave velocity propagation is difficult, we apply numerical procedures.
As it has been already mentioned, exact solution of the joined discrete problem governed by relation A 1 j, τ will serve us as the standard solution used for estimation of other solutions. Results of computations are shown on Figures 6 a -6 c and Figures 7 a -7 c . One can conclude that above described improved continuous models give possibility qualitatively taking into account effect of media discreteness: i occurrence of oscillations, being gradually damped at x const; ii infinite velocity of perturbations propagation owing to assumption of instantaneous masses interaction during their draw near ; iii occurrence of the quasi-front domain, where the stresses increase relatively fast, but without jump, as on the wave front set; velocities and deformations exponentially decay while the quasi-front τ − |z| increases finally becoming negligibly small. Equation 8.10 has sixth order in spatial variable, 8.12 and 8.14 have second order in spatial variable, so, they are simpler for practice applications. Equations 8.12 , and 8.14 give for case of wave motion qualitatively the similar results.

2D Lattice
In order to analyse the 2D case we use 9-cell square lattice Figure 8

9.1
Here u i·j , v i·j is the displacement vector for a particle situated at point x i , y j , x i ih, y j jh.
The standard continualization procedure for 9.1 involves introducing a continuous displacement field u x, y , v x, y such that u x m , y n u m,n ; v x m , y n v m,n , and expanding u m±1,n±1 , v m±1,n±1 into Taylor series around u m,n , v m,n . Second-order continuous theory in respect to the small parameter h is 87

9.2
Mathematical Problems in Engineering 21 Naturally, one can construct equations of higher order, but, as it is shown in 87 , generally it is not possible to create asymptotic theories that do not possess extraneous solutions in the nonscalar context.
Using staggered transformations one obtains in anticontinuum limit.

9.4
As in 1D case see 5.1 equations 9.4 do not contain parameter h 2 . For this type of 2D lattice equations in anticontinuum limit are decoupled.
The existence of the continuous approximations 9.2 and 9.4 gives a possibility to construct the composite equations, which is uniformly suitable in the whole interval of the frequencies and the oscillation forms of the 2D lattice of masses. Let us emphasize, that the composite equations, due to Van Dyke 118 , can be obtained as a result of synthesis of the limiting cases. The principal idea of the method of the composite equations can be formulated in the following way 118, page 195 .
i Identify the terms in the differential equations, the neglection of which in the straightforward approximation is responsible for the nonuniformity.
ii Approximate those terms insofar as possible while retaining their essential character in the region of nonuniformity.
In our case the composite equations will be constructed in order to overlap approximately with 9.2 for long wave solution and with 9.4 for short wave solution. As a result of the described procedure one gets

Mathematical Problems in Engineering
For 1D case from 9.5 one obtains 6.14 for u and the same equation for v. For small variability in spatial and time variables equations 9.5 can be approximated by 9.2 , for large variability in spatial and time variables equations 9.5 can be approximated by 9.4 .

Nonlinear Case
Fundamental methods of analysis of linear lattices are those based on either discrete or continuous Fourier transformations 44, 45, 119 and the integral Fourier operators method proposed by Maslov 67 . In a nonlinear case the similar methods, in general, are not applicable. However, there exists class of piece-wise acting forces-springs extensions relations, opening the doors for analytical tools application 120 .
Besides, one may achieve even a solution to a nonlinear problem using the following observation 61, 62, 119, 121 . Let nonlinear springs in 1D problem have potential energy U y , where y is a difference between the neighborhood masses regarding their equilibrium positions. In this case the stretching force is −∂U y /∂y. The chain dynamics is governed by the following infinite system of ODEs: where r n y n 1 − y n . Assume that solution to the ODEs is a traveling wave, r n R z , z n−vt. Then 10.1 is cast into the following form: where · · · d · · · /dr. Introducing notation Φ z U R z and applying the Fourier transformation to 10.2 one gets In the conjugated space one obtains where f q 4 sin 2 q/2 /q 2 . Term f q can be approximated by PA or two-point PA in the following way: We show also a simultaneous application of PA matched with a perturbation procedure using an example of the Toda lattice 122 m d 2 y n dt 2 a exp b y n−1 − y n − exp b y n − y n 1 . 10.8 In continuum limit one obtains for details see Appendix B where: t 1 ab/m t; y 1 b/2 y. One can obtain from 10.9 the Korteweg-de Vries KdV equation 122 . Here, we illustrate how one may get regularized long-wave equation 123-125 . Using PA one obtains 1 1 24

10.10
Then, 10.9 yields up to the highest terms 1 − 1 24 On the other hand, there are nonlinear lattices with a special type nonlinearities allowing to achieve exact solutions as soliton and soliton-like solutions Toda, Ablowitz and Ladik, Langmuir, Calogero, and so forth, 119, 126-130 in the case of infinite lattices or in the case of occurrence of physically unrealistic boundary conditions. In many cases nonintegrable systems like a discrete variant of sine-Gordon equation possess soliton-like solution.
Lattice discreteness allows the existence of new types of localized structures that would not exist in the continuum limit. Nonlinearity transfers energy to higher wavenumbers but it can be suppressed and balanced by the bound of the spectrum of discrete systems. Such balance can generate highly localized structures, that is, the discrete breathers 100-102, 131 . Improved continuous models give possibility to construct such type of localized solutions 100-102 .

Mathematical Problems in Engineering
On the other hand, it is more appropriate changing a system of Fermi-Pasta-Ulam not by KdV equation, but by Toda lattice. In the latter case one gets an asymptotic regarding nonlinear parameter instead of applying direct substitution of a discrete system by a continuous one. Occurrence of exact solutions of a discrete system allows applying the following approach; fast changeable solution part is constructed using a discrete model, whereas a continuous approximation is applied for slow components. The latter observation is very well exhibited by Maslov and Omel'yanov 130 , who analysed Toda lattice with variable coefficients

10.12
They construct soliton solutions in the following way: rapidly changing part of soliton is constructed using Toda lattice with constant coefficients, and for slowly part of solution continuous approximation is used. Then these solutions are matched.

Some Generalization
Let us show that PA for constructing of improved continuous models can be used iteratively. For three terms in expansion 3.9 one has PA ∂ 2 /∂x 2 Now let us consider the continuous models of the chain with two different particles in primitive cell, that is, the chain with alternating masses Continuous approximation for 11.6 is rod with localized masses Improved continuous approximation for 11.6 can be written as follows: y tt x, t .

11.9
Homogenization and other asymptotic approaches can be use for this continuous system 132, Section 3 .

Modified Navier-Stokes Equations
The Navier-Stokes equations have a long and glorious history but remain extremely challenging, for example, the issue of existence of physically reasonable solutions of these equations in 3D case was chosen as one of the seven millennium "million dollar" prize div u 0, k 1, 2, 3.

12.2
Equations 12.1 , 12.2 differ from the classical ones only in large-velocity gradients and coincide with them for ν 1 0. For 12.1 and 12.2 Ladyzhenskaya proved global unique solvability under reasonably wide assumptions.
As it is mentioned in 59 , "From the kinetic point of view these equations represent the long-wavelength limit. However, it often happens that the dynamics predicted by the hydrodynamics equations involves relatively short-wavelength scales. This, formally at least, contradicts the conditions under which the macroscopic system was derived." Note that the modified Navier-Stokes equations 12.1 , 12.2 allow for taking into account higher order harmonics influence via artificially separated terms.
More investigations in direction of improvement of the Navier-Stokes equations are carried out in references 133, 137-139 .
It should be emphasized a role of mathematical approaches in proving various physical theories 141 . Although nowadays a physical experiment plays a crucial role in verification of many novel theories, but it is difficult sometimes to realize a sure and proper experimental verification and validation of the developed theories. In those cases it seems that the proper support of mathematical rigorous approaches may serve as a tool of "experimental" validation of the introduced theories. In other words the rigorous mathematical approaches validating the physical theories can be treated as experimental ones for the so called "theoretical mathematics." Therefore, the modified Navier-Stokes equations are more suitable for fluid dynamics description for large Re from a point of view of theoretical mathematics. For small Re they are more close to the Navier-Stokes equations. In the latter case one gets a unique global solutions to the initial-boundary value problems as well as the solvability of the stationary boundary value problems for arbitrary Re numbers.
Note that the modified Navier-Stokes equations 12.1 , 12.2 allow for taking into account higher order harmonics influence via artificially introduced terms. It is clear that it will be very challenging in getting the modified Navier-Stokes equations starting in modeling from a molecular structure, as it has been demonstrated so far using simple examples of discrete lattice. Here we can mention paper by Rosenau, who applied regularization procedure to the Chapman-Enskog expansion 59 . This approach, however, requires further investigations.

Conclusion
Classical molecular dynamics simulations have become prominent as a tool for elucidating complex physical phenomena, such as solid fracture and plasticity. However, the length and time scales probed using molecular dynamics are still fairly limited. To overcome this problem it is possible to use molecular dynamics only in localized regions in which the atomic-scale dynamics is important, and a continuous simulation method e.g., FEM everywhere else 142 . Then the problem of coupling molecular dynamics and continuum mechanics simulations appears. In this case, one can apply a concept of bridging domain 142-144 . We think that in the latter case an application for describing bridging domain of improved continuous models can result in efficient computational time saving.
Let us compare an impact of the methods illustrated and discussed so far on the improved continuous models. Continuous models based on the PA and two-point PA can be applied in the case of 1D problems in spite of some artificially constructed examples 51 . Intermediate continuous models and BVP obtained on the basis of the composite equations can be used for problems of any dimensions.
It will be very interesting to use for study investigation of 2D lattices 2D PA 145 . Let us finally discuss problems closely connected with brittle fracture of elastic solids 120, 146-150 . Observe that a continuous model, which does not include material structure, is not suitable since any material crack occurs on the material structure level. This is a reason why some of the essential properties of damages of the classical continuous model theory are not exhibited. For instance, in a discrete medium the propagated waves transport part of the energy from the elastic body into the crack edges. This is why application of nonlocal theories during explosion process modeling e.g., for the composite materials looks very promising. As it is mentioned in 151 , "The continuous development of advanced materials ensures that a "one size fits all" approach will no longer serve the engineering community in terms of predicting and preventing fatigue failures and reducing their associated costs."

A. Padé Approximants
Let us consider Padé approximants, which allow us to perform, to some extent, the most natural continuation of the power series. Let us formulate the definition 152 . Let

B. Continuum Limit of Toda Lattice
Here we describe construct of the continuum limit of the Toda lattice m d 2 y n dt 2 a exp b y n−1 − y n − exp b y n − y n 1 .

C. Correspondence between Functions of Discrete Arguments and Approximating Analytical Functions
Let us suppose system mσ jtt t c σ j 1 − 2σ j σ j−1 s j t − s j−1 t , j 0, 1, . . . , n − 1, n C.1 s j t is the external force acting on jth point. Further we follow 44 .
In continuum limit we must construct a function s x,t representimg a continuous approximation of the function of discrete argument s j t . Observe that it is defined with an accuracy to any arbitrary function, which equals zero in nodal points x j, j 0,1,2,. . .,n. For this reason, from a set of interpolating functions one may choose the smoothest function owing to filtration of fast oscillating terms. As it is shown in 44  For h → ∞ function C.4 turns into Dirac's delta function. So for momentous theories of elasticity one must replace Dirac's delta function to the function C.4 .