LIBRARY OF THE UNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGN 510.84 no. 582-587 cop. Z Digitized by the Internet Archive in 2013 http://archive.org/details/energytransportp586miur 676. YH ??u^4 uiuCDCS-R-73-586 THE ENERGY TRANSPORT PROPERTIES OF ONE DIMENSIONAL ANHARMONIC LATTICES By KEN'ICHI MIURA July 20, 1973 ijiucdcs-r-73-586 THE ENERGY TRANSPORT PROPERTIES OF ONE DIMENSIONAL ANHARMONTC LATTICES PY KEN'ICHJ MIURA July 20, 1973 Department of Computer Science University of Illinois at Urbana-Champaign Urbana, Illinois 61801 ♦Submitted Ln partial fulfillment of the requirements for the degree of Doctor of Philosophy in Computer Science in the Graduate College of the University of Illinois at Urbana-Champaign and supported in part by the Department of Computer Science and the Advanced Research Projects Agency of bhe Department of Defense and was monitored by the U.S. Army Research Off J le-Durham under Contract No. DAHC0-U-72-C-0001. THE ENERGY TRANSPORT PROPERTIES OF ONE DIMENSIONAL ANHARMONIC LATTICES BY KEN 'I CHI MIURA B.S., University of Tokyo, 1968 M.S., University of Illinois, 1971 THESIS Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Computer Science in the Graduate College of the University of Illinois at Urbana-Champaign , 1973 Urbana, Illinois THE ENERGY TRANSPORT PROPERTIES OF ONE DIMENSIONAL ANHARMONIC LATTICES Ken'ichi Miura, Ph.D. Department of Computer Science University of Illinois at Urbana-Champaign , 1973 The energy transport properties of the solitary waves are studied both analytically and numerically. It was found that the piecewise linear function model clearlv explains the characteristics of such waves. Numerical experiments were performed extensively to clarify the various properties of the solitary waves in the one dimensional anharmonic lattice These results were used to interpret the heat conduction problem of the crvstal lattice. ill ACKNOWLEDGMENTS The author wishes to express his deepest gratitude to his advisor, Professor E. Atlee Jackson of the Department of Physics and the Coordinated Science Laboratory, University of Illinois, for his constant guidance, encouragement and the most fruitful discussions. He also wishes to thank Professor D. L. Slotnick of the Department of Computer Science and the Center for Advanced Computation, University of Illinois, for letting him pursue this work and for letting him use the computer facilities at the Center for Advanced Computation. This work was done while the author was a research assitant in the Center for Advanced Computation, and its publication was supported by the Department of Computer Science. The author would like to thank Mrs. Judith Rudicil for the most excel- lent job of tvping and Mrs. Gayenne Carpenter for her kind and valuable suggestions in preparing this thesis for publication. Finally, the author is thankful to his wife, Reiko, for her encourage- ment and also for drawing some of the figures. IV TABLE OF CONTENTS CHAPTER Page L. INTRODUCTION 1 1.1 Historical Introduction 1 1.2 Problems of Heat Conduction 5 1.3 Outline of This Thesis 7 2. THE HARMONIC LATTICES 9 2.1 Normal Mode Representation of the Vibration of the Harmonic Lattices 10 2.2 SchrBdinger ' s Solutions 16 2.3 An Analogy between the Harmonic Lattices and Electrical Circuits 21 2.4 Heat Flux and Temperature Calculation for Semi- infinite Harmonic Lattices 21 3. THE ANHARMONIC LATTICES- -CONTINUUM APPROXIMATION AND TODA'S EXACT SOLUTIONS 29 3.1 Polynomial Type Force Function 30 3.2 Quadratic Nonlinearity (^.=0, ^0) 32 3.2.1 • •- 34 3.2.2 £ 39 3.3 Cubic Nonlinearity (^0) 43 3.3.1 h* - (Hard Nonlinearity) 43 3.3.2 h (Soft Nonlinearity) 51 3.4 Piecewise Linear Approximation and Explanation of the Existence of Solitary Waves 58 3.4.1 v 2 • C \ (Subsonic Waves) 62 3.4.2 C Q 2 C Q 2 (! + '-) (Supersonic Wave) . .... 66 3.4.3 Summary 63 3.5 Toda ' s Solution for the Exponential Lattice 70 4. NUMERICAL METHODS AND RELATED TOPICS 7 5 4.1 Numerical Integration of the System of O.D.E 76 4.2 One Step Methods for the Special Equations 78 4.3 Heat Reservoirs 83 4.4 Computer Facilities and Programs 86 4.5 Description of a Typical Program 38 5. DYNAMIC PROPERTIES OF ANHARMONIC LATTICES 90 3.1 Basic Dynamic Properties of JPW Lattices 91 5.1.1 Creation of the Solitary Wave by the Impulse . 91 5.1.2 Energy Concentration 92 5.1.3 Tail Part 97 5.1.4 Mode Spectrum oi Solitary Waves 101 5.2 Boundary Effects ]_gi 5.2.1 Fixed Boundary 10L 5.2.2 Free Boundary 105 5.2.3 Heat Reservoir Type Boundary 105 5.3 Collision Process of Solitary Waves 109 5.3.1 Phase Shift in Two Wave Collision Process . . 109 5.3.2 Energy Loss Ill 5.3.3 Effect of Disturbances 114 5.4 Time Averaged Kinetic Energy 115 5.4.1 Harmonic Lattice (Dispersive Wave Packet) . . 117 5.4.2 Anharmonic Lattice (Solitary Wave Plus Tail)., 117 5.4.3 Temperature Holes 117 5.5 Piecewise Linear Force Function 120 5.6 Hard Sphere plus Harmonic Force Function 120 6. HEAT CONDUCTION PROBLEM IN ANHARMONIC LATTICES AND ITS PHENOMENOLOGICAL INTERPRETATION 127 6.1 Cooling Process of Anharmonic Lattices 129 6.2 Heating Process of Anharmonic Lattices 129 6.3 Heat Flux and Temperature Distribution under Different Boundarv Conditions 132 6.4 Phenomenological Interpretation of Heat Conduction in JPW Lattice 142 6.5 Hard Sphere Plus Harmonic Type Model and Other Possibilities of the Temperature Gradient 144 7. CONCLUSION 146 APPENDIX A. SOME PROPERTIES OF TCHEHYCHEFF ' S ORTHOGONAL FUNCTIONS 148 dw 2 B. COLLECTION OF SOLUTIONS OF I — I = P(w) 151 dz dw ^ C. ON THE TWO SOLUTIONS OF (— ) = P(w) 152 dz D. COMPUTER PROGRAM 154 LIS'l OF REFERENCES 163 VITA I 67 CHAPTER 1. INTRODUCTION 1.1 Historical Introduction The heuristic approach to solve mathematical problems in physics has attracted attention of physicists since the advent of digital computers Undoubtedly, J. von Neumann was the first to propose such approach. In 1944 he solved a simple hydrodynamical problem numerically by using the then available punch-card equipment at the Ballistic Research Laboratory which was developed by IBM for calculating firing and bombing tables. (Note that the first electronic computer ENIAC was yet to be completed at that time!) He solved the one dimensional hydrodynamical equation for shock waves in Lagrangian coordinates: x = F' (x )x tt a y aa where x(a,t) is the position of a moving fluid particle with Lagrangian coordinate a at time t, and F (x ) gives the pressure vs specific volume relation (homentropic equation of state) . His idea was to discretize the spatial coordinate a in order to reduce the above partial differential equation to a set of ordinary differential equations of the form: m d x /dt 2 = f(x -x ) - f(x -x ,). (1.1) n n+1 n n n~l This is nothing but the equation of motion for an array of identical- mass particles connected by springs. (Fig. 1.1) He observed a shock x 2 * wave formed in such mechanical system for f(x) = (1 - —) ,0 < x<2. This was the first numerical experiment ever made on a computing device. He further discussed that the extensive but well planned computational efforts would break the deadlock encountered by purely analytical approach i< Since there is no dissipative term in (1.1), the shock wave which von Neumann observed is essentially a solitary wave which is to be described in the next section. One of Von Nermann's results is rep !>;ced i n Y\ g. I . ? . |- / wM>^R^<>^rM>^H^ -^0?RPO^rjp-0 Fig. 1.1. Mass-spring system. "time 12. 3 * 5" 6 7 8 f tO 11 12 13 11 IS 16J7 i* If 20 U 22 23 2i 25 2Q 27 28 XJ % Fig. 1.2. Shock wave in the systerr (1.1) (2) in such problems as turbulence and shock waves. The next great step was made by Fermi, Pasta and Ulam, who studied the behavior of the vibra- tion in a nonlinear string on MANIAC I computer at Los Alamos scientific (3) laboratory in 1953. Their model was also a discrete one; the equation of motion is identical to (1.1) with the function f given as: 2 (1) f(x) = x+ax (quadratic nonlinearity) 3 (2) f(x) = x+f3x (cubic nonlinearity) (3) f(x) = 6 x jxj < d 6 x+c [xj > d (broken linear function) where a, j3, 6, , 6,, , c and d are given parameters. Their intention was to observe the mixing of energy among the normal modes of the system due to a weak nonlinearity. The result was, however, quite contrary to what they expected: starting with the initial energy given to the lowest mode, only a few modes were excited, after which almost all the energy was returned to the original mode. This result cast doubt on the long-believed hypothesis that only a slight nonlinearity will lead a dynamical system to a statistic- ally equilibrium state. This numerical experiment by Fermi, Pasta and Ulam will be referred to as FPU throughout this paper. The FPU experiment has stimulated extensive research during the last two decades, both computationally and analytically, among physicists and mathematicians. The research directions may be categorized into several areas: (1) ergodic (or non-ergodic) behavior of nonlinear vibrating systems, (2) theory of the normal mode coupling, (3) thermal conductivity in the anharmonic crystals (4) energy transport properties of nonlinear vibrating systems (especially in relation to the continuum approximation) . Oi course these areas are not independent of each other, but closely related. Here, particular attention should be paid to the last area because a special type of solution has been discovered in such dynamical systems; that is, the solitary wave (or the "soliton" as was named by Zabusky and (4) Kruskal ). As a matter of fact the solitary wave has been known in hydrodynamics since 1844, when Scott Russel reported a single elevation of water in a long canal which propagated without changing its shape. Theoretical investigations that followed by Lord Rayleigh, Boussinesq, Kortweg and deVries clarified the nature of such wave; especially Kortweg and deVries derived the equation which is now known as KdV equation: u + uu + Bu = (1.2) t X XXX where u is the velocity of the water. The maintenance of the wave is possible due to the balancing of two different effects: the nonlinearity that tends to sharpen the wave front, and the dispersion that tends to spread out the wave front. Recently the solitary waves have been discovered in plasma waves and ion-acoustic waves as well as in the anharmonic crystals, and much effort has been concentrated on the analysis of such waves. As for the crystal lattices, Zabusky and Kruskal showed that the equation of motion (1.1) may be reduced to KdV equation or its variations (9) by introducing new coordinates, and using an appropriate approximation. Toda found a particular solution of the solitary wave type for the so called "exponential lattice", which is described by a force function f(x) = (l-exp(-bx) ) /b . This is the only exact solution known for any an- harmonic lattice. As for the ergodic problems, many numerical experiments have been done under various conditions. The original FPU situation was re-examined by Tuck for longer periods of time, resulting in the "super recurrence". On the other hand, the system approached to the thermal equilibrium for a (13) strong nonlinearity or for an initial condition given to the higher (14) mode. The phenomenon of FPU recurrence was studied from the viewpoint of the normal mode coupling by Ford, } '^ Jackson^ ' and Pasta e_t ( 18 ) al . Jackson could quantitatively show the behavior of a few lower modes, which is quite in good agreement with the result obtained by FPU. The development of the area of the heat conduction will be dis- cussed more in detail in the next section. 1.2 Problems of Heat Conduction Among several published reports on numerical experiments of the heat conduction problem, the one by Jackson, Pasta and Waters (to be referred to as JPW throughout this thesis) is very well documented on the computational aspects of their model too. Therefore we will de- scribe the problem along the lines of their report. Let us consider a one dimensional lattice which consists of N atoms of identical mass m connected by anharmonic springs. We choose free boundary conditions on both ends. The force function is of polynomial type: f(x) = x - Xx 2 + "x 3 (1.3) where A and H are the nonlinear coefficients for the quadratic and cubic 2 k 1 terms, respectively. The special relation h = 2 \ /3 or h = — is used to assure that the force is always attractive. X may be varied from one spring to another (JPW called this "defects"), but we will only be concerned with a constant \ in all cases. Next we define a heat reservoir in the I'ollowing way; every time an end particle hits the reservoir, it loses its velocity and is kicked back toward the lattice with a new velocity v sampled from the probability density function p(v) = (v/T)exp(-v 2 /2T). (1.4) This definition was given by JPW and is very well suited for numerical experiments on digital computers. The lattice is located between two heat reservoirs of given temperatures and is subject to random motion. If two temperatures are different we expect an energy flow from the high tempera- ture side to the low temperature side. In other words we can define the heat flux J numerically by the time-avaeraged energy transferred through the lattice. We can also define the temperature of each atom, either by time-averaged kinetic energy or by time-averaged total energy. Usually the former is chosen for its simplicity. Fourier's law of heat conduction states that the heat flux J is proportional to local temperature gradient so that the thermal conductivity can be defined as follows: J =-hVT. Computationally, ■" can be calculated from the averaged heat flux J and the linear-interpolated temperature distribution. JPW performed several types of numerical experiments on the lattice-reservoir system just described. As they mentioned in the general description of their results, pulse-like waves are observed to propagate through the lattice without much distortion in spite of the irregular motions of atoms and defects (Fig. 2 and Fig. 3 of (21)). This fact suggests that the phenomenon of heat conduction might be considered in terms of localized waves (such as those we introduced in jfl.l). This leads to such questions as: (1) Is it possible to create solitary waves by impulse type initial conditions ? (2) If so, do such waves propagate through the lattice without appre- ciable dispersion of energy? (3) How do such waves contribute to energy transfer? (4) What are the effects of boundaries? (5) What will be the possible sources of the temperature gradient for anharmonic lattices which was observed by JPW and other, authors? And so on. Our main purpose in this thesis is to establish a dynamical picture of the heat conduction through the lattice in terms of the localized waves and also to find the limitation ot such approach. We will also dis- cuss some of the related topics herein. 1.3 Outline of This Thesis This thesis consists of two parts, one being on the mathematical derivations and the other being on the numerical experiments. In Chapter 2 we will study the dynamical properties of the harmonic lattices. An emphasis will be put on the solutions discovered by Schrodinger, which are very convenient for describing localized solutions. Some of the thermal properties of the harmonic lattice will be also described. In Chapter 3 we will discuss the properties of the anharmonic lattices; first in con- tinuum approximation and then about Toda ' s exponential lattice. The soli- tary wave solutions and the periodic wave solutions are obtained in both cases. The special features of such solutions will be extensively discussed and they will be simply explained in terms of elementary functions. Chapter l \ is about the numerical methods and some related problems encountered in the actual computation. We will show that Nystrom's integration formula for the special equations of the second order is quite advantageous for our purpose. In Chapter 5 we will analyse the computational results on the basic dynamical properties of our model lattice. Specifically we will show that the solitary waves can be created when the lattice interacts with the heat reservoirs. Since no exact solution is known for our model lattice, the results obtained in Chapter 3 are applied whenever applicable. In Chapter 6 we will discuss the heat conduction problem, considering the results in Chapter 5, our own computational results in the thermal situa- tion as well as those obtained by JPW and other authors. Some of the questions raised in 21.2 will be answered. CHAPTER 2. THE HARMONIC CRYSTAL LATTICES We consider the one dimensional crystal lattices in this chapter which consist of either an infinite or finite number of atoms of identical mass and connecting springs which obey Hook's law, f(x) = Kx. The equation of motion lor the displacement y of the nth a tom from its equilibrium position is then d 2 y m dt Boundary conditions and initial conditions are to be described later. 4 The linear theory of lattice vibration has a very long history in the theoretical physics, or we should rather say that a considerable part of mathematical physics was developed around this very subject in its early stages. Newton was the first to investigate this problem. He speculated that "the phenomena of Nature ••• -may all depend on certain forces by which the particles of bodies* • -are either mutually impelled towards one another, and cohere in regular figures or are repelled..." (Principia). it is along this line that Newton considered how the sound wave propagates in the air. By discretizing the continuous media such as air into particles of mass m he obtained the equation of motion identical to (2.1). Note that the theory of the partial differential equations was not developed at this time. Newton's idea was developed by Bernoulli, Lagrange, Cauchy, Hamilton and many other well known mathematical-physicists. For more details on the historical development of the theory, we refer the (22 ) reader to L. Brillouin. The basic idea behind the development of the theory oi harmonic lattices is the decomposition of the vibratory motion into the Fourier components, or the normal modes of the system. 10 2.1 Normal Mode Representation of the Vibration of the Harmonic Lattices Suppose that we have a finite number of atoms. The normal mode representation is obtained by putting the solution of (2.1) in the follow i ng form: y (t) = f exp(iuut). Substituting (2.2) into (2.1) yields a finite difference equation for f : n f J.1-C2- T^ 2 ) f +f i = ° > n+1 K n n-1 /K which we can solve bv putting f = exp(ikn). By defining U) = 2*/ — we n o jij obtain the dispersion relation x = uu sin (k/2) •_'■ k < n . (2.2) The system is dispersive for large k. The values of possible k's depend on the boundary condition, but in general there are N distinct k's (we call them the modes of the system) if there are N atoms in the system. The larger the number of atoms in the system, the greater the density of the allowed values along the k axis. Note that uu is always less than JJ . o This is why the mass-spring system is often called the mechanical low-pass filter. Actually we can show that there is a close correspondence between the mechanical filters and the electrical filters (section 2.3). Also refer to L. Brillouin. (22) Going back to the problem of finding the k's for a given boundary (23 ) condition, we take the following steps : (i) Put f = A exp(ikn) + B exp(-ikn) and find k's that satisfy the boundary condition, (ii) Express B (or A) in terms of A (or B) to obtain a unique repre- sentation for f . This gives a set of orthogonal functions from n ° ° 11 which we can obtain a set of orthonormal functions, (iii) Obtain w for each k. They are all located on the sine curve given in (2.2). Some typical boundary conditions are listed in Table 2.1 together with the values of k and uu, and normal coordinates. There is another type of solution for an infinite lattice which is obtained by putting y n (t) = f h -exp(-2(3t) (2.3) where f n - exp(2an). (2.4) Following the same procedures as before, we can obtain the relation between a. and j3 as: £3 = U) sinh (a) . (2.5) Because of its unboundedness this solution is seldom mentioned. If the amplitude growth is suppressed by some nonlinear effect, however, the resulting solution will have an exponential behavior in the region in which the amplitude is sufficiently small. This is actually the case for Toda's solitary wave solution. We will come back to this problem in §3.5. Going back to the normal mode representation, we can find the solution of (2.1) for a given initial condition by superposing the normal modes. As an example, let us consider the wave motion when the left end atom of a free-ended lattice is given an impulse of unit magnitude at t=0. The general solution is: p. B y (t) = A +B t + N [A cos(w t)+ —cos (a> t ) ]cos[mn(2n+l) /2N] (2.6) n oo m mOJ m <-- > m m=l Table 2.1 Normal Mode Representation periodic 12 y (t) = a (t)exp(2nimn/N) where a (t) = exp(id) t) m m 2 K 2 and a) = —sin (TTm/2N) m m m = 1,2, ■ • • ,N free ends l (t) = a (t)cos[Ti m (2n+l)/2NJ m - 5 1,---,N-1 where a (t) = exp(iJj t) m m ' 2 K 2 and a, = —sin (TT m /2N) m m O y * 5 K>TnfK>wK>w x iD fixed ends V (t) = a (t)sin( r imn/N) ' n m m = 1,2, • • • ,N-1 where a (t) = exp(iO- t) m r m 9 K 2 and ^!I^O^^<>^^> yTV t 13 where A and B are the Fourier coefficients of the initial displacement and velocity which will be given shortly. We can also obtain y (t) as ^n follows : N-l Y n (t) = B q + V ) [B m cos(uJ m t)-A m iU m sin(u) m t)]cos[mn(2n+l)/2N]. (2.7) m=l Initial conditions are N-l \ y (0) = A + ;■ A cos[mrr(2n+l)/2N] (2.8) n o m v ' m=l 3nd N-l y (0) = B + ; B cos[mH(2n+l) /2N] . (2.9) n o m m=l Or solving for A and B , m m' N-l e , -, A = -7 S y (0)cos[mn(2n+l)/2N] (2.10) m N n n=0 B = m N-l e m N y (0)cos[mn(2n+l)/2N] (2.11) n=0 , . 1 i f m=0 . , _ . where E =- . . A and B give the motion of the center of mass of m 2 if m?t0 o o ° the lattice. For the present problem, we have the following initial condi tions : y =0 for all n — -rr— ;■ A = for all m, n m E y = 6 n > B = -7cos(mn/2N) . ^n n,0 / m N Substituting A and B into (2.7) yields m m N-l y (t) = 1/N + 2/N cos(mn/2N)cosLtsin(mn/2N)]cos[m(2n+ln/2N)]. (2.12) m=l 14 This is the solution for the given initial condition. Note that the amplitude distribution among the normal modes is not uniform but shows a monotonous decrease as the mode number increases. The energy distribution is obtained from the square of the amplitude distribution, since initially there is no potential energy and the total energy of each mode is strictly conserved. Both distributions are shown in Fig. 2.1 and Fig. 2.2. For the anharmonic lattices with the same boundary and initial condition we can also use this distribution approximately if time is very close to t=0. This will be discussed in §5.1. Next let us consider the case when N goes to infinity. From (2.12) we have 4 > n/2 v (t) = — cos(x)cos[t sin(x) ]cos[ (2n+l)xJ dx +0(1/N) - n TT J o (N -» °°) where where mrr n X = 2^ ' dx = 2^ ' 4 ^ = ~ J cos(wt) cosL (2n+l )arcsin(w) J dw w = sin x, dw = cosx dx , 1 = "| j cos(wt)(-l) n U 2n+1 (w) dw, (2.13) where U 1 (w) is the second kind of solutions of Tchebychef f ' s equation: 2 2 2 2 (1-w )d U (w)/dw -wdU /dw+n U (w) - 0. Making use of the following identity n n n for an odd n, . . , , , o 1 , U (w) J (t) - j (-i) ~ cos(wt) dw = — - — , (A. 7) o See Appendix A for the derivation. AmpMuJe l"-- r «* 15 1 / X 3 f r ' 7 * -m. Fig. 2.1. Amplitude distribution (N=9) Ener^ i >— - 1X1*5* 7 t m Fig. 2.2. Energy distribution (N=9) 16 we can obtain the expression for y (t) when N J 2n + l (t) y (t) = 2(2n+l) . (2.14) n t We will obtain this solution much more easily in the next section. This solution is the response of the velocity of the n ttri atom to a unit impulse at the left end.; it shows how the energy propagates through the harmonic lattice. In a thermal situation initial impulses are given at random. We will come back to this problem in §2.4. 2.2 Schrodinger ' s Solutions A different type of solutions was found by Schrodinger in 1914 (24) for infinite harmonic lattices. Since his solutions are very convenient for describing the energy transport through the lattice we will summarize some of their features. Let us introduce the new coordinates T j according to Schrodinger; TL - v'm y (2.15a) 2n n \m - ^ 0, we will restrict ourselves to such solutions. Therefore the general solution is given as a linear combination of Bessel functions as follows: k where a, 's are determined from the initial conditions. Let us consider k an initial impulse given to the atom when all the atoms are in their equilibrium positions. Eventually we will consider the example discussed in §2.1 in terms of Schrodinger ' s solutions. We have, for the present example, a. = vm &, ^ , or k k,0 T] = Vm J(f) n n Going back to the original coordinates, we have y„ - J o ^o^ (2 - L9) n 2n o V y n+1 ■ '' WK J 2n+l ( V> -i" J 2»rt- (2 - 20) O In general Schrodinger ' s solutions are very convenient for studying the propagation of local disturbances through the lattices. Fig. 2.3 shows how an initially localized disturbance propagates with time. As is clear from Fig. 2.3 , the harmonic lattice is very dispersive and energy spreads out rather quickly. From the asymptotic formula for Bessel functions, the , ,"1/2 (25) amplitude of the wave is known to decay as (w t; In order to satisfy the free boundary condition at the left end of the lattice, we add another impulse of the same magnitude to the atom 18 ■3T» ■ ' T si J— « — •— » — ' — ' — "- <*> SO ' i t £ 1_J I I i I I '" f I" o •r-l c o B i-t pci (u ■u 0) cti J2 too •.-4 fa ■3^ I ' ' o to i vo I , Q 1 H OJ Ml C •r-l SO u si u CO 00 fa 19 on the left of atom n=0. Then the solution is *n = J 2n ( V } + "W^o * (2 ' 21) In particular the atom n=0 and n=-l oscillate exactly in the same way. Initially there is no energy stored in the spring that connects them; therefore atom n=0 actually feels no forces from its left neighbor. Con- sequently we can eliminate all the atoms with negative indices, and the resulting lattice is semi-infinite. By making use of a recurrence formula < J n-l +J n + l )/2 ' " V X ' < 2 ' 22 > (2.21) may be written as y = 2(2n+l) J_ J _. (cu t)/w t . (2.23) n zn+1 o o We can also obtain the response of y -y , to a unit impulse in the same r J n n+1 r way : y -y ,, = 2/oj [j_ . (J t) + J. -(cu t)] y n ^n+1 ' o 2n+l v o 2n+3 v o = 8(n+l)/x [j (a, t)/0) t] . (2.24) o 2n+2 o o Fig. 2.3 shows how each atom responds to an impulse; we can see that each wave decays somewhat faster than in the previous example. This is due to the fact that there is no force to pull atom n=0 back to the original -3/2 position. (2.24) asymptotically behaves as (u> t) . In the similar way we can obtain the velocity response function when the initial impulse is applied to the k th atom. The results are *„«> " J a(k-n> ( V> + J 2(n +k ) + 2 ( V> 0 n = k ^n (t) = ^(n-k)^ + ^(n+WV* " >k ' It is easy to see that the free boundary condition is satisfied at n=0. Parenthetically, the solution we just obtained is also a solution to the infinitely long lattice with atom n=0 twice as heavy as the other atoms, to which an impulse is applied. Schrodinger ' s solutions are very simple and clear as long as the lattice is infinite, but in the case of a finite lattice the reflection of waves must be taken into account. Consequently we have an infinite series of Bessel functions instead of a single term, even for the response func- tion to an impulse. For example, the velocity impulse function for a lattice with N atoms (n=0 and n=n-l being free ends) is: F (oj t) - J (U) t)+J. ,„0 t)+J (ou t)+J. M (uj t) n o 2n o 2n+2 o 4N-2-2n o 4N-2n o + ^N^ntV^WWV^SN^^n^o +J QM (a) t) 8N-2n N o + = ^ J [j 4kN+2n (u, o t)+J 4kN+2n+2 (Ui o t)+J 4(k+l)N-2n^o t) k=0 + J 4(k + 1 -2 ] - (2 ' 25) We can easily show that the free boundary conditions are satisfied (the method of mirror images). The response function (2.25) is identical to the 21 solution which we obtained in §2.1, (2.12). As we can see from this example, the normal mode representation and Schrodinger ' s solutions are complimentary to each other. 2.3 An Analogy between the Harmonic Lattices and Electrical Circuits There is a one-to-one correspondence between the harmonic lattice and electrical low pass filter (L-C ladder circuit) as has been (22) discussed by many authors, for example L. Brillouin. It is often very convenient to consider an equivalent electrical circuit of a dynamical system; especially the notion of impedance in the former enables us to skip some of the lengthy calculation. As an example we will consider Schrodinger ' s solution here. Table 2.2 summarizes the correspondence. The normal modes of an electric circuit such as a low pass filter may be obtained by looking at the zeroes of the impedance. Hence the normal modes of the corresponding dynamical system are easily obtained by simple substi- tution. The impedance is also very useful when we consider the wave reflection at an impurity in the harmonic lattices. This is a problem of impedance matching in electrical circuits. But we will not go any further ( 7 ft ^ in detail on this subject. Refer to van der Pal 2.4 Heat Flux and Temperature Calculation for Semi-infinite Harmonic Lattices As a simple application of Schrodinger ' s solutions we calculate the heat flux and the temperature distribution of the semi-infinite har- monic lattices which we discussed in §2.2. We suppose that only the atom at the end of the lattice interacts with the heat reservoir; it receives random impulses from the reservoir. The approach we take here follows (27) ol niilenbeck and Ornstein in their classical work on Brownian motion. Table 2.2 22 Dynamical System Electric Circuit difference of the displace- \ V : voltage of the n capacitor merit of two neighboring atoms (i.e. , y n _ 1 Y n ) S = 1 dy n ' : velocity of n*-" atom I K dt J ; n current through n th inductance m: mass of the atom L: inductance k: spring constant boundary free end on the left C : (reciprocal of) capacitance short circuit on the left a) = 2v'K/m o 1 0) = 2/ v / LC o initial condition y o (t) = V Q /2K6(t) ! V = E 6(t) ( o o equations 1/Kdu /dt = S . -S n n-1 n m dS /dt = u -u n n n n+1 solution S (t) = v /KJ_ (uj t) n o 2n o CdV /dt =1 -I ! n n-1 n Ldl /dt = V -V , n n n n+1 I (t) I n 2E /L J_ (oj t) o 2n o or y n^ t} = V o J 2n Table 2.2 (continued) 23 impedance .Uj Z(io) = v mK J L-(— ) o Z(uu) = c ?o 2 J L + 1 , = A/uu ds o z 25 4A(2n+l) 2 r(2n+l/2) (^ o 2 2 2 r(3/2)r(3/2)r(2n+l+3/2) A(2n+l) 2 /ra 2 o = 2 , (2.27) (2n+l) -1/4 where the phase average is taken at time t over an ensemble of lattices (a large number of similar but independent lattices), which have started at t=- oc with identical initial conditions given to the corresponding atoms (25) therein, and we have also made use of the following formula J J m ( ° V° m Y(\) ((m+n-W/2) ° t 2 K r((\4tn-n+l)/2) r((\+n-m+l)/2) f( (m+n+X+1) /2) (2.28) We can determine the constant A from the condition y (t) = kT/m, that is, 2 A = 3hju kT/4m or o f( T )f(T') = 3thju kT/4m6(i-T') . o Fig. 2.4 shows the temperature distribution. There is a slight decrease of the temperature near the end atom, then the temperature becomes almost uniform. This tendency is not changed if we define the temperature of each atom by the total energy instead of kinetic energy. The heat flux (instantaneous flux at each point of the lattice) is obtained by using the (21) following formula ; J H (n) - -K(y n +y n+1 >(y n+1 -V/ 2 ■ (2 ' 29) We can obtain the first term immediately from (2.23): ( VW /2 ■ f f ('» J 2„ ( "o t)+2J 2n + 2 ( "o t)+J 2„rt (a, o t))dl/2 - (2 - 30) 26 it., m n Fig. 2.4. Temperature distribution for semi-infinite lattice. 27 To calculate the second term we have to define the response function G (t) n for y ,,-y . From (2.24) we define y n+l n G O t) = - N ^ """• — - — . (2 3D no U) a) t ^ . ~ l ; o o ((D t) = . 8(n + l) ^n+Z^o^ O 'X 1 O) t O o Unlike F (u) t) this function has the dimension of time. Using G (a t) we n ° no can obtain the second term as a convolution form: t Vn^o " y n ( V } = J" f( T )G n (^(t-T))dT. (2.32) — OC Substituting (2.30) and (2.32) into (2.29) yields m Slaa3*(f f(T)(J 2n («» o t) + 2J 2n+2 («. o t) + J 2 ( . t))dT) o _ot> X Cl f(^> ' ( °_ T) dT') (2.33) where we have neglected the effect of initial conditions as before. Taking the phase average and changing the variable x (t-T)=s yield ,, ., v ... a J (s)J ^(s) 2 J. , (s) J ^ (s)J ^ (s) ~ — 7~"T 4 (n-fl)aAK p , 2n 2n+2 2n+2 2n+2 2n+4 J u (n) = — ^ 1 ( ) + + ds H 3 " s s s ju o ° 2 _ 6nkT(n+l)K • 2n+2 ^' ■ ds ma J s o o 3nkTK 2 ma o 3nkTai o 8 (2.34) where we have made use of the formula (2.28) for \=l , or more precisely, " J m ( ° J n (t) 2 sin((m-n)n/2) A . , dt = *"r t 1 — ' — =0 if m-n: even o (m -n ) tt and 28 « J 2 (t) ^ m At- —L J ~1 dt = 2m • o Hence a constant heat flux exists. Since the temperature distribution becomes more and more uniform as n increases, it is obvious that Fourier's law of heat conduction does not hold in this case: we cannot define the thermal conductivity! Similar results have been obtained by several authors /TO \ /O Q \ for different boundary conditions and different thermal situations. The results obtained in this section entirely depends on the fact CO r 2 that the integral ( F (t)dt converges. That is, the energy fed by the o random force to the left atom is transferred to other atoms at a proper rate. As for classical problem of the Brownian particle (Langevan equation) the energy fed from the random force is dissipated by the viscosity so that there is no accumulation of energy in the Brownian particle. If we consider the Schrodinger ' s lattice, for which the velocity response func- f* h tion is given as F (t) = J (cu t) (suppose that the atom interacts with the heat reservoir), the temperature of each atom rises indefinitely since 00 r 2 the integral ( J_ (^ t)dt diverges for all n. Such catastrophe also o occurs for the harmonic lattice of finite length if there is no dissipation of energy. 29 CHAPTER 3. THE ANHARMONIC LATTICES --CONTINUUM APPROXIMATION AND TODA'S EXACT SOLUTIONS We can approximate the system of nonlinear ordinary differential equations (1.1) by a nonlinear partial differential equation, provided that the relative displacement of any two neighboring atoms is not large. In the first part of this chapter we will treat such approximation and summarize two special solutions. These solutions have their complete analogues in the water wave theory , in which they are called the cnoidal wave (periodic solution) and the solitary wave (localized solution). ' A new name "Soliton" is frequently used for the latter because of its (4 ) particle-like behavior (due to Zabusky ). We start from the equation of motion ,2 y„ = f(y^,-y J - f(y -y„_,> • (3-D , 2 J n VJ n+l 'n' v/ n 'n-1 dt We assume that f(x) ~ x when x ~ 0, that is, when the lattice is near its natural length. Historically, Zabusky and Kruskal were the first to introduce ( 9 ) the continuum approximation to (3.1). They reduced (3.1) to a partial differential equation for y(x,t), which is equivalent to (3.5). Then they reduced it into the KdV equation by considering the waves which propagate in one direction (say, in the positive direction). Since the equation (3.5) has solutions which are similar to those in the KdV equa- tion and (3.5) can also describe the waves which travel in the opposite direction, we will use (3.5) as our basic equation. Recently, Wadati and Toda showed that the two solitary wave solutions of both equations give a very good agreement. Readers are referred to (30) and references therein for the KdV equation. 30 In §3.1, §3.2 and §3.3, we will consider the polynomial type force function, specifically the quadratic nonlinearity and the cubic non- linearity, for which we can find the analytic solutions. We will show that both the periodic wave and the solitary wave can exist. For the heat con- duction problem, the periodic solution is quite irrelevant. But it will also be treated in full since it seems to be the fundamental solution, of which the solitary wave is a special case. For example, only the periodic wave can exist as a travelling wave solution under certain conditions. In §3.4 we will simplify the force function by introducing the piecewise-linear approximation. We will show that the essential features of the solutions obtained in §3.2 and §3.3 are still retained in this approximation. By this model we can explain why and when the solitary wave solution is possible. The contents of §3.4 has not been discussed before. In §3.5 we will summarize Toda ' s solutions for the exponential lattice. 3.1 Polynomial Type Force Function In the following sections we will consider a cubic polynomial f (x) = x + ax 2 + £x 3 (3.2) for the force function. It is more convenient to introduce the relative displacement of two neighboring atoms u = v , -y instead of y itself. ° n - n+1 J n n Then the equation of motion for u is obtained from (3.1) as H 2 -^u = f(u ..) - 2f(u ) + f(u ) . (3.3) , I n n+1 n n-1 at We will denote the spatial coordinate by x in the continuum approximation. Then we have ,u/ du \ h i ° u n , " / o u N n ,ou, Vi = u x=n± h ± T^T$ + 2T ( 74 } + ' • • — x=n ox x=n ox x=n ox x=n ,_ , N (3.4) 31 where h is the interatomic distance in equilibrium. Substituting (3.4) into (3.3) with (3.2) yields u = h u" + 2oh 2 (uu 1, +u* 2 ) + 3ph 2 (u 2 u"+2uu* 2 ) .4 + -u«" + ••• where • means the time derivative and ' means the spatial derivative. Putting h 2 = C 2 , t :C 2 = 2ah 2 , 36h 2 = p-C 2 and h = £- = -2- , o o ^ o 12 12 ' u = C 2 [u" + e(uu')' + u.(u 2 u')'] + hu ,,m . (3.5) In order to consider the travelling wave solution, we put u(x,t) = w(x-vt) - w(z) where v is the propagation velocity of the wave (the phase velocity in case of periodic solution) which is to be determined in relation to the amplitude and other factors. Then (3.5) yields 2 2 ? ? *C u.C (v -C )w" = hw"" + — 2_ ( w z )m + — 2_ (w j )m _ o L J Here ' means the derivative with respect to z. Integrating both hand sides twice yields 2 2 2 2 • ■ o 3 o 2 , o C D' w" = - -r w - — w + w - — z - — (3.6) 3h 2 'i ^ K K v where C and D' are the integration constants. We eliminate the term C -— z since it contains z and is regarded as a spurious term which appeared in the process of integration. This can also be verified by approximately (3.1) by a nonlinear partial differential equation for y and comparing its integrated form (in which w corresponds to -tp(z)) with (3.6). So we adopt (3.6) with C = 0. Multiplying (3.6) by w 1 and integrating with reg.ard to 32 z once yield 2 2 2 2 9 hlC EC v -C . , , N 2 o 4 o J o 2 „ (w } = " ~eZ~ w " T5T w + — ~ w + Dw + E (3 - 7) 2D' where D = - and E = 2E' are the new integration constants. The right hand side of (3.7) will be denoted by P(w). We will find the solutions of (3.7) in the following sections. An emphasis will be put, whenever the periodic solution is obtained, on the wave for which the average over its period vanishes (w = 0) . Such wave motion takes place around the equilibrium positions of atoms when the crystal lattice is of natural length. If the external force is applied to the lattice so that it is either compressed or expanded, a periodic wave is set on around the new equilibrium length. This corresponds to considering a solution for which w ^ 0; w > for the expanded lattice (to be called "positive background") and w < for the compressed lattice ("negative background") . The propagation velocity of the waves under such circumstances differs from that of the waves in the lattice of natural length due to the nonlinearity of the force function. As for the solitary wave solution, we are most interested in the one for which the lattice is of natural length at z = + °°. This means that the energy is concentrated at the wave front. (Fig. 3.1) We call such wave zero-back ground wave. 3.2 Quadratic Nonlinearity (u.=0, e^O) The equation (3.8) yields, in this case 2 2 2 dw 2 KC o 3 V o 2 (^) = " -3^- w + -f- w + Dw + E = P(w). (3.8) 33 equi'libn'm I \ x ^ \ \ \ \ \ \ ' \ \ \ v \ N v V \ \ \ \ \ \ \ \ — ^ 30 h tatv u/q. ve. Fig. 3.1. Solitary wave on zero background, 34 3.2.1 e > Periodic Wave We consider the case in which P(w) has three distinct real roots: 2 P(w) = P (a-w) (w-b) (w-c) (a > b > c) 2 2 2 r, EC ^ _ V "C L I o 2 , , N o p =— ,P (a+b+c) — . Then we have a bounded solution in [b,a] (Fig. 3.2) which may be expressed in terms of the Jacobian elliptic functions as where w(z) = a-(a-b)sn \~~^ — Pz,k) = c+(a-c)dn (^— ^ — pz,k) 9 o — K •>'•■ where k = is the modulus of the elliptic function. Going back to a-c r ° x,t coordinates, we have 2 2K u(x-vt) a a-Asn ( — (k x-u)t),k) (3.9a) = c+ ^ dn 2 (^(k x x-cut),k) (3.9b) k where and A = a-b (amplitude TTC ' 3h k x "" 4Kk v (wave number) u) = vk x (frequency) K = K(k) (the complete first kind) . The velocity v is given in terms of a,b,c from the quadratic term in (3.8): 2 2 2 V o 2 o jf- = P (a+b+c) = —^~ (a+b+c) or v 2 = C o 2 (l + j (a+b+c)) . (3.10) See Appendix B for collections of solutions 35 . 3.2. Quadratic nonlineari ty (e > 0) : periodic wave 36 Note that v > C if a+b+c > and vice versa. The periodic solution is o completely determined by three zeroes of P(w). If we further assume the zero background wave, we have the two parameter family solutions, for which A, k , uo and k are explicitly related. First we should note that -^ J sn 2 xdx = -^ J -|(l-dn 2 x)dx = -|(1 - |^) 2K 2K J . 2 V y . 2 ok k K(k) n 2 E (k ) since Z(x) = Jdxdn x - ^ x is periodic of period 2K (Zeta elliptic func- o tion of Jacobi) . Here E(k) is the complete elliptic integral of the second kind. Using the above relation, we have or Therefore u = a-A sn = a- — r(l- — ) = k 2 K A n E ^ a = -2<1- -) u(x-vt) = — (dn (— (k x x-u)t)) - -) k A_ d_ , 2 dZ Z(z) z= — (k x-a)t) rr x From the definition of A and k , we also have b = a-A c = a- 2 " Substituting a,b,c into (3.10) yields 2 ? e A v = C Z (l+ ^(3a-A- ^r)) k = c 2 (1+ 1(3A _ 3AE _ o U+ 3\2 ,2 V A k k K r 1 (^+ £A a i 2 3E = C o (1+ „2 (2_k " K )} k 2>) (3.11) 37 Hence the velocity of the periodic waves can become either "supersonic" (v > C ) or "subsonic" (v < C ) depending on the degree of deformation o o (i.e. the modulus k) . The critical modulus k is given by solving 3E(k) 2 ■ ' ), ( = 2 - k , for which v = C holds independent of the amplitude A. Approxi' K(k) o r r ^ mately k ^ . 98. J c The dispersion relation for zero background wave is 2 2, 2 „ 2 ,, £ N1 2 AK 2 ,3E _ N1 4 uu=vk = C (1- -A)k - k(— -) (— - -2)k x o 3 x n /v K x A 2 2 2 4 If k » and — r finite, 'Jj = C k -4nk (dispersive harmonic wave) and 2 o x x r if k > 1 uu = C k (H — r) (approaches to solitary wave). o x 3 Solitary Wave If b approaches to c, we have the solitary wave solution. Since k = • 1 , dn + sech . In this limit (3.9b) reduces to a-c 2 u(x-vt) = b +A sech (ctx-pt) (3.12a) where A = a-b v'a - b . / >--A a = J 3* V 12h v 2 = C q 2 (1+ ^(a+2b)) (3.12b) and 3 = u*v . Consider the zero-background solitary wave, i.e. b = 0. Then 2 u(x-vt) = Asech (ax-pt) (3.13a) where v 2 = C 2 (1+ |A) (3.13b) o J and \i = a* v . Note that the solitary wave is always created toward the positive w direc- i Lon. Since u = y , "V > 0, this means that only rarefactive solitary n ^n+1 J n wave is possible for t > 0. 38 Fig. 3.3. Quadratic nonlinearity (£ > 0) : solitary wave Fig. 3.4. A typical profile of solitary wave. 39 3.2.2 e < The derivation of solutions for the negative £ is almost identi- cal to 3.2.1. The equation (3.7) yields (-j^O = P (a-w) (b-w) (w-c) where o K|C v -C Z 2 ' o 2 . ,, . o P = -3^— , P (a+b+c) = — . Obviously bounded solutions exist only in [c,b] (Fig. 3.5). The periodic solution is: w(z) = c+(b-c)sn (r~^ Pz,k) 2 \I 3. — c = a-(a-c)dn (f— - — Pz,k) or . 2 b-c k = a-c 2 2K u(x-vt) = c+Asn ( — (k x-wt),k) (3.14) = c+ ^rfl-dn(%k x-0)t),k)3 where k 2 v n ,~ x - w,-/ j A = b-c (ampli tude) k 2. x 4kK sj /■$* (wave number) co = vk X (frequency) K = K(k) and v 2 = C 2 (1- -Lfi (a+b+c)) . (3.15) o J Obviously the velocity becomes "supersonic" if a+b+c < 0, and vice versa. We can obtain the more explicit (two parameter) solutions for zero background waves (u=0) by following the same procedure as before. 40 Fig. 3.5. Quadratic nonlinearity (& < 0) : periodic wave, 41 Results are is u(x-vt) = - — {dn (— (k x x-U)t),k)- |) k 3k 2 K The critical modulus k is also " .98. The dispersion relation c 2 _ 2, 2 .4K/,. 2, 3E ox , 4 a- = C k - h( — ) (k + — - -2)k ox v tt /v Kx Hence A 2 2 2 4 if k -► and -r- finite, jj = C k -4nk and ,2 o x x k if k - 1 2 n 2,.^ 'EJA N , 2 jo = C (1+ — )k o 3 x Solitary Waves We also have the solitary wave solution in the limit a ->b (Fig. 3.6) u(x-vt) = b-Asech (ax-|3t) (3.16) where A = b-c a = yiil and P = < i • v 2 2 v = C (1- o _llA (2b+c)) (3.17) For the zero-background wave, (b=0) u(x-vt) = -Asech (ax-; it I (3.18) where A = -c > /HE 42 pc«» Fig. 3.6. Quadratic nonlinearity (e < 0) : solitary wave and o 3 43 v2 " C , 2 ( 1+ -T A )- (3-19) Only the compressive solitary wave is possible if £ < 0. Com- bining this result with the one obtained for £ > 0, we can summarize that the solitary wave is always created in the region in which the nonlinearity becomes hard . This is one of the common features of the solitary waves. 3.3 Cubic Nonlinearity (u- ^ 0) We can expect that the higher order nonlinearity will produce more variety in the solutions. We often put £=0 in this section for the sake of simplicity, because the quadratic nonlinearity does not change the nature of the solution. That is, the cubic term in (3.7) can be always absorbed in the quartic term by a simple transformation of variable: 2 2 2 2 ^o 4 ' C o 3 V " C o 2 P(w) = - — w - w + w + Dw + E x ' 6h 3h H 2 2 P.C 4 C 2 2 2 o + D* (w+ ^— ) + E* l\i where D 1 , E' are new constants. It is clear that the existence of the cubic term results in the shifting of the solution and of the velocity. 3.3.1 u- > (Hard Nonlinearity) ? 2 2 ^ C o Let P(w) = -P (w-a) (w-b) (w-c) (w-d) where p = . a,b,c,d are the four real roots of P(w) = and they satisfy the following equations: « a + b + c + d = — ^ = - — (3.20) 3kp 2 * 44 2 2 v -C 2 ab+bc+ca+ad+bd+cd = -j_ = - ~(~r - 1). (3.21) Kp 2 ^ C 2 o Periodic Wave There are two bounded regions in w, for which the periodic solu- tions exist. They are [d,c] and [b,a] (Fig. 3.7). The periodic solutions are expressed in terms of the Jacobian elliptic functions (see Appendix B) or and or where and a c_d 2 / y(a-c)(b-d) d + a sn (-^ f-* L pz,k) / \ a - c z w (z) = 1 + £lA jti 2 ^/(a-c)(b-d) 2 PZ ' k) d < w < c , , c-d 2 .2K,. . . . d + a sn ( — (k x-uut) ,k) a-c tt x "l (X " Vt) " , 4. <=-<> V*a ^T^ ' (3 - 22) 1 + sn ( — (k x-uot),k) a-c rr x j / a "b s 2 ^(a-c) (b-d) , . a + d(— ) sn (^ ^ *- pz,k) w (z) = ~ b < w < a a-b 2 , V (a-c) (b-d) , s 1 + — sn (^ ^ *- pz,k) a + d(^) sn 2 (%k x-out),k) U lI (x " Vt) - , a-b 2,2K 1 + b^d" sn ( ~ (k x *-«*>& .2 (a-b) (c-d) , , ., N k = (a-c) (b-d) (modulus) k = ^° / M.(a-c)(b-d) x 4K V 6h uu = vk. X K = K(k) (the elliptic integral of the first kind) 44a Fig. 3.7. Cubic nonlinearity (m- > 0): periodic waves I 45 The velocity is given by (3.20) and (3.21) v 2 = C [1 - ^(ab+bc+ca+ad+bd+cd)] o 6 = C q 2 [1 - ^-((a+b+c+d) 2 - (a 2 +b 2 +c 2 +d 2 ))] 2 = C Q 2 [1 - |- + Jjj(a 2 +b 2 +c 2 +d 2 )] . (3.24a) Or, eliminating d completely yields v 2 = C q 2 [1 + |(a+b+c) +^-{(a-b) 2 + (b-c) 2 + (c-a) 2 }] . (3.24b) If £ = 0, the velocity is always greater than C . For a nonzero e, the 2 condition e — 3 j-L < must hold. This condition guarantees the stability of the lattice; the interaction force is always attractive when the relative displacement is large enough. Solitary Waves In the limit b -> c , both solutions (3.22) and (3.23) approach to the solitary waves: , b-d .2, n N b(a-d) a(b-d) .2, d+a — - tanh(ox-pt) N / s l sech (ax-pt) v*- vt > - ,;:, -r - ■ l'-\ b-d ? R , — < 3 - 25 > 1+ — r tanh (ax-Bt) — sech (ax-Bt) a-b K a-b a-b K and a+d^ tanh 2 (ax-Bt) ^^X - *g£± sech 2 (ax-Bt) i , a ~ p , , 2 , _ , a-d a-b . 2 . _ . 1+ r~j" tanh (ax-Bt) r— r ~ 777 sech (ax-Bt) U II^> ■ ,.l, .2. ~ ■ l-i a-b ,2 < 3 - 26 > where a = J T^T (a - b)(b - d) and B = va . Hence the solitary wave solution in general is an infinite series of 2 sech (ax-Bt) even in this approximation. If we choose the special value 46 b = c = 0, the above solutions become simple. In this case we can obtain the zero background solutions. Since a+d = 0, and with t -\ 1-tanh (qx-Bt) , ,„ u ] .(x-vt) = -a- ^ u-* 1 - = -a-sech(2ax-2Bt) 1+tanh (ax-Bt) u (x-vt) = a-sech (2ax-2Bt) (3.28) and B = va 2 2 e v = C Z (l+ - a o J {?•'> (3.29) Periodic Waves around the Natural Length of Lattice Since P(w) is a quartic polynomial with the negative leading coefficient, bounded solutions are possible even if P(w) does not have four real zeroes. As a simple example, we will consider the periodic wave which is symmetric around the natural length of the lattice and also t=0. Such a wave can be realized by putting P(w) = (a 2 -w 2 )(w 2 +b 2 )p 2 . (Fig. 3.8) The equation is .dw. 2,2 2 . . 2 , 2 . (— ) = P (a -w )(w +b ) , dz solution of which is /~~2 2 2 (z) = acn(V(a +b ) pz,k), with k = 2 u 2 a +b (See Appendix B ) Or, in terms of x and t ,2K u(x-vt) = acn(— (k x-^)t),k) (3.30) (3.31) 47 PM — M> Fig. 3.8. Cubic nonlinearity (u. > 0) wave II (symmetric). periodic 48 where rr k x = 2K hlC 6h (a 2 +b 2 ) = f- C /£- ^ y 2K o J 6h k uu = vk K = K(k) . The velocity is given by the quadratic term in (3.30) 2 r 2 v -C o 2 2 o . L I o (a -b ) = 6>i (a - — +a ) = _, (2 - -) k k or v 2 = C 2 (1 + ^- (2 - ij)) . o 6 , I k Hence the wave is "supersonic" if k > — and "subsonic" if k < — . The >fl dispersion relation is 2 „2. 2,.,| 2, ,2K N , 4 a. =C k (1+ t 1 a ) - * ( — ) k o x 3 *i x (3.32) 2 2 2 4 As a -0, (3.32) reduces to uu = C k - *k , the dispersion relation for ' o x x r the linearized equation of (3.5) ii = C 2 u" + -hi"" . o Algebraic Solitary Waves ^l ' All the solitary wave solutions discussed so far (both quadratic and cubic case) decay exponentially as z ► +°°. But there is also an extra -2 solution for the cubic nonlinearity which behaves like z as Z * + 00 . Such a solution exists if P(w) has a triple zero (Fig. 3.9). Namely, (~) = p 2 (w-a) 3 (b-w) - ; P(w) (a < b) dz (3.33) has a solution w(z) = a + (3.34) l + (2bpz)' 49 Fig. 3.9. Cubic nouLinearity (u. > 0) : algebraic solitary wave, 50 To simplify the formulation, we assume that £=0. Since P(w) = - p (w -(3a+b)w +3a(a+b)w -a (a+3b)w+a b) , we obtain the following equations from the cubic and quadratic terms 2 o-(3a+b) = - -~- = 2 2 2 ^ - C o > 3ap (a+b) = - f— . From the first equation b a- - j • From the second equation 2 2 r 2 ^fo_ ,2 ,2. V " C o = ~ ( 3 b } ' The solution (3.34) may be rewritten in terms of b alone: / n b b w(z) = - - + ~ . l+(2bpz) If we define the amplitude of the wave as A = a-b = !>■ w(z) = ff-1 + 3 C 2 , o , 2 2 1-K-r A z 8*< The velocity-amplitude relation is v 2 = C 2 (1 +T ^A 2 ) 16 The algebraic solitary wave differs greatly from the ordinary 2 (sech type) solitary waves. The differences are summarized as follows: (1) decay is much slower (Lorenzian) ; (2) the amplitude is large and the condition to allow such wave is very strict (triple zero required for P(w)). 51 (3) the wave must cover both the positive and the negative w regions. In order to obtain the zero background wave, we must introduce the nonzero e, in which case the solution yields w(z) = - — 1 + — z z 3 ho. Hence the amplitude is determined uniquely. (4) The velocity is slow compared with the ordinary solitary wave of the same height. Since the algebraic solitary wave exists under very limited conditions, we will exclude this type of solution in the future consideration. 3.3.2 u. < (Soft Nonlinearity) We can put 2 P(w) = p (w-a) (w-b) (w-c) (w-d) where 2 M C o 2 a+b+c+d = -At (3.36a) IJH 22 2 v -C p (ab+bc+cd+da+ac+bd) = — . (3.36b) Periodic Wave There is only one region of w in which bounded solutions exist (Fig. 3.10). The solutions are: a /hz£.\ 2 / V(a-c) (b-d) . . c-d.(^)sn Q" '± pz, k) w ] .(z) = — ^ ^(a-cHb-d) pzk) 52 PM Fig. 3.10. Cubic noniinearity (m- < 0): periodic wave, 53 or b-a < a-c l , P-c N 2 /v/Ca-c) (b-d) , . b-a( )sn (^ ^ <— pz,k) w n (z) = /b-c. 2 , y(a-c)(b-d) 1 - (~ )sn <~ ^ *- Pz,k) As is discussed in Appendix C, these two solutions differ from each other only by phase. That is, we can match w and w by shifting one of the lem by 2K/ps/(a-c)(b-d) . Going back to the x,t coordinates, w and w are, respectively, c-d(^7)sn 2 (— (k x-uut),k) b-d rr v x U I (X_Vt) = : ,b-c, 2,2K„ ~ (3 ' 37a) and 1 (— )sn (— (k x x-u)t),k) b-a( )sn ( — (k x-0)t),k) u n (x - vt) - ; 'b-c, 2,2k" — < 3 - 37b > 1 - ( )sn (— (k x-o)t),k) a-c rr x ' ' ' where C TT _ _o_ / u,(a-c)(b-d) x " 4K V 6h 0) = vk x K = k(k) . The velocity is given by (3.36a) and (3.36b): 2 „ 2 r . e 2 LL 2 L 2 2 ,2^ v = c o [1 + iR " "ri (a +b +c +d )] : Q 2 [1 + j(a+b+c) - -^-{(a-b) 2 +(b-c) 2 +(c-a) 2 }]. (3.38) 2 e 2 We can see that the velocity is always less than C (1+ r-i — r) . In par- o 3 | p. I ticular the velocity is less than C if e = 0. We should note the periodic o r wave solution is always unstable for a large enough amplitude when u- is negative. Therefore the only amplitude that yields the right hand side of (3.38) positive is admissible. 54 Periodic Wave around the Natural Length Since we consider a symmetric wave around w=0, we put b=-c > 0, hence a=-d > (Fig. 3.11). We also assume that e=0 for the sake of simplicity. Therefore . 2 4ab k = 2 (a+b) C TT i x 4K J 6k and = C 2 [1 - f {(afb) 2 +2b Z -2ab }] = C Z [l - £{(1- \) (a Z +b Z )+2b Z } ] < C 2 , r , H- r, k 2 w „, 2 ^ „ : v O " t> " ' O DZ O From the above three equations, we have 2 - C 2 Cl - *(1- t-><«=) S5.k2.ttb2] o 6 v 2 V C tr u x 3 o . C 2 (1 - f b 2 ) - «A (I" ^T)k 2 , (3.39) O J ' ' Z X or 2 2 „ 2, 2 ^ C o ,4K 2 . k 2 4 „ , m ud = C k - — — - h(— ) (1- — )k . (3.40) ox 3 n I x Solitary Wave If either a approaches to b in (3.37a) or d approaches to c in (3.37b), we have the solitary wave solution on the positive (rarefactive) background or on the negative (compressive) background, respectively. 2 2 Since sn * tanh , these solutions are, respectively, c-de^Hanh (ax-pt) u(x-vt) = ^ 5 (3.41a) 1 - (^)tanh (0K-pt) and b-a(— )tanh 2 (ax- ( I I u(x-vt) = ^ ~ (3.41b) 1 - (— )tanh (ox-^t) a-c 55 Fig. 3.11. Cubic nonlinearity ((J, < 0) : periodic wave (symmetric). 56 where r C Q 2 (b-c)(b-d) a = -* ^ *- P = / — for (3.41a) (a-c) (b-c)(j,C 2^ — for (3.41b) P = VO! and o v = C (L- £(b-c) ) for both cases o 6 2 u 2 Note that there is no solitary wave solution on zero background: if a=b were 0, it would give a+b+c+d - c+d=0. But this is contradictory to d < c * 0. On the other hand, if c=d were 0, it would give a+b=0 which contradicts the condition a > b -> . Hence the above proposition is proved. This corresponds to the fact that — , — , where f(x) is the dx force function, takes the maximum Extra Solution There is also an extra solution for K ■ 0. But this solution is quite different from the algebraic solitary wave for ,- > 0: it requires the positive background on one side of the lattice and the negative back- ground on the other. Since we are not interested in this solution, we will give only a brief description. We assume that P(w) has two double roots a and -a (Fig. 3.12). The equation is 2 ,dw. 2 . . 2 . .2 (— ) = P (w+a) (a-w) dz 2 2 2 2 = P (a -w ) 57 Fig. 3.12. Cubic nonlinearity (|i < 0) wave solution. non-solitary 58 the solution which is w = a • tanh (apz) v 2 = C 2 (1 - ^ ) o 3« Table 3.1 summarizes all the solutions discussed in »3 . 2 and <3.3. 3.4 Piecewise Linear Approximation and Explanation of the Existence of Solitary Waves So far we have treated the polynomial type force functions exten- sively. We have made use of the special properties of the elliptic func- tions in deriving the solutions. Some of the common characteristics of the solutions are: (1) the periodic solution can be either supersonic or subsonic depending on the wave form; (2) the solitary wave solution is a special case of the periodic solution; (3) the solitary wave is always formed in such a way that the interaction force becomes stiffer as the displacement increases; (4) the propagation velocity of the solitary wave is always higher than the ambient sound velocity; (5) the wave form of the solitary wave may be approximated by the exponen- tial function in the region where the displacement is small. On the other hand we learned in Chapter 2 that the harmonic lattice has the exponential solutions as well as the sinusoidal solutions. This fact together with (3) and (5) suggests the lorce described by two linear func- i inns as a simple model: 59 x> H C c« CM c •H CO C o •l-l 4-1 a r-i o CO CO CO o a. >i V-l CO O •i-l u O o II ri c73 b o u = C 2 (l+A)u M + m"" u < b . o We assume that the solution u(x-vt), its first derivative and its second derivative are continuous at u=b . Putting w(z) = b-u(x-vt) and integrating the equations for w(z) twice with respect to z yield 2 2 hw" = (v -C )w + D.z + E. (w 0) o 1 1 H W " = [v -C (1+A)] W + D z + E (w > 0) (3.42a) (3.42b) where D 1 , E , D„ , and E ? are the integration constants. As we discussed in .^3.1, we can put D = D = 0. From the condition that w" is continuous at w=0, we have E. = E = E. Multiplying the both hand sides of (3.42a) and (3.42b) by w' and integrating with respect to z yield 2 ^w* 2 = (v 2 -C 2 )t + Ew + 2* o '2 w ■ M 2 r 2 9 1W rw 1 = [v -C "(1+4)]- + Ew + F w ■ I o I I (3.43a) (3.43b) where F 1 and F_ are the integration constants. From the condition that w' is continuous at w=0, F = F_ = F. Dividing the both hand sides of the above equations by ~ yields 2 -r 2 2 V o 2 9 E 2F w' = w + — w + — = P, (w) w < (region I) 1 (3.44a) v 2 -C o 2 d^) 2 2E 2p . n ^ . Tn w' = w + — w + — - P ? (w) w > (region II) (3.44b) C o 2 A Obviously, P (w) = P,(w) - -^ — w 2 i. I I • < 61 1 X) 1 1 1 JLj I 1 1 / ^r ! / -N / I 1 > 1 ' / 1 if /l / 1 / 1 / 1 / 1 / 1 Fig. 3.13. A typical piecewise linear force function (one break point). 62 The solutions of the differential equation of the forms given in (3.44) are known to be either the sinusoidal functions or the hvperbolic functions, the exponential function being included in the latter as a special case. More explicitly, if the quadratic polynomial P (w) or P„(w) has the positive leading coefficient, the solution is +cosh, sinh or exp depending on whether the discriminant of the polynomial is positive, nega- tive or zero, respectively. On the other hand, if the leading coefficient is negative, the solution is the sinusoidal function, provided that the dis- criminant is positive; otherwise there is no real solution. In solving (3.44) we can select one of the above mentioned functions in each region, depending on the leading coefficient and the discriminant of P (w) or P-(w), then connect them in such a way that the derivatives of the selected functions match up to the second order at w=0 . There are two different cases of our interest. 2 2 3.4.1 v C (Subsonic Waves) o Both P (w) and P (w) have the negative leading coefficients. Obviously, we have the periodic solution in (w ,w J where w is the smaller zero of P„. (w) and w is the larger zero of P. (w) . (Fig. 3.14.) The solution may be obtained bv connecting two sinusoidal waves at w=0 . The resulting wave is not symmetric with respect to the z-axis due to the difference between P (w) and P (w) . In order to maintain the wave form as the composite wave propagates, the phase velocity of the wave in each region must be the same. This is best explained graphically by using the dispersion curves (Fig. 3.15). Since the force function is steeper in region II, the disper- sion curve for this region is always located above that for region I. For a given phase velocity v , we draw a straight line in the x-k plane, which P x 63 Fig. 3.14. Subsonic wave in piecewise linear lattice. u)=V- p & x Fig. 3.15. Dispersion relation for subsonic waves 64 crosses two dispersion curves at k and k . This is possible for a suitably chosen v . Then we put the solution in region I as w (z) = A cos (k z) + B (3.45a) and the solution in region II as w i;L (z) = +A 2 cos(k 2 z+0) + B 2 . (3.45b) We have the following six equations to determine A , B , P, A 9 , B , z from (where w (z*) = w (z*) = 0): A + B = w > A cos(k z ) + B = B 2 - A 2 = w 2 A 2 cos(k 2 z'+^) + B 2 = A k sin(k z ) = A k sin(k z +d) 2 * 2 A k cos(k z ) = A k cos (k z +p) (3.46a) (3.46b) (3.46c) (3.46d) (3.46e) (3.46f) From these equations we can eliminate all but A and A_ , result- ing in the following two equations, each representing a hyperbola: l - A : 2 2 2 1 w 2 (k 2 -k L ) 4 2 2 2 K ■;■ w 2 (k -k )n ± : a i ± 1 — 2 2 2 2 • ? 2 k, V (L -k, ) ^ k/ 1 2 v 2 1 = 1 (3.47a) 4 2 2 2 k w i^ k 2 ~ k l ^ 2 2~ 2 . 2 ! A l + k 2 W] _ (k 2 -k i k J- W l ^ 2 1 ' 1. (3.47b) We can solve (3.47) for A. and A„ graphically (Fig. 3.16). Since there is only one intersection of the two graphs, A. and A_ are determined uniquely for the given set of k , k. , w and w . We can obtain all other unknowns in a straightforward manner from (3.46). By 65 Fig. 3.16. Graphical solution of (3.47) 66 changing v , it is possible to obtain the relations among such quantities as the amplitude (w +w ) , the wave number and frequency, but we will not go into the detailed discussion, since we are interested in showing the existence of the solutions similar to those in §3.2 and the feasibility of the piecewise linear model rather than the numerical calculation. 3.4.2 C < v ,n + it)) (3.51) whi i - i = sinha. (3 . 52) 71 Note that the relation (3.52) is identical to the one given for the unbounded solution of the harmonic lattices (§2.1). This is not a coincidence. Toda's solution (3.51) approaches to the exponential func- tion as an + (3t -* co . That is, the outskirts of the solitary wave behave like the wave in the harmonic lattices, for which B = sinhct holds exactly. Since Toda's solution is a traveling wave, the same relation between a and 6 should hold all over the wave. Hence, (3.52) seems to hold quite generally for the solitary waves in the anharmonic lattice. If we expand (3.51) and (3.52) in the Taylor series, the first order terms are u = -R 2 sech 2 (an + St) (R 2 « 1) n and 3 2 2 S 2 = (crt-^) ~ a 2 (l+^") . From the latter equation the velocity yields 2 2 v 2 i. i- 3 , a 2 which matches (3.13) if we put c =1 and £ = 1. As for the periodic wave, the exact solution is u n = -log[l + (2S2) {dn 2 ( 2 ^(k x niu)t) - |3] (3.53) where UJ and k are related in the following way: x s n ( — k ) 2KC0 5L ^ tt ^x TT (3.54) 7l+ (|> l)sn 2 (^k x ) The relation between the phase velocity and the wave number is obtained from (3.54) and is shown in Fig. 3.21. The curves are very close ^ 72 1.1 k=.W9 j i i i 7C to 10 30 40 50 60 7o ?o Jo Fig. 3.21. Phase velocity vs wave number relaf.ion for Toda ' s lattice . 73 to that for the harmonic lattice if k is small. As k increases the velocity starts to take values greater than C for the lower wave numbers. Clearly the periodic wave for Toda ' s exponential lattice can be either ? 2 2K(J0 supersonic or subsonic. The Taylor expansion is valid when k ( ) « 1 2 E 2 since the amplitude of (dn - — ) is k . Therefore either k or ()) should be K small compared with 1. Then the first order approximation is u n - -(^) [dn 2 {^(0)t+k x n)} - |] . (3.55) The dispersion relation yields 2 ( 2K 2 _ 1*£ ( 2K < 2K0J l _. • tt V 3 ' TT V , 2K , ,2 .likl.E 1W _V W ( ~ k x } " ( ~l~ + i- 1)( — > 4 = (^ k x ) 2 + }(2 - k 2 - ^fX—^) (3-56) where we have made use of the following approximation: 2 sn(x,k) ™ x - — — x (k.x « 1) . o We can easily show that (3.56) matches (3.11) if we put C = 1, h = — . It is expected that the Taylor expansion of (3.53) up to the second order matches the solution for the cubic nonlinearity (§3.3), pro- vided that the latter solution (which is in the form of the rational func- tion) is also expanded in the Taylor series and the terms with the order of sech and higher are truncated. Detailed calculation will not be made here . Toda also discovered the two solitary wave solution for the exponential lattice. He could show that there is a phase shift after the 74 two such waves collide each other. The phase relation is such that if two waves make a head-on collision the phase shift is positive for both waves, whereas if one wave takes over the other, the phase shift for the faster wave is positive and the phase shift for the slower one is negative. The important thing is that these waves preserve their identity even after the collision, unlike the conventional shock waves in the dissipative medium. Recently Wadati and Toda showed that both the original equation of motion (3.46) and its continuum approximation (3.5) with ^=0 give exactly (32) the same phase shift for the two wave collision problem. Therefore, we can expect that the continuum approximation gives quite good quantative results as well as qualitative results. Zabusky discussed in his recent report that a travelling sine wave in the anharmonic lattice breaks into many small solitary waves of different amplitudes and velocities ("L Soliton"). Since Zabusky chose the periodic boundary condition, we can expect that these small waves are actually the superposition of the periodic waves of different amplitudes ind wave numbers, each of which are similar to Toda ' s periodic solutions (33) 75 CHAPTER 4. NUMERICAL METHODS AND RELATED TOPICS As we introduced in §1.2, our model consists of two essential parts: the crystal lattice described by a set of ordinary differential equations and the heat reservoirs described by random sampling of the velocities. We need only the former part to investigate the transmission properties of the lattice (deterministic problem). As for the lattice we have chosen the free boundary conditions on both ends, thus obtaining the following equations A, dt d 2 y. 2 = f (y 2 -y].) (left end) dt d 2 y, dt 7= f(y ,,-y )-f(y -y -, ) (n=2~N-i) 2 ^n+1 n n n-1 (4.1) = -^vW (right end) where y is the displacement of the n 1 ^ atom from its equilibrium position and f is given as 2 3 f (x) = x + ax + px (4.2) We can also write (4.1) in a vector form: d 2 —7 y = F(y) dt where r> i s y N ! -J f(y 2 -y L ) y 9 | » F =: f(y,-y 9 )-f(y 9 -y,) | 3 J 2' w 2 J V •^vW 76 We also use the vector notation in this chapter when we mention the system of ordinary differential equations in general, in which case F may depend on y,t as well as on y. If F depends on y only we call such system of equations "special equations" according to Henrici. Through- out this chapter O.D.E. stands for the ordinary differential equations. In §4.1 and 24.2 we discuss the numerical integration schemes. 34.3 we discuss the heat reservoirs more in detail. In ?4.4 we discuss the computer system which we used and also the semi-interactive features of our programs. In §4.5 we show the typical structure of our programs by giving an example. (Also in Appendix D.) 4.1 Numerical Integration of the System of O.D.E. Generally speaking, there are two different types of numerical integration schemes for the system of equations such as given in (4.1): (1) One step methods (Runge-Kutta type) (2) Multistep methods (Predictor-Corrector Methods). As for the special equations, we can solve such equations directly without reducing them into a system of first order O.D.E. with 2N variables. There- fore we have four different types. First, we compare the one step methods with the multistep methods. (1) The one step methods are self-starting, that is, we only have to specify the initial conditions. On the other hand, multistep methods require the information of several time steps before since in these methods we solve the finite difference equations. Especially in our problems, as will be clarified in §4.3, we have to restart the integration everytime the end atoms interact with the heat reservoirs. Therefore the mu lii step methods are very difficult to use. 77 (2) One step methods usually require more number of function evaluations . per time step than the multistep methods do. A typical fourth order one step method requires four function evaluations per time step, whereas a comparable multistep method requires only two function evaluations per time step. In most cases this is the major obstacle against one step methods, since the function evaluation is the most time-consuming part. But in our present problems, as Waters extensively (35) discussed, the function which we deal with is the polynomial (cubic at most) and is very easily evaluated. Hence the one step methods are quite promising. Gear also mentioned this. As for Toda ' s poten- tial, however, this is not the case; if we need to do a long time computation with Toda ' s potential, the mixture of the multistep method (for the inner atoms) and the one step method (for those atoms which are near the boundary) will yield faster computation. (3) One step methods can be programmed very simple. (4) Error bounds" for the multistep methods are more explicitly given than those for the one step methods. For the one step methods we have to evaluate the error by, for example, the step reduction technique. From the above consideration we will solely discuss the one step methods for the integration of the equations (4.1). * Remark given by C. W. Gear, There are two different types of errors in the numerical integration of O.D.E. One is the truncation error which stems from the approximation formula for integration. The other is the round-off error which stems from the finiteness of the length of the registers in the computer. Here, we are solely concerned with the former. 78 4.2 One Step Methods for the Special Equations Usual technique to solve the second order O.D.E. is to reduce them to a system of first order O.D.E. by introducing extra variables z =y . n ^n There are quite a few integration schemes for the system of first order O.D.E. Most popular ones are Runge-Kutta ' s fourth order method (to be referred to as RKC4 in the following discussion) and its variations by Ralston and by Gill. These methods are listed in Table 4.1. Ralston's method was derived as a result of minimizing the truncation error of order h , where h is the step size of integration. This method gives the truncation error which is about 1/2 of that for RKC4. At the same time, however, the simpleness of the integration scheme is lost and more compu- tation time results. Gill's method was derived from the viewpoint of /TON reducing the size of temporary storage rather than an error consideration. Blum later showed that the same reduction of temporary storage is possible (34) for RKC4 too. Therefore Gill's method is quite irrelevant nowadays. (Note that the simpleness of integration scheme is also lost in Gill's method.) Hence we can naturally select RKC4 as the most optimal one among * the fourth order one step methods. For a higher precision computation Waters and Pasta et al_. recommend Butcher's fifth order method. It (18 ) is reported that the truncation error could be reduced as small as the round off error with a choice of h = .01. This is a "super accurate" example. Usually the truncation error is several orders larger than the roundoff error. We do not require such a high accuracy in our computation According to (39), p. 363, the Runge Kutta subroutines in most computer libraries still employ Gill's method. Probably this is just because the original (and the best!) method was derived in 1901, whereas Gill's method was derived in 1951, the latter giving an impression that it is quite i mproved . 79 Classical Runge-Kutta (1901) Table 4.1 (34) W - *n + t ( «l + 2t 2 + 2lt 3 + *4> K. = F(y + hK„) 4 n 3 (37) Runge-Kutta-Ralston (1965) y , = y + h(.17476028K 1 - .55148053 K_ •'n+l n 1 2 + 1.20553547K- + .17118478K. ) 3 4 K L = F(y n ) K_ = F(y + .4hK. ) z n 1 K = F(y + ,29697760hK„ + .15875966hK ) 3 J n 2 2 K. = F(y + .218100381^. - 3 .05096470hK o + 3 .83286432hic o ) 4 n 1 2 3 Runge-Kutta-Gill (1951) ^ 34 ^ ?n+l = V |(^+(2W2)K 2 + (2-V2)K 3 + K 4 ) *2 = ^n + 2 K V K 3 = F(y n + IC-l-K/2)^ + |(2- V / 2)K 2 ) Table 4.1 (continued) Nystrom's method for the special equation (1925) , 2 (34) Vfl = y n + h 'K + I^ (23K 1 + 75K 2 ~ 27K 3 + "V 80 ?Ul = K + I9T (23K 1 + 125K 2 " 8LK 3 + 125K 4 } h = f ^n } K 2 = F(y n+ fh?; + ^ h 2 K L ) h = ? V 3 " Z n + 3*1 *3 " ?(7 n + ¥ 3 ' V 4 " Z n + 5 (K 1 + V -► "*,-> 4h-* „ K 4 " ' " f (n=2 ~ "- 15 at ^N „, . % dt 2 = " f(y N" y N-l } " ^N~dT+ B N (0 ( ri § ht end > where and B 1 (t)B 1 (f ) = 2T 1 6(t-t') V t)B N (t,) = 2 V (t - t,} 86 Nagazawa chose |i and |i that yield the optimal contact between the lattice and the heat reservoirs for the highest mode of the lattice. ' 4.4 Computer Facilities and Programs The computer which we used for our numerical experiments is Burroughs B6700 system located at the Center for Advanced Computation at the University of Illinois. The word length is 48 bits with 39 bit mantissa. Its characteristics are implied by the following figures: addi- tion time .20o.s, multiplication time 2.00jj,s and division time 10.8|j.s. B6700 has a special architecture designed suitably for processing ALGOL language. It is not primarily designed for high speed scientific computa- tions, but its flexible operating system (MCP) allows users to manipulate remote terminals, on-line monitoring, etc., without requiring much exper- iences in system programming on users' part. Fig. 4.2 shows the B6700 system as far as our computation is concerned. Our programs are stored on a 9 track magnetic tape and are easily read on the disc memory. Editing of programs can be done from the remote terminals, from which an immediate access to the disc memory is possible. The remote terminal which we used is the Hazeltine 2000. It is a character display terminal (27 rows x 68 columns), having the cursor addressing feature. This latter feature allows users to write programs which plot curves (although very crude since its resolution is 25x60 points) while the main program is being executed. This is very convenient when we monitor what is going on in the program. By making use of the terminal just described, we wrote our pro- grams in a semi-interactive way. Some of the features of our programs are: 87 REMOTE TERMINAL LINE PRINTER Fig. 4.2. Computer facilities 88 (1) A command to execute a program is given from a terminal. This terminal will be tied up to the program until the program is terminated. (2) Input parameters are fed from the terminal. There are also some options for printing only the desired variables such as the coor- dinates, velocities, temperature and flux. (3) Intermediate results are displayed on the terminal at an option. Time interval can also be set at which we require displays. (4) Programs can be terminated at any time by sending a command from the terminal if the user does not need it any more. (5) An interactive use of the terminal is possible only when the program reaches to the points at which it requires some information from the user. Therefore our programs are not interactive in the strict sense. But even semi-interactive programs of this level still enable us to use the computer system quite efficiently for the heuristic research. Instead of creating a huge, complicated program that handles everything, we have made separate programs for considerably different tasks. Table 4.3 gives a list of programs which we have created. 4.5 Description of a Typical Program As a typical example of our programs we describe SOL/XXIV that was written for computing the heat flux and the temperature distribution. A simple flow chart is given in Fig. 4.3. In this program the heat flux and the temperature are computed and printed (also on the terminal) at every 50 seconds. The temperature gradient is computed by the linear least- square interpolation. The entire program is in Appendix D. 89 Table 4.3 List of Programs SOL/XXIV ZABUSKY IMPACT IMPACT2 DISORDER DISORDER2 DALT BRKLIN RESPONSE vise RKC4 NYS5 SPHERE FPU FPU2 HAZELTINE FPU3 heat conduction (JPW) periodic boundary condition (Zabusky) creation of solitary wave (JPW) creation of tail (JPW) one impurity at n=25 (harmonic) impurities at n=25~50 (harmonic) impurities alternatively (harmonic) piecewise linear force cooling process of the lattice (JPW) viscous ends (Jackson's problem) Runge Kutta Classical test Nystrom test hard sphere •+- harmonic recurrence phenomena (cubic force) recurrence phenomena (quadratic force) test program for the terminal recurrence phenomena (free ends) INITIAL CONDITIONS, PARAMETERS 89a INTERACTION IN LEFT RESERVOIR INTERACTION IN RIGHT RESERVOIR NYS5 INTEGRATOR (ONE STEP) TEMPERATURE TERMINAL TOTAL ENERGY , FLUX TEMPERATURE LEAST SQUARE YES PRINT I * Fig. 4.3 Flow chart of a typical program (SOL/XXIV) END 90 CHAPTER 5. DYNAMIC PROPERTIES OF ANHARMONIC LATTICES This chapter deals with the wave motions in the anharmonic lattices. Specifically JPW lattice (which is characterized by the force function 2 3 f(x) = x-lOx +66. 7x ) is considered in full, but the piecewise-linear force model and the hard sphere plus linear force models are also treated. Since no analytical solutions are known, the numerical methods are employed and the results are compared with the analytical solutions in the con- tinuum approximation whenever possible. The solutions in Chapter 3 were considered in the special cases; the solitary waves and the periodic waves in the infinitely long lattices, or equivalently , under the periodic boundary condition. Most of the numerical (33) experiments have been performed along this line. In the heat conduc- tion problem there are new factors. Firstly, the lattice is necessarily bounded by the heat reservoirs. Therefore the effect of the boundaries must be taken into account. Secondly, the initial condition is different from that for a pure solitary wave solution; we assume the impulse type interaction between the lattice and the heat reservoirs. In fact, an impulse produces an oscillatory wave packet together with a solitary wave. The main purpose of this chapter is to clarify the basic energy transport properties of the anharmonic lattices from the numerical experiments. Most of the results in this chapter have not been reported. From §5.1 to §5.4 we consider solely the JPW lattice. In §5.1 we consider the creation of the solitary wave by the impulse. In §5.2 we will compare various boundary conditions and the reflection of the soli- tary waves therein. We will show that the fixed boundary is a perfect reflector but that the free boundary destroys the solitary waves. In §5.3 91 we will consider the two wave collision process. We discovered that the solitary wave slightly loses its energy in the collision process. We will also mention the wave propagation through the disturbances. In §5.4 we will discuss the behavior of the time averaged kinetic energy of the atoms when the waves propagate through the lattice. This quantity defines the local temperature in the thermal situation. In §5.5 and in §5.6 we will consider the piecewise linear force function and the hard sphere type force function, respectively. We will numerically show that the solitary wave can exist in these models. The results obtained in this chapter will be used in Chapter 6 in understanding the heat conduction phenomena. We should note that we will employ the coordinate space which is defined by the atom number n and the time t rather than the spatial coor- dinate x and the time t. The former definition corresponds to the Lagrangian coordinate system and the latter to the Eulerian coordinate system in hvdrodynami .s . The advantage of the former is that the atom number n is h h fixed all the time whereas the n atom itself moves in the one dimensional space. This difference is very important when we consider the collision process. The approach we took in deriving the partial differential equation in Chapter 3 relates the (n,t) space in the lattice with the (x,t) space in the continuum approximation. 5.1 Basic Dynamic Properties of JPW Lattice 5.1.1 Creation of the Solitary Wave by the Impulse We discovered that an impulse applied to the end atom of the lattice produces a solitary wave and a tail. By the latter we mean a small wave packet which consists of high frequency waves. A typical profile of such 92 wave is given in Fig. 5.1a. Separation of the two parts was observed at an early time for the initial velocity greater than .08. If the initial velocity is small (<. 05) separation was not observed within the limit of our lattice size (N=100) . Such an example is also given in Fig. 5.1b. In this case the entire wave propagates as a (dispersive) wave packet. Presumably separation will take place after the wave propagates for a long distance. Zabusky also found a similar type of the wave as a solution of ft-1) the KdV equation. Unfortunately, there is no analytical solution known for such composite (solitary wave plus tail) waves. First we look at the solitary wave part. We performed numerical experi- ments with different initial velocities. Results are in Table 5.1. Fig. 5.2 shows the relation between the amplitude of the wave and the propagation velocity v. The solid curve is the theoretical prediction from the con- tinuum approximation (3.29): 2 2 e u, 2 v=C(l+-A+ fk) o Jo 2 with C = 1, E = 2x10, M-=3X66.7 and a=-A. Experimental data are always located below the theoretical curve. This tendency becomes manifest as the amplitude increases. Apparently the truncated terms in the continuum approximation are responsible for this effect, since the wave form becomes sharper as the amplitude increases. Fig. 5.3 shows the relation between the initial velocity V. and the propagation velocity v. We can derive an empirical formula v = 1+1. 85V. in the range [.1, .4"]. 5.1.2 Energy Concentration Since the propagation velocity of the solitary wave is much faster than that of the tail, we can calculate the energy contained in each part separately. The energy of the n atom is defined as Table 5.1. Dynamic Properties of Solitary Waves 93 V. 1 Amp . V E . sol tail E onl CD sol h Error .08 .0364 1.13 .00292 .00028 91.1 .02 .5 10" 8 .1 .0465 1.19 .00467 .00033 93.4 .05 .9 10" 5 .2 .0900 1.39 .01960 .00040 98.0 .025 <. 5 ID" 8 .3 .1252 1.57 .04455 .00045 99.0 .05 .2 lO" 6 .4 .1550 1.75 .07952 .00048 99.4 .05 .6 lO" 6 .5 .1819 1.89 .12450 .00050 99.6 .02 5 .8 lO" 7 .6 .2066 2.02 .17940 .00060 99.67 .05 .3 l0 - 5 .8 .2482 2.27 .31931 .00069 99.8 .05 .15 lO" 4 1.0 .2840 2.51 .49920 .00080 99.99 .025 .1 .1 .6 lO" 5 1. 2 .3199 2 . 67 1 . 5 .3640 2. 97 V.: initial velocity l Amp: amplitude of the wave V: propagation velocity of the wave E , : total energy in the solitary wave sol E . • total energy in the tail tail E, .(%): E /(.5 V 2 ) sol sol 1 h: integration step Error: numerical error in total energy means that such quantity was not computed. 94 .OX -.02 - -.0*- -.06 - (a) V. = .2, t = 15 a. Tl Fig. 5.1. Typical profiles of waves created by impulse 95 Amplitude z.o 1.0 • 6 .1 .1 •Of • 02h •01 i i — i i i i i t 1 — i — i — r—r Solid Curve V*=1+6.67A + 33.3SA* 1.0 -t — i — i i i 1 1 1 ■ i ■ . ■ ■ -i i i i ■ » ■ ■ .i .♦ 1.0 2.0 10 V-i Fig. 5.2. Velocity vs amplitude relation. 96 ix V Fig. 5.3. Initial velocity vs propagation velocity, 97 E n = ly/ + *[v(y n+L -y n ) + V^-y^)] (5.1) where V(y -y ) is the potential energy stored in the spring connecting the n and n+1 st atoms. Results are in Fig. 5.4. Numerical data for the JPW lattice are also given in the fourth column of Table 5.1. The stronger the nonlinearity is, or the larger the amplitude is, the more fraction of the energy is concentrated in the solitary wave, although the energy in the tail itself also increases. Fig. 5.5 is a typical example of the energy distribution in the JPW lattice. Fig. 5.6 shows the relation between the energy and the propagation velocity. From this graph we can obtain an 2 empirical formula which is similar to E = mv /2 : E = M (v-l) 2 /2 (5.2) s s 2 where M = .24 (L+.35(v-l) ) (effective mass!). It is obvious from the above results that the anharmonic lattice can transmit energy more efficiently than the harmonic lattice because of the nondispersiveness and the high propagation velocity of the waves therein. 5.1.3 Tail Part As we mentioned in 5.1.1 an oscillatory wave is always created with a solitary wave. Its amplitude is very small compared with that of the soli- tary wave which it accompanies with. We can assume that the tail is the harmonic dispersive wave. This was verified by the following numerical experiment : (1) Create a solitary wave and a tail with an impulse and let them propagate until they are separated from each other; (2) erase the solitary wave; (3) reverse the direction of the time and perform the numerical integration backward until t=0; 98 &»erj/ i%) Ito % 90 70 50 T JPW to -<* Fig. 5.4. Energy concentration in solitary wave (7<>) . 99 10-' - I0* v - 10 1 r 1 1 1 1 1 — i 1 1 1 T- T A / / / / \ \ \ \ f 3 — i \ \ ■■ / \ - / / i - / \ / / \ \ m / / / / i \ i \ r* / \ — / / i • 1 l i . * i - 1 i 1 1 i - 1 l 1 1 l 1 I - -j 1 l 1 ™ " 1 i ~ ■ 1 1 \ \ - i \ ■ I 1 \ - 1 i \ \ / \ / / l \ \ - r*, 1 i 1 / • i i i i i i I • i i i i 1 1 - ... 10' - 5^ 5J 56 $7 5J SI lo n Fig. 5.5. Typical energy distribution in solitary wave, V. = .2, t = 37, h = .05, AE =.5xl0 -7 . 100 Eherjy »•» Society U/ave .01 ool T 1 ■ I I / I I I I I I I -i — i— i — r-r t / / / j 1 — i i i i i | 5oli' time 10 20 Fig. 5.9. Reflection of solitary wave at fixed boundary, 105 5.2.2 Free Boundary (Fig. 5.10) The solitary waves decay at the free end when the amplitudes are small, whereas those with large amplitudes can be reflected. This difference may be explained in the following way. When the solitary wave reaches the end of the lattice, atoms near the end tend to travel in the same direction, resulting in the expansion of the lattice there. This newly created expan- sion wave can propagate through the lattice back again as a solitary wave only when the force (in the positive region of the displacement) becomes stiff as the amplitude increases. In the case of the JPW lattice, however, the force at first becomes limp due to the negative quadratic nonlinearity , then it becomes stiff again due to the cubic nonlinearity which dominates the quadratic one. Hence there is a threshold of the amplitude below which the solitary waves become unstable. To make this point clear we also per- formed a numerical experiment with the negative impulses applied to the JPW lattice. Fig. 5.11 shows the responses of the lattice to various inputs. The critical velocity turned out to be V. = .3. 5.2.3 Heat Reservoir Type Boundary In this boundary condition, the coordinate of the end atom is checked at every m integration steps. If the end atom is found to be inside the heat reservoir a new velocity is assigned. We performed the numerical experiments with various m. Since we are concerned with the boundary effect in this chapter rather than the thermal properties, we chose the predeter- mined values for V. and V , where V is the velocity with which the end l r r atom is returned. Specifically, V.=l and V =.1. Results are in Fig. 5.12 a-e. The amplitude of the reflected wave decreases as m increases. If m=l the end atom is effectively frozen at the boundary, therefore the boundary condition is very similar to the fixed boundary condition. As m increases, V L = 1.0 106 -.ol v = .a t»6* -ol - 10 15 20 IS 4 1 1 71 -.01 - Un. 03 - ■ol .01 f** ^r 10 U -n Fig. 5.10. Reflection of solitary waves at free boundary, 107 •OS • VL--.1 VL— * V : =-.3 n v t =M Fig. 5.11. Responses to negative initial impulses (t=20) 108 w (b) (C) (d) Yl (e) m 1 U n /0 eafer" reflection (#) 65.% 57.2 26.6 7.8 6.0 u, Fig. 5.12. Heat reservoir type boundary- Broken curves: incoming wave Solid curves: reflected wave v, = 1, \' - 1 - 109 the end atom can receive more energy from its left neighbor before it inter- acts with the heat reservoir, hence more energy is exchanged. We also examined the boundary condition (called "soliton eater") which supposedly yield the most efficient energy exchange for the solitary waves; the end atom does not interact with the reservoir until its velocity reaches the maximum value. The result is shown in Fig. 5.12e. The effect of the boundary condition on the heat flux will be discussed in Chapter 6. 5.3 Collision Process of Solitary Waves 5.3.1 Phase Shift in Two Wave Collision Process It was theoretically shown that the phase shift occurs in the two wave collision process in Toda's lattice (23.5). In the JPW lattice the phase shift was observed experimentally; both waves were advanced after a head-on collision, whereas the faster wave was advanced and the slower wave was retarded after a take-over collision (Fig. 5.13). The amount of the phase shift was larger for the small amplitude wave in both types of collision. This seems to indicate that the motion of "the center of the mass" is not affected by either type of collision, the theoretical result (42) obtained by Wadati and Toda for the KdV equation. We should note that the above mentioned phase shift is considered in (n,t) space. Since the lattice itself is contracted toward the center, the phase shift is less obvious in (x,t) space. The phase retardation for the harmonic lattice, (21) reported by JPW, is due to this difference ' ; there is no phase shift in (n,t) space. The phase shift observed when the solitary wave hits the fixed boundary (5.2.1) can be explained in terms of the two wave collision process. If two solitary waves of equal amplitude travel in the opposite direction in such 110 (a) Head-on collision. time. 70 ZO fO (b) Take-over collision Fig. 5.13. Collision process of two solitary waves Ill a way that a collision takes place exactly at the location of, say, m atom, the m atom never moves. Hence the fixed boundary condition is strictly satisfied (Fig. 5.14a). We can further conclude from the numeri- cal results in 5.2.1 that the solitary waves lose no appreciable energy in a collision process which is equivalent to the fixed boundary condition. 5.3.2 Energy Loss We can extend the above discussion to consider a different initial phase relation of two waves. Suppose that a collision takes place in the middle of two atoms. Such type of collision is equivalent to the reflection of the solitary wave at the fixed boundary, except that the interaction force at the boundary is twice as strong as that at all other places (Fig. 5.14b). Since the lattice becomes inhomogeneous , we can expect that the solitary wave may be deformed. We performed a numerical experiment on the two wave collision process along this line, and discovered that some of the energy is actually left at the site of a collision (Fig. 5.15). There are 50 atoms in the lattice, two end atoms having the initial velocities, .5 and -.5. A collision takes place between n=25 and n=26. The amount of the energy taken from the solitary waves is very small (.37o of the initial energy), but this is far above the numerical error. This numerical experi- ment was performed twice with different integration steps (i.e. .02 and .01); both gave the same results within the numerical errors. Hence the energy loss is not due to the numerical integration scheme. This newly created wave consists of high frequency components and is assumed to spread out in the entire lattice at t = °°. Obviously the continuum approximation cannot treat such waves properly. This fact may seem contradictory to the theoretical results because in (32) Toda ' s lattice the only effect of a collision is a phase shift. In 112 m-i m m*i m-i Fig. 5.14a. Fixed boundary Fig. 5.14b. Another type of collision process 113 IS -.1 X - 20 j is t-H t = 16 t=16 U n E n : h = .02 h = .01 h = .02 h = .01 ' la -.01736568 .01736566 . .00009400 .00009400 ' 19 -.15784123 .15784116 .02532067 .02532063 ; 20 -.10197037 .10197044 .09260386 .09260386 ! 21 -.00731625 .00731626 j .00607652 .00607653 i : 22 -.00756556 .00756556 .00004459 .00004459 I I 23 -.00259986 .00259987 .00005328 .00005328 i 24 -.01865858 .01865858 .00011343 .00011343 25 -.01257805 .01257805 .00015234 .00015234 1 ^E . — L .000000056 .000000002 ', Fig. 5.15. Energy loss in a collision process 114 fact, the theory states that two solitary waves are shown to exist asymp- totically (t=+°°) . But there are some other terms in the expanded form of Toda's solution, which vanish asymptotically. These terms are neglected in the theoretical argument. Note that these terms arise in the process of approximating the solution, not the equation. In our numerical experiments two solitary waves are created at both ends of the lattice, between which there is no disturbance initially. Therefore the initial conditions are different in these two cases. As we mentioned before the wave created in a collision process will disappear as t -*• °°, if we assume the lattice of infinite length; asymptotically there are only two solitary waves. The solitary waves in all. collision experiments propagated after the collision without any decay. No decrease of the propagation velocity due to the energy loss was noticeable. We will show in §5.6 that the energy loss in a collision is more obvious for the hard sphere type force function In the thermal situation collisions can take place at any location of the lattice with equal probability. Also the energy loss, although exists, is negligible in the small sized lattice. Therefore we can summarize that the only important effect of the collision in our model lattice (JPW lat- tice) is that the propagation velocity of the solitary waves is effectively increased due to the phase shift. 5.3.3 Effect of the Disturbances There is a question as to the propagation of the solitary wave when the disturbances (that is, the incoherent motion of the lattice) exist. Since the solitary wave is in the sea of disturbances under such circum- stances, it is very difficult to isolate the energy in the solitary wave part and to measure its change. The only possible way is to observe the propagation velocity. As we discussed in §5.1, the solitary wave is 115 uniquely characterized by its velocity. Therefore any decrease of the velocity means that the energy is being lost. Fig. 5.16 shows an example. this figure was obtained by JPW, but it was not publicized in their * (21) report . As we can see in Fig. 5.16, the velocity change is mainly due to collisions; the original velocity is regained after each collision. Hence we conclude that the interaction between the solitary waves and the disturbances is negligible for the lattices of order 100 atoms or so. This does not mean, however, that no such interaction is possible. For example, if the lattice size is quite large the solitary wave may lose significant amount of its energy due to the interaction with the distur- bances as well as due to many collisions with other solitary waves. As for the solitary waves of very small amplitude, we do not have any way to distinguish them from the disturbances. 5.4 Time Averaged Kinetic Energy We define this quantity as follows: M K (t) - 7 > y 2 (knt) (5.3) n t ^_ , n k=l where M is the number of time steps and t = MAt. This quantity gives the local temperature when the lattice interacts with the heat reservoirs and M ■* ". We calculated K (t) for the impulse type initial condition. We n considered both the harmonic and the anharmonic lattices. The author is indebted to E. A. Jackson for this figure. ** r 1 Saito e_t_ al_. reported that a significant decrease of the velocity was observed for the solitary wave in Toda ' s lattice . \^^> 116 p., C ■r-l W E O 4-1 O o 4-i m •r-l Z 1—1 00 •r-l fa 117 5.4.1 Harmonic Lattice (Dispersive Wave Packet) As we discussed in Chapter 2 the created wave packet is very dis- persive. That is, each component of the normal modes propagates at its own group velocity. Consequently K (t) is a decreasing function of n for a fixed time t (Fig. 5.17a). 5.4.2 Anharmonic Lattice (Solitary Wave Plus Tail) The solitary waye gives a uniform distribution as far as it has propagated. The tail, on the other hand, gives a decreasing contribution to K (t) due to its dispersiveness . The resulting distribution is shown n ° in Fig. 5.17b (V.=.l) and in Fig. 5.17c (V.=.05). The latter is a non- separating wave packet case. A sharp drop of K (t) in Fig. 5.17b shows the location of the solitary wave. 5.4.3 Temperature Holes We discovered that K (t) for a two wave collision process has a dip n r r around the site of the collision (Fig. 5.18a,b). This is partly due to the concentration of the spring; the kinetic energy is transfered to the potential energy when two waves collide. Obviously this effect is inde- pendent of the type of the force function, and in fact the same thing happens in the harmonic lattice (Fig. 5.18c,d). For the purpose of removing this spurious dip from K (t) we experimentally defined the temperature in terms of the total energy. Although the temperature hole disappears in the harmonic lattice with the new definition of the temperature, it still remains in the anharmonic lattice (Fig. 5.18e,f). Hence we must conclude that the phase shift of the solitary waves is also responsible for the tem- perature hole; both waves are effectively accelerated during the collision process, contributing less to the local temperature. In the thermal situa- tion the temperature holes can occur at any location of the lattice with an i£0"O Hi IS 1.0 \ (a.) Hwnwn/c (b) Vi=.l Q ' ■ ■ ' ' ' J L 3o TL Fig. 5.17. Time averaged kinetic energy, TiirK Av*n>jfJ Kinetic Energy ■■••••>•»»•»» t-y (tt) U9 t=10 10 10 30 *0 Jin 5? Time Aferaaed Kinetic E»«j/ 10 Zo 30 40 (b) t=£0 *o 71 Time A*n»|ed V. <— | 1 JO i "\ ~* J Tottl Enegy .20 i i * (e) t=zo Fig. >.18. Time averaged energy of each atom: (arbitrary unit). JPW lattice Tme Avcnyed Kinetic Energy to XO %o ¥0 (C) 119; so Time Avenged Kinetic Energy (d) 1 -b (5.4) f(x) = 2x+b x • -b . It is obvious that the amplitude must exceed the threshold value b; other- wise the wave motion would be pure harmonic. The numerical results are given in Table 5.2 for various b and the initial velocities. A typical profile is given in Fig. 5.19. The two wave collision process was also examined for a typical value of V. and b. Two solitary waves were verified to exist even after the collision. There was no phase shift nor energy left at the site of the collision. 5.6 Hard Sphere Plus Harmonic Force Function The hard sphere plus harmonic force function first introduced by North- . (13) cote and Potts may be considered as a special case of 55.5. This model combines the harmonic lattice and a very strong nonlinearity , hence it is theoretically very convenient. We define the force function as: f (x) = x x ' -b f(x) « -« x -b where b is the collision distance (Fig. 5.20). Computationally the hard 121 Table 5.2 Solitary Waves in Piecewise Linear Force Model b V. i Amp . V .01 2 .524 1.35 .05 2 .628 1.33 .1 2 .731 1.29 .2 2 .851 1.25 .4 2 .981 1.16 .6 2 1.025 1.086 .2 1 .491 1.167 b : break point V. : initial velocity l Amp.: amplitude of the solitary waves V: propagation velocity 122 Fi . 5.L9. Solitary wave in piecewise linear force lattice f(x) = x x > - .2 f(x) = 2x+ .4 x ■ - .2 123 sphere type interaction can be realized by interchanging the velocities of two atoms when a collision takes place. We can obtain, on one hand, the harmonic lattice by increasing b, on the other hand, we can obtain the hard sphere gas model by letting b close to 0. In our numerical experiments we specifically chose b = .5. Judging from the results in the previous sec- tion we can expect the solitary wave solution in this model too. The numerical results are given in Table 5.3. A typical profile of the wave is also given in Fig. 5.21. We found that the solitary wave is quite stable, provided that the initial impulse is sufficiently large. Note that the wave motion is harmonic if the amplitude is less than b. For some critical values of amplitude the once formed solitary waves decay as they travel along. We also examined the two wave collision process in this model. Fig. 5.22 shows an example. In this example, the initial velocities are .9 for both ends. The following facts are remarkable: (1) the energy is left at the site of a collision. The amount of this energy is considerably large (~4%) ; (2) both waves are advanced after the collision; (3) the velocities become smaller after the collision due to the loss of the energy. It is notable that the energy loss in the collision process is large enough to affect the propagation velocity in this model. Presumably the same thing could happen in the JPW lattice; the nonlinear effect is too small to be noticeable. It is quite conceivable, however, that the soli- tary waves encounter many collisions in the thermal situation, if the lat- tice size is sufficiently large. Then the loss of the energy will be significant . Table 5.3 Solitary Waves in Hard Sphere plus Harmonic Lattice 124 •J- b V. l V E . sol E ., tail E . col h .5 .5 .19'" .80*' .025 .02 5 1.19 .5 .5 .90 1.0 1.60 1.82 .4043 .4999 .0007 .0001 .01502 .02 5 .025 .5 1.5 3.07 1.12500 -5 .2 10 .00248 .025 .5 2.0 4.0 2.0000 '.5 10~ 7 .025 Collision of the two solitary waves of the same size. ** Decay at n=12 Vf-'c-'r Decay at n=38, b: break point V.: initial velocity (initial energy = .5-V. ) i ■ i V: propagation velocity E , : energy in the solitary wave part sol E . • energy in the tail tail E , : energy loss at a collision col h: integration step 125 Fig. 5.20. Hard sphere plus harmonic force, v\A " \ ' -.1 -.a 20 30 —•— • r- b - -.5 V. - .9 r -.3 Fig. 5.21. Solitary wave in hard sphere plus harmonic force lattice, 126 "time n Fig. 5.22. Head on collision for hard sphere plus harmonic force model. 127 CHAPTER 6. HEAT CONDUCTION PROBLEM IN ANHARMONIC LATTICES AND ITS PHENOMENOLOGICAL INTERPRETATION We will consider the heat conduction problem in this chapter. Our main purpose is to investigate the applicability of the idea of the soli- tary waves to this problem. As we discussed in the previous chapter the JPW lattice and the hardsphere plus harmonic lattice (to be called H+H lattice) both have the solitary wave solutions. The heat conduction problem has been numerically considered by several authors and the temperature (21) (44) gradient has been observed in both types of lattices. ' We can obtain a better insight into this problem by locating these two types of lattices between two extreme cases; the harmonic lattice, the JPW lattice, the H+H lattice and the hard sphere lattice. They are put in the increasing order of nonlinearity , in a sense. The last model will need some introduction. In this model each atom of the lattice can travel freely in the interatomic space until it collides with its neighbors. The colli- sion is assumed to be elastic as is the case for H+H lattice. Since the velocities of two atoms are exchanged in each collision, the energy can propagate through the lattice as if it were carried by one atom which passes through others without any interaction. Therefore, if this lattice inter- acts with two heat reservoirs of different temperatures, the internal temperature distribution will be uniform, although there is a heat flux from the high temperature side to the low temperature side. This model is equivalent to Knudsen gas model, in which each atom literally propagates (45) from one reservoir to another without collision. The harmonic lattice, on the other hand, also has the same thermal characteristics as the hard sphere lattice. In this case each Fourier component (or the normal mode) carries the energy independently of other modes at its group velocity, hence 128 there is no temperature gradient. We also derived the same tendency for the semi-infinite harmonic lattice (§2.6). It is interesting to note this conceptual resemblance in these two models. If we are to consider the solitary waves as the media which carry the heat in the JPW lattice and also in the H+H lattice, we encounter a serious paradox: the solitary waves are conceptually identical to the individual atoms in the hard sphere lat- tice or the Fourier components in the harmonic lattice. Hence there would be no temperature gradient. This, of course, is a contradiction to the observed temperature gradient both in numerical computation and in the real crystals! We will propose in this chapter a phenomenological theory to solve this apparent paradox. Our discussion will be mainly qualitative due to the following drawback in the present stage: (1) No concrete theory has been established for the solitary waves in the JPW lattice. Therefore we must use the empirical rela- tions such as those obtained in Chapter 5. (2) The stochastic nature of the present problem greatly limits the accuracy of the solutions due to the fluctuation therein. A tremendous amount of computation time is required to obtain the numerical results which are accurate enough for the quantitative discussions . In the first part of this chapter we will analyze the numerical results for our model lattice (N=50, JPW type force function). In §6.1 we will con- sider the cooling process. The heat reservoirs have zero temperature, so that they work as heat sinks. We will compare our heat reservoirs and the viscous type heat sink introduced by Nakazawa. In §6.2 we will adversely consider the heating process of a cold lattice. Only the anharmonic lattice 129 will be treated here since our purpose is to observe the formation of the temperature gradient. In §6.3 we will consider the heat conduction process which we briefly introduced in Chapter 1. We will compare different boundary conditions. In §6.4 we will interpret the heat conduction phenom- ena in terms of the solitary waves. In §6.5 we will discuss H+H model and we will imply the importance of the collision process with regard to the formation of the temperature gradient. 6.1 Cooling Process of Anharmonic Lattices We computed the residual energy in the lattice when two reservoirs have zero temperature to demonstrate that our heat reservoirs actually have viscous effect. We also examined the viscous type heat reservoir used by Nakazawa for comparison (§4.3). In the latter case the end atoms lose their energy no matter where they are located (§4.3). The results are in Fig. 6.1. The curve I corresponds to our model heat reservoirs and the curve II corresponds to the viscous type. The initial conditions are. the same for both cases. The energy loss is faster in the latter case due to the continuous interaction. Our model may be described as having the infinite viscosity inside the reservoirs but none outside. 6.2 Heating Process of Anharmonic Lattices The purpose of this numerical experiment is to observe how the tempera- ture behaves in the early stage of the heating. The only difference between this experiment and those in the next section is the initial condition. The temperature distributions are given in Fig. 6.2. The reservoirs have the temperatures .02 and .01, respectively. At t=0 two solitary waves were produced. Fig. 6.2a, for example, is the temperature distribution at t=50. 130 Resi. • i • . i ( ure o I 2 r )t li atom Fig. 6.3 (continued). Internal temperature gradient and T^^ (Data No. 2) 137 138 -vT(* ,0 " 5 J 2-1 J L j I I i_ -i I I "time looo Jlooo Sooo *aoo (b) Temperature gradient Temperature ( * \0" y ) l.» - '*» 1.5 - 1.1 o j i i i_ j I i 1 1 i _i I 1 L 10O0 MOO 3ooo time *O0O (c) Temperature of 25th atom Fig. 5.4 (continued). Internal temperature gradient and T^ (Data No. 5). 139 a vel ocity v 1 given to the atom n=l at t=t we can find a corresponding interaction in the other end of the lattice at t=t , where t -t is close to the time it takes a solitary wave to cross the lattice with the initial velocity v„ • We used the initial velocity vs propagation velocity rela- tion in Fig. 5.6. We also found that the interaction time is not distributed uniformly but more dense when a solitary waves arrive at the end of lattice. It is clear that the solitary waves introduce interaction with the reservoirs. (6) Histograms of the velocities of the end atoms, at the time of interac- tion with the reservoirs were obtained from the above data. They are in Fig. 6.5. The distributions after the interactions are merely the sampling from the probability density function (1.4) with T=.02 for n=l and T=.01 for n=50, respectively. We could not determine, however, the nature of the distributions for the incoming velocities due to fluctuation. This diffi- culty was also reported by JPW . (7) The collision frequency is higher on the low temperature side. This difference is necessary to retain the lattice between the two reservoirs. In other words the total momentum of the lattice is balanced to zero in the (21) time average. Consequently the pressure is the same on both sides. (8) The temperature distribution showed a periodic structure. This phenom- (44) enon was also observed by JPW and by Helleman. (Fig. 6.6.) (9) The heat flux J is considerably lower in the harmonic lattice. This result was also reported by JPW and by Helleman. Since the collision fre- quency is also low in the harmonic lattice we cannot single out the nonlinear effect from this result. We can explain the low collision frequency in the harmonic lattice as follows; firstly the anharmonic lattice is stiffer than the harmonic lattice, therefore the "sloshing" effect of the lattice is less a 'To 140 u 03 o O o o II o 2 Q sis n ir t*> -1* *2 ^ 1 1 1 ^^ < _) fr — r - •< — i • c E Hm 1 r U. 1 | ' I J— 1 i 1 © TO a a o en QJ •H 4J ■H U O ,—1 QJ > O e crj n o ■U CO •r-l m ton ■H 141 -2' Tewpera-We (*I0~ ) yo n Fig. 6.6. Periodic behavior of internal temperature distribution for harmonic lattice. 142 important in the anharmonic lattice. Secondly the solitary waves propagate faster than the sound velocity, resulting in more collisions per unit time. From the above reasons we can say that the thermal contact is better in the anharmonic lattice. In the present scheme of the interaction between the lattice and the reservoirs the heat flux seems to be limited by the boundary effects rather than by the properties of the lattice itself. Especially the "sloshing" effect of the lattice is not likely to happen in the real crystals. 6.4 Phenomenological Interpretation of Heat Conduction in JPW Lattice We consider the heat conduction problem in the JPW lattice in terms of the solitary waves and the tails. First we should recall that the inter- action between the lattice and the reservoirs produce either the solitary wave and the tail (when V. is sufficiently large) or a wave packet which does not yield the solitary wave and the tail within the limit of the lat- tice size (when V. is small) (§5.1) . The numerical results in the previous section as well as in Chapter 5 strongly indicate that the solitary wave carries the major portion of the heat . Because of the high concentration of the energy and the high propagation velocity, the energy is transfered very efficiently (§5.1, §5.2). The solitary wave is not absorbed completely in the reservoirs but partially reflected (§5.2). At the same time, new soli- tary waves are created due to the interactions induced by the incoming solitary wave. Therefore the solitary wave also works as the "exciter". Because of the high propagation velocity, the collision frequency between the lattice and the reservoirs is considerably increased ((9), 56.3). An important fact is that the momentum balance of the entire lattice ((7), §6.3) is maintained chiefly by the "catch ball" of the solitary waves between the two reservoirs. 143 The tail (including the wave packets for the low initial velocities) is characterized by dispersiveness , low energy and low propagation velocity. Its contribution to the heat flux is not significant, but the individual tail gives a decreasing contribution to the temperature in the direction in which it propagates (§5.4). We should note that there is a big difference between the waves in the harmonic lattice and the tails in the anharmonic lattice, although both waves are harmonic. In the former case, the harmonic waves are responsible for the momentum balance of the lattice as well as for the heat flux, whereas the tails are responsible for none of them in the latter case. In other words, there is a constant (small) flow of the tails from each reservoir in the anharmonic lattice, the flow rate being determined by the traffic of the solitary waves. We can expect that the unbalanced flows in two directions due to the temperature difference will result in some gradient from the high temperature to the low temperature side . This is what we propose to solve the paradox which we stated in the beginning of this chapter, since both the collision process and the inter- action between the solitary waves and the disturbances do not seem to be significant in our model lattice. Unfortunately we do not have enough numerical evidence as yet to support this proposition quantitatively. The disturbances in the lattice are the collection of dispersed tails. Since the tail loses its coherency as it propagates, the disturbances are quite incoherent and do not contribute to the heat flux. The energy of the disturbances does not increase indefinitely but is maintained at a certain level on the average due to the interaction with the reservoirs. The heat flux and the temperature gradient are related to each other through the nonlinear coefficients; the propagation velocity, fraction of 144 the energy in the solitary waves and the stiffness of the lattice are all dependent on the nonlinear coefficients. The effects of the boundary conditions on such quantities as the heat flux, the temperature gradient and the thermal conductivity are not clear. Presumably the energy is more efficiently exchanged at the reservoirs if the end atoms have more freedom (i.e. m is large). Therefore we can expect higher flux in such cases. At any rate the dependency of the thermal conductivity on the boundary conditions implies that this quantity is not intrinsic to the lattice. Finally, we will summarize the results of this section. (1) The solitary waves carry the major portion of the heat. (2) The solitary waves considerably increase the collision frequency with the reservoirs, due to high propagation velocity. (3) The solitary waves also maintain the momentum balance of the lattice. (4) The tails are one possible source of the temperature gradient in the JPW lattice. For other possibilities of the source of the temperature gradient, refer to the next section. 6.5 Hard Sphere Plus Harmonic Type Model and Other Possibilities of the Temperature Gradient As we discussed in the previous section, the JPW lattice has the characteristics of Knudsen type model. We next consider the H+H type model, Helleman investigated the heat conduction problem in this model and dis- covered that the temperature gradient becomes more obvious than in JPW lattice. We showed in §5.6 that H+H type lattice also has the solitary wave solution and we discovered that a non-negligible amount of the energy 145 was lost from the solitary waves in the two wave collision. From these results we can consider the collision process as very important in the strongly anharmonic lattices. Namely, the solitary waves lose significant amount of energy as they collide with other solitary waves. In this picture the relation between the heat flux and the temperature gradient becomes inherent to the lattice. Another possibility is the interaction between the solitary waves and the disturbances. In the JPW lattice the effect of the disturbances is not important. But if we increase the number of the atoms enormously, such effect may not be negligible any more. In this case, too, the solitary waves will decay as they propagate through the lattice. No research has been done in this area, however. We will need a very high speed computer of the Illiac IV class to perform numerical experiments with super large-sized lattice. 146 CHAPTER 7. CONCLUSION In this thesis we considered the properties of the solitary waves both analytically and computationally. In the first approach we introduced the piecewise linear function model, by which most of the special characteris- tics of the solitary wave solution and also the periodic solutions can be explained with the elementary functions. As for the second approach, we performed numerical experiments on the anharmonic lattice introduced by Jackson, Pasta and Waters. We discovered that the impulse applied to the end atom of the lattice produces the solitary wave very efficiently. By this method we performed various kinds of experiments as to the properties of the solitary waves. One remarkable thing is that the solitary waves thus created lose their energy in the collision process. This has never been reported. Although the amount of loss at each collision is very small, this fact indicates the picture of the decaying solitary waves in the an- harmonic crystals. We could also show that the piecewise linear function model and the hard sphere plus harmonic type model also have the solitary wave solutions. Our numerical results were applied to interpret the heat conduction phenomena. Numerically we could show that the major part of the heat flux is carried by the solitary waves. The origin of the temperature gradient is still unclear, but the dispersive wave packets which are created with the solitary waves while the lattice interacts with the heat reservoirs might be a possible source. The collision process and the interaction with the disturbances seem to be negligible in our model lattice. Our numerical results indicate that the heat flow in our model lattice is near Knudsen 147 type. It is conjectured that the thermal conductivity will become the intrinsic quantity by increasing the lattice size considerably. Finally, we mention that the numerical integration method due to Nystrom was found to be superior to the classical Runge-Kutta type inte- gration method. 148 APPENDIX A SOME PROPERTIES OF TCHEHYCHEFF' S ORTHOGONAL FUNCTIONS Tchehycheff s equation ,1 2 ^ d ^ dy 2 (1-w )2" W dt +ny = (A<1) dw . . (25), (26) has two solutions and y = T (w) = cos(narccos w) first kind (A. 2) y 9 = U (w) = sin(narccos w) second kind . (A. 3) Note that the above solutions differ from the usual solutions by the factor 2 (according to van der Pol ). T (w) is a polynomial of n order, whereas U (w) is the product of Vl-w and a polynomial of n-1 order . Explicit solutions are given as follows for n = ~ 4 : T Q (w) = 1 U Q (w) = T (w) = w U (w) = Vl-w T (w) = 2w 2 -l U 2 (w) = 2w7l-w' 3 2 7 2 T (w) = 4w -3w U 3 (w) - (4w -l) v / l-w L. 9 3 I 2 T (w) = 8w -8w +1 U 4 (w) = (8w -4w)Vl-w . Note that T (w) is an even function of w when n is even, whereas U (w) n n is an odd function oi w when n is even. TT Since arcsin w = — - arccos w, cos[ (2n+l)arcsin w] = cosf. (n+2)n- (2n+l )arccos w] = (-1) sin[ (2n+l)arccos w] " (' 1)nU 2n + l (w) ' (A ' 4) 149 The relation between Bessel functions and Tchehychef f ' s polynomials is best explained by using Hansen's integral representation for the Bessel (25) function of integer order : • " n n -4- ft , X 1 f-ltCOSO Ann • . r- N J (t) = —— J e cosn9 d9 . (A. 5) o By introducing a new variable w = cos9, . -n -1 .. /x i f itw r -1 -i , . Ax-1 , J (t) = e cosLncos wj(-sino) dw n TT « .-n 1 T (w) . 1 ? n iwt, .. ,. I ■ e dw (A. 6) ! -1/7-2 1-w The above relation shows that the Bessel function has no higher frequency 2 components than w > 1. This corresponds to the fact that we are dealing with the low pass filter (§2.1). Furthermore , - J V~ = ^(J ,l(t)+J l(t)) t zn n+1 n-i . -n+1 1 iwt = —: f (-T ,(w)+T . ) ) 6 dw 2nn J v n +l v ' n-l y/ / ~ •J 1-v -w . -n+1 1 1 iwt I U (w) e dw , nn *J n v since T . (w) - T (w) = cos[ (n-l)arccos6]-cos[ (n+l)arccos9] n-1 n+1 = 2sin (narccos9)sin(arccos9) = 2U (w)-U n (w) n L -ifi ' w U (w) . n 150 Hence, J (t) 2 .l-n 1 — = z~ U (w)cos(wt)dw t nn CJ o u > CO CO CO » — ' ^> > > CN CN X •*■•> ^ S~^ M s~\ X /— \ X /— N X X « ^ •> ^ »> ^ « ^ «■> ^5 + + N « N « N ^ N ^ N « CN CN 00 N 00 N GO N 00 N 00 N , S S— - v — ' GO " — ' GO v — ' GO s — ' GO ^— ' 00 CJ CJ CN ^-^ CN 1 — ' CN " — ' CN •^s CN " — ' 1 1 c CN C CN c CN c CN c CN cO X co c en c CO c CO c CO c ^~s ■~~s u o CO X a CO TJ -a CO CJ CJ CO X X CO II II 1 i U o i i X CJ 1 i X) XJ 1 1 CJ cj 1 1 T3 TD CN CN X CH l 1 cO CO l 1 U X 1 1 X CO i 1 cO X 1 1 N ^ >w> N 3 X ^ * *• " M N 00 ,Js! y : •— V / — ., •"" > C " * ^ ,~ N X •-N ^ / — s ^ 1 N N ^ X ~ ^ «■• ^J u ! OO GO N rv N ^ N ** CQ ^— N v — ' — ' GO N 00 N 00 N cO CQ CN CN ■w* 00 ^^ GO ^^ 00 1 1 C C CN ^ CN "•— ' CN ^-^ X <: en CO C CN c CN c CN CJ i 1 u X l 1 X CJ CJ X X S-i L ~' 1 1 1 1 1 i cd V ^ •^\ ^ CN XI XI T3 X X 1 1 + 3 J3 3 CN * ^ ^ ^ x CJ o u u 3 CJ CJ 1 V ■-> 1 ^3 CJ 3 i 3 3 X) /~\ /— s s-~^ ^^ ,^-s •— N 3 3 | X 3 J3 X X X *> t X I 3 ^3 3 CM ^ ^-^ > ~ / ^"^ i T: t 3 3 3 3 1* 1 03 cO cfl I cO cO i ! 152 APPENDIX C d 2 ON THE TWO SOLUTIONS OF (—) = P(w) The differential equations which we are interested in is of the following form: , 2 O - P(w) (C.l) where P(w) is a cubic or quartic polynomial. We showed in Appendix B that the solutions of the above equation are expressed in terms of the Jacobian elliptic functions, sn, en and dn. Corresponding to the fact that the circular functions sine and cosine are both the solutions of ,dw 2 ( dz } = 1_W with different initial conditions at z=0, there are also two solutions to (C.l) with different initial conditions at z = 0. These two solutions can be matched by shifting one of them by a quarter period, as is the case for the circular functions: TT cos z = sin(z+ — ) . For example, , 2 (j 1 ) = (a-w)(w-b)(w-c) dz has two solutions, / u \ 2 ,va-c , . .2 a-b w = a-(a-b)sn ( — — z k) , k = f 2 a-c and , a-b 2 , \a-e . b - c s n (-* — r— z , k ) a-c L w 2 = — l-(— )sn (*-y z,k) 153 By using the following formula (47) sn 2 (x+K) = ^Sl = lZ§ f J f dn (x) 1-k sn (x) we can show that W L ( z + ~ < k ) = a-(a-b) Va-c 1-sn , a-b 2 1- sn a-c a (a-b) 2 2 a- — * sn - a+b+(a-b)sn a-c a-b 2 1- sn a-c , (a-b)(a-a+c) 2 b- J " L sn a-c , a-b 2 1- sn a-c a-b 2 b - c -a^c~ Sn , a-b 2 1- sn a-c w 2 (z,k) . Hence we can choose one of the solutions to (C.l) to discuss the nature of the periodic solutions without losing generality. 154 APPENDIX D COMPUTER PROGRAM r i ^e i mil/xmv 155 tjF. CjIN t% THIS PROGRAM COMPUTES THE. HEAT FLU* ANT TEMPERATURF GRADIENT *% OF PME UlMENSIUNAL AmhARMIJNIC LATTICE wTTH FORCE FUNCTION t% Fs**ALPHAX**2 ♦ RI*X**3 X(AND IMPACT TYPE HEAT RFSERvOIRS AT bOTH ENDS.* %% wrITTEn 2/24/73 Ano rEvIsl) 3/25/73 k«mI|JrA FILE STATK)N(KlNF) st 3#MA'HZSCL'H2 , H4»HIVHl5»ZMXJ REAL M»*5. Kft, K7. K«, K9.H. I3» 16. 19. I192.HI19?»hM23» H«o* HEAL AMP,wnTH,SUM,TMP,RFLUX.LFLllX,TMH.YiPl,Yli,Zjl,H?i3» REAL ARHA* Y,Z[0» 5 1] . TP . T<5 [ 1 « 50 . « 5l]» REAL ARHA Y K 1 , «2 , K 3[ 1 : 5 1 ] . K NTp[ 1 I 50) , LF. RFT * » ^ 00 1 » INTEGER AKRAY XC t) » YCflt 1 « 1 00 J J INTEGER LCNTR.RCNTR.TW035' INTEGER I.J,K,L.H,JJ.KK,LL»L1.FLG,FLG?.FLG3, INTEGER X; XFOR RANDOM NUMRERS INTEGER M.TIME.STP.ISTAHT»ILAST»HI> HEAL P H Xt=SECO (271 REAL P Reg EN n PRUCED VAL INT REA REA %%%%t%%%%ii^%%% RE RJCE EGIN NUrtO U26l Ni)» *i ■** ROCE IN *UNE * REAL TMP« MaXW ; % Hit, URE UE EGER L aR L %%%* GIn REAL INTE REAL i%t%%t%%%%%%ttt%ti*%l%%%%ii*l%%%%%t %%%%%%%%%%% UiIRE RNDXj l%ti%%i.lZ%%%tt%l*%i%%%%%%%%l%%%%t%%%t%%%%1%%Xt <* OUE TU D. K NUTH R0(0lNTEGER(3iai592653 MUX X)+OTNTFGfR «29)) M ni) T w 035» RNUX:»X/T^035; *RNnX X* *%%*%% tt%%l%%it%t%*%%%%%%t%t%l%l%%%%X%t%%%%%ti DURE MaXWelL(T)' REAL T» l%%t%%t% *%*%%%*%%%*% 1*%%%%%%%%%%%%%%%%%%%%%%%% OIMFNSIONAL MAXWELL-HOLTZMANN OlSTRlROTlON TMp: =-T*LN( 1-RNDX) » ELL,=S«RT(TMP + TMP); %%% j IS ALWAYS ASSUMfD T * mAxwEll"ROlTZmAnn i% i*,t>t%i%%%*t%%%*%%%*%%%% LSS^fA.IS.lL»INTL»lCPT»SLHP); IS.IL.INTL; IS, IL, INTLj RAV A[* J ; I C P T » S LO P J i%)k*%%%%%%l%l%%%lit%%%l Tl J l» E H I . J I SUMY, SUMZ. MAX. TEMP; SUMYt*SUMZ:»MAX:sTEMP|sOJ OOOO0100 00000200 00000300 oooooaoo 00000500 00000600 00000700 noooottoo 00000900 OOOO1O00 OOOO1100 00001200 00001300 OOOOHOO OOOOl^OO 00001600 00001700 00001800 00001*00 00002000 00002100 OOOO2200 00002300 00002400 00002500 00002600 00002700 OOOO2800 S***X**OOOO2900 00003000 ***s***oooo3ioo 00003200 00003300 oooo3aoo OOOO3500 00003600 *«**»**noOO3700 OOOO3800 *XXXX**OOOO3V00 OOOOUOOO 00004100 00004200 OOOO4300 00004400 RE POOOOO4500 00004600 OOOO4700 00004800 00004900 00005000 00005100 00005200 00005300 00005«00 OOOO5500 00005600 OOOO5700 00005800 156 fur i» = is step intl until il do begin J'=I/INTL' SUMYl.SUMy+I; TEMPJ=TEMP*1*A[ j j i SUMZ««SUHZ*Al Jl I MAXt=MAX*I*I» ENDj j:=( I| - I S)/ I NT L + 1 I II l=J*MAX"SUMY**2j , Kpt:s(sumZ*max-sumy*temp)/ti* S*«lNlTERCFPT pRncLnuw VALUE JNTE^FR INTEGER ^eg In FORMAT I) F 2( X« ) ,.f- F3?(# n " a » f 9 8 t 4"A6A7AH K « PI «PlJ I NT Ql=pUINT REPLACE WRITL(ST IF S C Al E RE'ilN SCALE WRITE(ST PflR I I «0 IF I EM„ F HAZE N.MS; fJ.MS; AHKAY SLnP,*(J*TE«P-SUMY*i)UMZ)/Tl| *»SL()PE Up L^ A ST SijyARE pITTINt,; l.PLUTCXCIl,YCC3»N#MS)» XCi), YLf)[ * i t UT F0U"A1 1UA1 1 |)Al 1FA1 1 1040200") . » > »U' ••»•*;• • • # t < » 4 > • fi , x 6 4 , n - n ,a M oo M ). ii,lJ(i«i'f«j"),a».()0 , »). ( j3»X2)'a"Ali200 , '>» GtH SCALE*- Lf,J INTEGER MARRAY X r o « ] » p. q; RaY ASCII 2u3J72n2E2F16US?SOROCOf)OEOFlOlll?133C3o32 3F2»'lClL)lEl f "** I """*J»8' ( )*♦•"./ 12345678 AriCUEKr,HljKL M N(lpQRSTlJVwXYZ[ ,,a, 'EO"«)-,'' u .a i<4..ai.M<.H/UMMoa.o..u->u.. u. OiUrouou.. . H i .dj(l (4 8b06rt< r M«tt99i9 i ' STEP 1 UNTIL 2A OU HEtilN m U L) b = U T H L -j tfRlTE/H£SCL> E WHI TE(ST AT IUN.F2) i END! TATIUN, f3 ) > tat I un • r a »FUN l»«o step 5 until 60 11 in TATIlJN.<4"AilFAllFU0">>l TATIUN. <4"Al 1 200 n > ) I TATIUN.< n TlML«".T e >.X3."ENRc.Y* H .R11.5.X3» ,, ALPHA« , '.En.4> MS»SuM. ALP H D ) J 1 STEP 1 UN UL n no YCUl I ) t*YCUl 1 ) + \S I REPLACE P HY u*XCurIJ EUR l.u+YCOII] EUR 1J 00005900 00006000 00006100 00006200 00006300 0000*400 00006500 00006600 00006700 0000* tt 00 00006900 00007000 00007100 S**X*X**00007200 00007300 00007400 00007500 00007600 * % * % % I % t 7 7 00007800 0000/900 00008000 00008100 00008200 0000*300 00008400 00008500 00008600 00008700 00008HOO O0O0A900 0000 9 000 00009100 00009200 00009300 00009400 00009500 0()009ftOO 00009700 00009800 00009900 010 00 00010100 00010200 010 3 010 4 00 010 50 001 06 00 00010700 00010800 00010900 0001 1000 0001 1 100 0001 1200 0001 1 300 0001 lion 00011500 0001 uoo 0001 1 700 0001 1800 0001 1900 00012000 157 WRITHSTATIUN* liX[*l)l 0001210C end» 00012200 write* stat ion. <4"a112a11900">>» 000 12 300 END i)F Ha-?ELTINE; 00012^00 PROCEDURE KINETICTEMP(N»X,TIME.N1 >' 00012600 value NtTiME.Ni> 00012700 i NTEtiLK N.TIME.N1I 0001280C HEAL AHRA Y *[*]> 0001290C *I«I*I«J»X*X%*«*******»«X**»**«»*4«»*4»*i»I«»X^8*I«Hil-«I«I«««I«««*«lO00 13000 fit&TN FORMAT out F51 1X5, w KINETIC TEHPERATORE AT". " TI*E «".14)» E^2( x20»" I «. 5(9(".")»" | "). ) » E56(xl'El3.6»X2M2»X2'"0 rt ')» F54(Xl.El3.6.X2.l2.X2,"o".*<''X''),), FS7(X20»"0"»Xtt»". FUR I J3 1 STEP 1 UNTIL N 00 TAXiBrtAxCCFUURlEKtl] l"X(I ]/TMH)*T*X)) wRITE(LINE »f 51 . TIME ) * WRITE(LINE .F57); *RITE(LINE »E52); Tlj=100/TAx; FUR !» = * ST E P * UNTIL N* nO BEGIN K: S T1*(TEMP: 3 F0URIER[I1 )l If a NEO THEN WRITE(LINE »F5a»TEMP»I.K5 El SE WRITF(LINt .F56.TLMP. I )> ENO, *7); rtRITE(LINE.'LFLUX»KFLUX) J l.SSQ( FOURIER* 1 . N . 1 » T 1 . T 2 ) > *RlTE(LlNE,<"SLUPL=".E12.5.X3."LEFT LSO TMP»". E12.5, X 3,nRl bH T LSQ TmPc".E1?,5>.t2,T1 + T?, Tl*N*T2)l WRITE(LINE» , LFL'Jx* . 5 • 00013100 00013200 00013300 OOO13A0C OOO1350C OOO1360C OOO137O0 00013800 0001390C 00014000 oooiaioo OOOH200 OOO143O0 OOOlA^OO OOOH500 0001460C OOO1470G 0001480G 00014900 00015000 00015100 00015200 00015300 0001540C 0001550C 0001 56 0t O001570C O00l580fc 0001S90T 000 16000 00016100 00016200 00016300 00016400 00016500 00016600 00016700 OO016800 00016900 00017000 00017100 OOO17200 0001^210 OOO17300 oooi7aOO RFLUX*.5,LCNTH.RCNTR.T2, FOURIER [25]) I END; «*XKINETIC TEMPE-RATUKL %%%%% %%%t%%%%ti,%%%%ti1,ii,%i%%%t%tii,%tn%%t,i,%i,%%%%%%i%%%t%%%%%%l%%%%t>%%%%%%%%%noOl750Q PROCEDURE PnLYGRAPH(P. N.w. SCALE) > 00017600 VALUE P.N, SCALE; OOO17700 INTEGER P,N. SCALE; OOO178O0 REAL ARRAY W[*.*l; 00017900 %%%*%%%% t *%%%%%% it tt%t.%%t%%%%%%%ttt%i%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t 000 I SOOO CqMMEnT THI S PRqCEDURE OI s PLAY s 10 GRaPH s U s InG A LINE 00018100 158 begin FflHM, printed. they shuuld not overlap. w contains the OATA.N IS THE NUMtiEH OF GRAPHS. P tS THE STARTING POINT UF ARSCISSA.j F«( In In wKITF.d AT OUT FO( x 4. 12( 13. x 7) ) . »« F2( Xl » l4»Xl . W( x* ."*")) i FJ(X6»il( ,, X , '»*("."))» , 'X ,, )» FFF ( I3.10CF12.8). I«)» [ii lU(Fli.lO)* I4)« F2(Xl.I4,Xl»X*»".".2(X*. n *"»X*."«i".X*. *•♦''. X*» "()"»X*'"»"))J TEtiER t»J#K.TMP»L' rt*»t"R ARRAY VdUOJi lNt,<« v l».,xl2»"v2 ".xll."OUAu M »Xin,«CU n E»,X' n » "SCALE xio , "stfpsUe">) ; "HlTE(LlNL.<6(Rll.a#X3)/>,AMP.rt0TH.ALPH0.HI,SCALF»H)! FUR l«sO STEP 1 UNTIL N 1 V 10 -1 DD BEGIN «RITE(LINE»FU'Ff)R JJ=0 STFP 1 IJNTTL 10 DO ►*♦ I * 1 0* J J » *RITE(LINE,F1 ); F OR J:s l STEp 1 ijNTI L N D" BEGIN FOR K l « 1 STEP 1 UNTIL 10 m IF TMPlsIMO + K LE« N TMFN Yt K ] : = w[ T M p. J)*SCALF» VRITE(LINL'F?»J .U'^Ytll-' + Yt^TYn]' 9*Y13 J-Y12 J»«+Yt4]-Yl1]»9+Y[5]-Yt4l« Vvyt 9]- Y l8J . V*y[ I (J ] "Y T V 1 ) S EN D , ^RITE(LINE.FJ) J «RITE(LINE tSMP ll.FO'FlIR JtxO STFP UII (P* J *10*.l) ) J FUB u»0 STEP 50 UNTIL N I V 50-1 On. MEr, IN wRITL(LINE»»FnR L I = 1 UNTIL 10 D0(P*l*10*L))> K)B Lj«l STL^ 1 UNTIL 50 no W H ITE( LlNf «FFF» J*L»FnR KJ = 1 e; T E P 1 IJNTI lu un wl k*i *in. j*l1 » J*L } ' WRMECLlNE.Fa .J*5l.FU» Kt»l STE P 1 UNT 10 0" rtlK*I*10. j*5l ] , j*5l ) • «RI TE(LINECSMP l J ) » Eno ! ENDJ ** e.no; *) J t UNTIL 11 STEP 1 tuss 1*4* ***i Eni>1 BEAL PRnCE HEGIu A*************************************************** UUHE Pi)T(X)| VALUF Xj REAL Xj EiMD ; PJT»*X«X*(.5*X*(ALI3*Bia*X)); 00018200 018 3 00018400 00018500 00018600 00018700 00018800 000l89()0 00019000 00019100 00019200 no01 9 300 00019400 019 5 0001 9600 0001 97 00 Q0019H00 00019900 00020000 O002O1O0 O0020200 D0020300 no 020400 0020500 00020600 0020700 00020800 00020^00 00021 000 00021 100 0021200 00021300 00021400 00021500 O0021600 00021 700 218 00021 900 0022000 2 210 00022200 022300 00022310 00022320 0022**00 0022500 0022600 0022700 ***00022800 0022900 * 1! t 2 3 00023100 00023200 00023300 O0023400 XIIOOO235O0 00023600 X * X 2 3 7 00023800 023900 00024000 00024100 IL 159 %%%%%%^ y %i>1>t%i>i*i>t^1>ti*%%*%X%t%%**%t%%%%%%%%%*%%%%r%t%%%%%%%t%tt%%%t%%%%0002'i200 HKOCt'UUNfc- NYS^(N»Y»Z)) 0002*300 VALUE N» 0002**00 INTEGER M ; 0002*500 REAL AHKaY Y,^[0]» 0002*600 %%%t%%*%1,%*it%*%thl%%%%%%%%%t%%t%%%%%%%%i%%*%%%%%%X%%%%%%X%%%%%%%t%%%%%0002(K9;sV(-(YTl« = YlPi) + {YlPll 3 Yri ♦ n ) )); KMsKbj XK6xYUI-l] K5l»YH*«4*((Zlll»ZCl3)*.?*K4*H)*H»*l ST IT K21I Jl»V(K5-K6) J X*K? TEMPORARY USE FOR *8»=K7> X Y2CI-11 k7 I «yIUh2*(ZI1/J*h*k'»/ 9 ) J X2N n ITERATION] K3t I 1 «*V(K7-K8) » eno; k b » s o ; K2[l ] Jsk3[1 j:a K 2[N+l ] jsk3[N + 1] ,sOj FUH I « * 1 STEP 1 UNTIL N 00 BEGIN Y I 1 * = k2[ T] :=k2[ 1+1 ]-k2[ I ] j K3tlJ!»K3[l*l]-KitI]> K6«=KS> %Y300028*00 K5l«Y[I]*,80*(ZtI]*.20*((?Ml«Kl[Tl)*Y!l)*H)*Hl00028500 Z[ I ] »=**H M 23*ZI1/I92j 00028600 Kl t I J J = V(Kb-K6) J *K1 TEMPORARILY USED FOR V00028700 LNll * 00028800 K 1 C 1 1 !*X1 CN+1 ] :*0; 00028900 FUR I:=l STEp 1 UNTIL N 00 00029000 BEGIN 00029100 KfttsKl [1*1 ]-Kl [ I ]> 00029200 Y[I)| 3 * + H*(Z[n + H*(75*(Yin = K2tTn 00029 300 -(K6tx2^*K3[l))*25*K4)/l9?)| 00029*00 Zm»=**H*(i25*(YIl+K4)-(K6*K6 + K6))/i9 2 ! 0002 9 5oO ENDJ 00029600 YUm-M J l»YtN] j 00029700 ENOJ XNYSTROm FREE HUUNQARY REVISED 00029800 %%%%%%%it*iit%%t%-i,ii.lH%%i%%%i%%t%t%*i>%i%%%%*i%%%%%%%%%%%%%%%%t%i%%%%%%%000299Q0 %*%%%%%%% t%%*t%**%*t%il%%l%l%%%%%%%%l*%%%%%%%%%%%%%%%%%%%%%%l%%%%%%%%%%XOOO 30000 00030100 WRITL(STATI0N. ) I REAn(STATIl)N./.J.AMP.wnTH,ALPHO.RI»TIME,STP.TMP,SlJM.FLG.ISTART. ILAST.h.Kl&3.L1»FlG2,hZsCL)> WRITE( STATION. )! wRI TE 'N»AMP' 3.Ll .Sum, t mP) . K*XX*****i.CUNSTANTS i t,%%t t it Hi 1 1, 1%%1 1 I Al I 3 « = ALPHu/ 3 » K 1 4 « ■ I? I / 4 I HM^J * B M*2 } ; HI « =1/H» HZ »=H«?J * * * H INVERSE Twiiib I mi** lb | tiX**Ji***/. it******4rilTlAL C UN i)I T 1 DNS I l t t l I it k t%%%%X%%Xt%l% ll%t%% I FOR I««l STEP 1 lmTIl 1000 ()() HI3I«RNI1XJ XHIDLING (JF /MX! t ( AMP*H|)TH) *5UM» IF STP E (j L 1 THEN t THERMAL SITUATION mjr 1 1 »i step i u,nt II n do HF.<-,IN Y[ I ] » = TMP*(RNOX-.S) I ZMllalF RND* GFU .5 THFN MAxWFlL(ZMX) Else -mAxwEllc zmx ) j End ELSE * transmission property HEuIN Z[ 1 ) ' ,AMp; Z[ N ] •«-«OTHI ENnl YC N*l 1«»Y[ N] » Marti FDR I I «t STEP 1 UNT IL N DO XCU[ I ] l»l*6| L« I IF n < MMt THEn HEbIN ooo3oaoo 00030500 00030600 00030700 00030800 00030900 00031000 00031100 00031200 00031300 00031^00 00031500 00031600 00031700 00031800 00031900 ^0032000 00032100 0003??00 O0O3?3OO 00032400 00032500 000 32600 00032700 00032800 00032900 00033000 00033100 00033200 00033300 00033400 00033500 00033600 00033700 00033710 00033800 000 3 3 9 00 034000 00034010 00034100 00034200 00034300 00034400 00034500 00034600 00034700 00034800 000349 00035000 00035100 00035200 00035300 00035400 00035500 00035600 00035700 00035800 00035900 00036000 00036100 00036200 00036300 161 COMMENT CuMMpNT CUmmFnT KKJ=K/50+l t WRITECLINEISKIP l3>» W HITE ( LINE.< X 20»"C0LLISIUNS BET W FEN TIME»"»I*« « AND TImEs", U,//x4."( LEpT)" . X4 . " T I mE" , X5 , "VELnClTY" , X25>'"(RlGHT)'"X6'"TtME' , 'X«i ,w vELnCTTY* ,/ <23,"' , »Xl0.fT0>' , //> ( K» K + 50) j rUR jtxl STEp 1 u NT k SO „0 h e r, i n Fnrt JJ» = 1 STE^ 1 UNTIL HI DO «EGIN IF JJ MUU LI E«l THEN REGIN IF Y[i]+H*(ZMXl a Z[l])<..oo01 THEN «E(iIN LCNTR1»LCNTR*1 ; Z[ 1 J J=M A XWp|_L( AMP) > WRlTE(LlNE,FlF,K*J*H,ZMx,7tl>» LFLUXJ=LFLIJX*(ZI1 *ZMX)*(ZT1 -ZMX)» ENOJ IF Y[N]+H*(ZMXl*Z[N] >>,0001 THEN BEGIN RCNTR JaRCNTH+1 : zin = ZIN J »="MaX»«ELL( WDTH) J RFLUXl»RFLUX-(Zn *ZMX)*(Zll -ZMX)» C U M M FN T WRlTE(LlNE»FRF,K + J*H#ZMX,Zini ENnj EN n Uf COLLISIONS *IT H H^ AT RESErv^^I ■JYSalN »Y»Z>J FOR I:=l STEP 1 UNTIL N DO K NTp[ I ] :=*♦/[ I]**2; end of jj loopj if flg eol 1 then 'print y and z begin Y( N-f 1 ] tsY[N] ; SOMjsO; ZMX J=0{ FOR 1» = 1 STEP 1 UNTIL N DO BEGIN *ii «»tQ[ Jiii * = z r 1 1 » Sum;s**,5*( zn*zi UZmx + czmx »=poT(TPr j» n i Yl 1 + 1 J-YC I 1 )) ) J ENOJ TP[J.NJI»Y[1J> TPC J»N+1 3 JcYtN] X TQ[ J.N+1 1 JsSUMJ X TOTAL ENFRGY END Of SAVING DATA r OR PRlNTlNr.J IF FLG2 N E « ™E N I F J M0D FLr >2 F 0L BEGIN XX HAZFLTINE FOR I»=l STEP 1 UNTIL N On YCOt m ( Z[I ]*HZSCL> HAZELPLUT(XCO. Y CO»N,K*J); END HAZELTINEJ END OF J LOOPJ TH E N 00036400 00036500 00036600 00036700 00036800 00036900 0003 r 000 00037100 00037200 00037300 00037400 00037500 00037600 00037700 00037800 00037900 00038000 0003*100 00038200 00038300 00038400 00038500 00038600 00038700 00038800 O0038900 00039000 00039100 00039200 00039300 00039400 00039500 00039600 00039700 00039800 00039900 00040000 000*0100 00040200 00040300 00040400 00040500 00040600 00040700 00040000 00040900 00041000 00041 100 00041200 00041300 00041400 00041500 00041600 00041700 00041800 000«1 9 00 00042000 00042100 00042200 00042300 000*2^00 00042500 162 LF[KK j »=LFLUX; «FIKK ] JsRfLUX J If FLG3 NEQ THEN IF K MOD FLG3 EQL THEM KINFTICTEHP(N.KNTP.K*50 .N ){ IF FLG EQL 1 THEN SPNINT Y AND Z » E G T N polygraphs, n .tp.istart)! % coordinate polygraphs, n . t« . i l ast ) > **velnctty eno» IF (jl=K+5(>) MOO 1000 EQL THEN hEg i N WHI TLCLINF (SKIP i 1 ) I w *ITt(LlNE..J)J pl)K L,=100 STEp 100 UNTIL J-500 DO HEUl N LSStHLF.L. J ,50»SUM,ZMX) J LSS^(RF.L.J ,b0.SUM,Zi i ) • "KITE (I. INE.<2(Xb»I5.X2.Ell.a)>»L.ZMX* t S, L*2ll*t 5} J END (tF FLUX LEAST LINLaW F I T T I N G J w»lTE(STATIUN.)I KEAD(STA1 IUN./.M) J lp M NEy 1 T H EM GO TU Lb» tNDJ **** MOO 1000 K tiK + Su 1 GU TO L4; EN () |gft*s* lip k L00P(50 SEC)) L5» END. 00042600 00042700 0004? 8 00 00042900 00043000 00043100 00043200 00043300 00043400 00043500 00043600 00043700 00043^00 00043900 00044000 000^4100 00044200 00044300 00044400 00044600 00044700 00044*00 0044900 00045000 00045010 00 045 100 00045200 00045300 00045400 00045500 00045600 00045700 00045600 00045900 00046000 00046100 00046200 00046300 00046400 00046500 00046600 163 LIST OF REFERENCES (1) von Neumann, J. "A New Numerical Method for the Treatment of Hydro- dynamical Shock Problems", AMP Report 108. 1R (1944), also in Collected Works , Vol. 6, Macmillan, New York, 1963. (2) von Neumann, J. and Goldstein, H. H. "One the Principles of Large Scale Computing Machines", in Collected Works , Vol. 5, Macmillan, 1963. (3) Fermi, E., Pasta, J. and Ulam, S. "Studies of Nonlinear Problems I", Los Alamos Scientific Laboratory Report LA-1940 (1955), also in Collected Works of E . Fermi. (4) Zabusky, N. J. and Kruskal, M. D. "Interaction of Solitons in a Collisionless Plasma and the Recurrence of Initial States", Phys. Rev. Let., Vol. 15, pp. 240-243 (1965). (5) Russel, S. "Report on Waves", Brit. Ass. Rep. 1844, also in Lamb, H., Hydrodynamics , Dover Publications, New York, 1945. (6) Lord Rayleigh. "On Waves", Phil. Mag. (5), pp. 257-279 (1876), also in Scientific Papers , Vol. I, Dover Publications, 1964. (7) Boussinesq, J. "Theorie de 1 ' intumescence liquid appele' onde solitaire ou de translation, se propageant dans un canal rectangulaire", Compte Rendu des seance de l'academie des sciences, June 19, 1871, pp. 755-759, also "Theorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en communicant au liquide contenu as ce canal des vitesses sensiblement pareilles de la surface au fond", J. Math. Pure. Appl., Vol. 17, pp. 55-108 (1872). (8) Korteweg, D. J. and deVries, G. "On the Change of Long Waves Advancing in a Rectangular Canal and on a New Type of Long Stationary Waves", Phil. Mag. (5), xxxviii, pp. 422-443 (1894). (9) Kruscal, M. D. "Asymptotology in Numerical Computation: Progress and Plans on the Fermi-Pasta-Ulam Problem", in Proceedings IBM Scienti fie Computing Symposium : Large Scale Problems in Physics , Dec. 1963. (10) Toda, M. "Wave Propagation in Anharmonic Lattices", Jour. Phys. Soc. of Japan, Vol. 23, No. 3, pp. 501-506 (1967). (11) Toda, M. "Waves in Nonlinear Lattice", Prog. Theor. Phys., Suppl. No. 45, pp. 174-200 (1970). (12) Tuck, J. L. Unpublished Numerical Calculations at Los Alamos Scien- tific Laboratory (1961). (13) Northcote, R. S. and Potts, R. B. "Energy Sharing and Equilibrium for Nonlinear Systems", Jour. Math. Phys., Vol. 5, No. 3, pp. 383-398 (1964). 164 (14) Ooyama, N. , Hirooka, H. and Saito, N. "Computer Studies on the Approach to Thermal Equilibrium in Coupled Anharmonic Oscillators II. One Dimensional Case", Jour. Phys. Soc. of Japan, Vol. 27, No. 4, pp. 815-824 (1969). (15) Ford, J. "Equipartition of Energy for Nonlinear Systems 1 ', Jour. Math. Phys., Vol. 2, No. 3, pp. 387-393 (1961). (16) Ford, J. and Waters, J. "Computer Studies of Energv Sharing and Ergodicity for Nonlinear Oscillating Systems", Jour. Math. Phys., Vol. 4, No. 10, pp. 1293-1306 (1963). (17) Jackson, E. A. "Nonlinear Coupled Oscillators, I. Perturbation Theory; Ergodic Problem", Jour. Math. Phys., Vol. 4, No. 4, pp. 551- 558, and "Nonlinear Coupled Oscillators II. Comparison of Theory with Computer Solutions", Jour. Math. Phys., Vol. 4, No. 5, pp. 686-700 (1963). (18) Pasta, J. R. , Metropolis, N. and Bivins, R. L. "Nonlinear Coupled Oscillators: Modal Equation Approach", Los Alamos Scientific Labora- tory Report LA-4934 (1972). (19) Payton, D. N., Ill, Rich, M. and Visscher, W. M. "Lattice Thermal Conductivity in Disordered Harmonic and Anharmonic Crystal Models", Phys. Rev., Vol. 160, No. 3, pp. 706-711 (1967). (20) Nakazawa, H. "On the Lattice Thermal Conduction", Prog. Ther. Phys., Suppl. No. 45, pp. 231-262 (1970). (21) Jackson, E. A., Pasta, J. R. and Waters, J. F. "Thermal Conductivity of One-Dimensional Lattices", Jour. Comp . Phys., Vol. 2, No. 3, pp. 207-227 (1968). (22) Brillouin, L. Wave Propagation in Periodic Structures , McGraw-Hill, New York, 1946. (23) Hemmer, P. C. "Dynamic and Stochastic Types of Motion in the Linear Chain", dissertation (unpublished). (24) SchrHdinger , E. "Zur Dynamik elastisch gekoppelter Punktsysteme" , Ann.d, Phys, Vol. 44, pp. 916-934 (1914). (25) Watson, G. N. A Treatise on the Theory of Bessel Functions , Cambridge University Press, London, 1966. (26) Van der Pol, B. Operational Calculus Based on the Two-sided Laplace Integral , Cambridge University Press, London, 1964. (27) Uhlenbeck, G. E. and Ornstein, L. S. "On the Theory of the Rrownian Motion", Phys. Rev., Vol. 36, pp. 823-841 (1930), also in Wax, N. (Editor), Selected Papers on Noise and Stochastic Processes , Dover Publications, 1954. 165 (28) Saito, N. and Hirooka, H. "Long Time Behavior of the Vibration in One Dimensional Harmonic Lattice", Jour. Phys. Soc. of Japan, Vol. 23, No. 2, pp. 157-166 (1967). (29) Teramoto, E. "Heat Flow in the Linear Chain of Harmonically Coupled Particles", Prog. Ther. Phys., Vol. 28, No. 6, pp. 1059-1064 (1362). (30) Scott, A. C, Chu, F. Y. F. , and McLaughlin, D. W. "The Soliton: A New Concept in Applied Science" (to be published). (31) Zabusky, N. J. "A Synergetic Approach to Problems of Nonlinear Dis- persive Wave Propagation and Interaction" in Nonlinear Partial D_if_- ferential Equations , Academic Press, New York, 1967. (32) Wadati, M. and Toda, M. "A Soliton and Two Solitons in an Exponential Lattice and Related Equations", Jour. Phys. Soc. of Japan, Vol. 34, No. 1, pp. 18-25 (1973). (33) Zabusky, N. J. "Solitons and Energy Transport in Nonlinear Lattices", Computer Phys. Comm., Vol. 5, No. 1, pp. - , (1973). (34) Henrici, P. Discrete Variable Methods in Ordinary Differential Equations , John Wiley and Sons, Inc., New York, 1962. (35) Waters, J. F. "Methods of Numerical Integration Applied to a System Having Trivial Function Evaluation", Comm. A. C. M. , Vol. 9, No. 4, pp. 293-296 (1966). (36) Gear, C. W. Numerical Initial Value Problems in Ordinary Diff erential Equations , Prentice Hall, Englewood Cliffs, N.J., 1972. (37) Ralston, A. A First Course in Numerical Analysis , MacGraw-Hill , New York, 1965. (38) Gills, S. "A Process for the Step-by-Step Integration of Differential Equations in an Automatic Computing Machine", Proc. Cambridge Phil. Soc, Vol. 47, pp. 96-108 (1951). (39) Carnahan, B., Luther, H. A. and Wilkes, J. 0. Applied Numerical Analysis , John Wiley and Sons, Inc., New York, 1969. (40) B6700 System Reference Manual. Burroughs Corporation, Detroit, Michigan. The author is indebted to Messrs. D. Mclntyre, T. Layman, and J. C. McMillan for this information. (41) Zabusky, N. J. "Solitons and Bound States of the Time Dependent Schrodinger Equation", Phys. Rev., Vol. 168, No. 1, pp. 124-128 (1968). (42) Wadati, M. and Toda, M. "The Exact N-Soliton Solution of the KdV Equation", Jour. Phys. Soc. of Japan, Vol. 32, No. 5, pp. 1403-1411 (1972). 166 (43) Saito, N. and Ooyama, H. "On the Stability of Lattice Solitons", Prog. Theor. Phys. Suppl., Vol. 45, pp. 201-208 (1970). (44) Helleman, R. H. G. "Heat Transport in Harmonic and Anharmonic Lattices", dissertation, University of Rochester (1972 ) (unpublished) . (45) Lebowitz, J. L. and Frisch, H. L. "A Model of Nonequilibrium Ensemble: Knudsen Gas", Phys. Rev., Vol. 107, pp. 917-930 (1957). (46) Van der Pol, B. and Weijers, T. J. "Tchehycheff Polynomials and Their Relation to Circular Functions, Bessel Functions and Lissajous- Figures", Physica, 's-Grav., Vol. 1, pp. 78-96 (1933), also in Selected Scientific Papers , Vol. II, North-Holland Publishing Company, Amsterdam, Netherlands, 1960. (47) Byrd, P. F. and Friedman, M. Handbook of Elliptic Integrals for Engineers and Physicists , Springer-Verlag, Berlin, West Germany, 1954. 167 VITA. The author, Ken'ichi Miura was born in Futatsui-machi, Akita-ken, Japan on 25 July, 1945. He received his B.S. in Physics from the Univer- sity of Tokyo in March, 1968. Upon graduation, he entered the Graduate Collage in the Faculty of Science, University of Tokyo. In September 1968 he began his graduate work in the Department of Computer Science in the University of Illinois at Urbana-Champaign. He has been a research assistant since then, first in the Illiac-IV project and then in the Center for Advanced Computation, working for the Application Group of the ILLIAC-IV. He received his M.S. in Computer Science in October 1971. He is a member of the Phi Kappa Phi Society and an associate member of the Sigma Xi Society. jibliographic data ;heet 1. Report No. UIUCDCS-R-73-586 Title and Subtitle The Energy Transport Properties of One Dimensional Anharmonic Lattices 3. Recipient's Accession No. 5. Report Date July 20, 1973 Author(s) Ken 'ichi Miura 8. Performing Organization Rept. No - UIUCDCS-R-73-586 Performing Organization Name and Address Department of Computer Science University of Illinois at Urbana-Champaign Urbana, Illinois 61801 10. Project/Task/Work Unit No. 11. Contract /Grant No. 2. Sponsoring Organization Name and Address 13. Type of Report & Period Covered Ph.D. Dissertation 14. 5. Supplementary Notes 6. Abstracts The energy transport properties of the solitary waves are studied both analytically and numerically. It was found that the piecewise linear function model clearly explains the characteristics of such waves. Numerical experiments were performed extensively to clarify the various properties of the solitary waves in the one dimensional anharmonic lattice. These results were used to interpret the heat conduction problem of the crystal lattice . 7. Key Words and Document Analysis. 17o. Descriptors Anharmonic Lattice Solitary Wave Thermal Conductivity One Step Methods Runge-Kutta-Ny strom Formula 'b. Ideniif iers/Opcn-Ended Terms 'c. ( OSAT1 Fie Id /(.roup 9. Availability Statement Release Unlimited 19. Security Class (This Report) UNCLASSIFIED 20. Security Class (This Page UNCLASSIFIED 21. No. of Pages 172 22. Price ORM N TI3--JS ( 10-70) USCOMM-DC 40329-P7 1 s s pp> JVJV. 5\tf*