% (fh)--------------------------------------------=  ii.tex        =--
%
%   This file was
%   written  in TU Budapest QTG by= Varga Imre       at= 10-05-89 11:14
%
%   modified in TU Budapest QTG by= M rk G‚za        at= 11-05-89 21:19
%
% (fh)----------------------------------------------------------------
%
%       P l Pacher, G‚za M rk, L szl¢ Udvardi and Imre Varga:
%       Teaching Physics by Means of Computer Modeling,
%       Eighth Summer School on Computing Techniques in Physics,
%       Skalsky dvur, Czechoslovakia, September 1989
%       Proceedings editor: J. Nadrchal, Institute of Physics Praha
%
\magnification=1200
%\nopagenumber
\hsize=15truecm
\vsize=22.0truecm
\baselineskip=0.90truecm plus 0.1truecm
\parindent=1truecm
\centerline {\bf TEACHING PHYSICS BY MEANS OF COMPUTER MODELING}
\vskip 0.2truein
\centerline {P\'al Pacher$^{\dagger }$,
    G\'eza M\'ark$^{\ast\ddagger}$,
    L\'aszl\'o Udvardi$^{\star }$ and Imre Varga$^{\star }$}
\vskip 0.2truein
\centerline {\it $^{\dagger }$Department of Physics, $^{\ast }$
     Department of Atomic Physics,}
\centerline {\it $^{\star }$Quantum Theory Group,}
\centerline {\it Institute of Physics, Technical University, 1521
    Budapest, Hungary.}
\vskip 0.4truein
\noindent {\bf Abstract}
\vskip 0.3truein
    Computer modeling in teaching physics has a growing importance. It
    helps to solve and understand complex physical problems using
    numerical methods, computer simulation and animation. Typical
    examples developed in our institute are presented: the 3--body
    problem, the Van der Pol oscillator, the recursive graphs and the
    problem of the ideal gas simulation. Both the numerical and the
    modeling difficulties are discussed.
\vskip 3truecm
\noindent $^{\ddagger }$proofs should be sent to this author.
\vfill
\eject
\vskip 0.4truein
\noindent {\bf 1. Introduction}
\vskip 0.3truein
    Traditionally physics is divided into two main branches:
    experimental and theoretical physics.

    The experimental physicist tries to understand nature by
    performing experiments  on some representative samples. He records
    sample properties which can be precisely measured against
    variables under his control. In order to understand his results
    and to generalize them to other samples he needs the help of a
    theoretician.

    The theoretical physicist is not interested in the irrelevant
    details of the experimentalist's samples. By assuming a theory for
    the behavior of real materials, he works out those same properties
    the experimentalists had measured as exactly as possible.

    But some properties of the sample might be easy to measure, but
    too difficult for the theoretical physicist to evaluate without
    unreasonable approximations. In other cases theories predict
    interesting consequences that are impossible, or at least very
    hard, to measure.

    In both cases a product of our age, {\it computer modeling}, can
    help in solving the problem. The computational physicist uses a
    finite size model, which is a compromise between the accuracy of
    his model and the size of available computer storage and computer
    time. He runs the model and makes ''measurements'' on it.
    These results can be compared with the experiments (if there are
    any), or they themselves are data to compare with theoretical
    predictions. In contemporary physics a number of subjects rely
    more and more on computational physics.

    Computer modeling is, of course, not less important in teaching
    physics. Computer models operating within the framework of
    classical physics can show the students a close-up picture of
    materials. They can reveal why we see various sample-averaged
    properties measured by experiments or calculated by theory.
    Maxwell invented a ''demon'' who could watch the atoms going by.
    Nowadays students can observe the motion of particles on the
    screen of a PC and can open or close the shutter to
    separate atoms with different velocities or other properties and
    follow the entropy change of the process.

    At the Faculty of Electrical Engineering of the
    Technical University
    Budapest a new branch, Computer Engineering, has been brought to
    light two years ago. In the third semester the students have to
    write a complex program for practicing the computer languages they
    have already learned, such as Pascal and C. Many students solve
    different problems in physics, most of them by computer modeling.
    Besides developing skills, it also helps them to understand how
    nature works. In what follows we describe some of these programs
    written for IBM compatible PC-s.
\vskip 0.4truein
\noindent {\bf 2. Moons of Saturn}
\vskip 0.3truein
    The simulation of three--body problems is always a big challenge.
    This time the task was the simulation of the motion of the
    co-orbital moons of Saturn discovered by Voyager 1 [1]. These
    moons of approximately the same size and mass follow the same
    orbit with slightly different distance from Saturn. The inner moon
    revolves obviously a bit faster than the outer one. Due to the
    gravitational attraction between the two moons, however, they
    never collide but exchange their orbits when they are close
    enough. The aim of the simulation was to show that the minimum
    distance of the two moons is nonzero.

    For the integration of the Newtonian equations, we have used the
    energy conserving algorithm by Greenspan [2]. The algorithm
    involves the solution of the Newtonian equations of motion for the
    interaction
    $$F=-{{G}\over {r}}+{{H}\over {r^n}}-\alpha v,$$
    where the dominant part is the attraction and the repulsive and
    the friction terms are present in order to balance the instability
    problems arising for taking nonzero timesteps. We have chosen
    $G=2250$, $H=1$, $n=3$, $\alpha =10^{-6}$ taken from [2]. We have
    found that Greenspan's algorithm was faster than usual integration
    procedures but the accuracy was considerably lower. The fine
    tuning of the parameters appearing in the force formula may,
    however, improve the results.

    As it is shown on Figure 1. the simulation has yielded the
    expected answer: the two moons under consideration never collide
    but approach each other to about $21.5\% $ of their largest
    separation. If we chose $H=0$ we may obtain
    approximatelly the same answer.
\vskip 0.4truein
\noindent {\bf 3. Van der Pol oscillator}
\vskip 0.3truein
    The solution of several physical problems, like the three--body
    problem in classical physics or the problem of
    turbulence, exhibits chaotic properties. One of the historically
    most important  example is the Van der Pol oscillator. It was the
    first system in which the existence of chaotic solutions could be
    proved [3]. The study of the Van der Pol oscillator had
    considerable contribution to the evolution of the theory of
    dynamical systems [4].

    We made the oscillator be chaotic using a tunnel diode having
    nonlinear U--I characteristic. See Figure 2.
    A set of coupled ordinary differential equations may
    be derived describing this oscillator [5], where a nonlinear
    function represents the characteristics of a tunnel diode
    appearing in the circuit.
    $$L{{dI}\over {dt}} = {{MS-RC}\over {C}}I + U - V,$$
    $${{dU}\over {dt}} = -{I\over C},$$
    $$C_1{{dV}\over {dt}} = I - f(V),$$
    where function $f(V)$ was chosen as
    $$f(V) =\cases {\alpha (e^{\gamma V}-1), &if $V \leq 0,$ \cr \cr
                         aV^3+bV^2+cV, &if $0 < V \leq V_0,$ \cr \cr
                         e^{\beta (V-V_0)}+d, &if $V > V_0$,}$$
    where the parameters were chosen requiring the existence of the
    first derivative of the function $f(V)$: $\alpha =0.0945A$,
    $\gamma =100V^{-1}$, $a=56.688AV^{-3}$, $b=-53.072AV^{-2}$,
    $c=9.45AV^{-1}$, $\beta = 2.7945V^{-1}$, $V_0=0.787V$, and
    $d = 0.78769A$. In order to solve the set of differential
    equations
    a third order Runge--Kutta procedure using analytical derivatives
    was applied. The waveform of the solutions and the strange
    attractor of the system were determined changing the parameters of
    the oscillator. As one can see on Figure 3. the attractor of the
    system tends to a limitcycle for the following set of parameters
    $$g = {{U}\over {I_m}}\sqrt{{{C}\over{L}}} = 0.05, \qquad
    {{C}\over {C_1}}g = 0.02, \qquad {{MS-RC}\over {\sqrt{LC}}} =
    0.1.$$
\vfill \eject
\vskip 0.4truein
\noindent {\bf 4. Fractal dimension}
\vskip 0.3truein
    Computer modeling is essential in teaching the behavior of
    nonlinear systems. One of the most interesting features
    in such systems is the appearance of fractals, which has also been
    discovered in many other fields of life. Fractals happen to be one
    of the most complicated and yet one of the most beautiful
    mathematical objects as anyone can verify it from the wonderful
    pictures in Mandelbrot's [6] famous book.

    The fractal dimension describes the volumetric structure of any
    set of points distributed in real space. It may be used to
    characterize special graphs, strange attractors, i.e. fractals in
    general.

    In this program two dimensional recursive graphs have been generated:
    the Sierpinsky carpet, the Sierpinsky graph, the
    Hamilton graph, and the Bethe lattice (see Figure 4.). The fractal
    dimension of these graphs have been calculated by means of the
    basic definition introduced by Hausdorff and Kolmogorov [6]
    $$d=\lim_{n\rightarrow\infty}{{\log (P(n))}\over {-\log(\epsilon
    (n))}},$$
    where $n$ is the fractal index, $P(n)$ is the number of squares
    necessary to cover the $n$-th fractal and $\epsilon (n)=0.5^n$
    is the length of the sides of these squares. Convergence based on
    the definition was fairly slow and the generation of the recursive
    graphs is also a computer consuming task. We have found
    $d=\log (3)/\log (2)$ for the Sierpinski carpet, $d=2$ for the
    Sierpinski graph and the Hamilton graph, and $d=1$ for the Bethe
    lattice.
\vfill \eject
\vskip 0.4truein
\noindent {\bf 5. Ideal Gas Simulation [7]}
\vskip 0.3truein
         Statistical physics is  a profitable field  for computer
    modeling.
    Here  one always  studies a  system consisting of
    many  particles.
    In  a  computer  model  the  microscopic -
    macroscopic metamorphosis is witnessed, i.e., the multi-faced
    behavior of the system is  built up by the motion  of several
    particles governed by simple laws.

         The   statistical   physical   simulation   methods  are
    classified  as  {\it  Monte-Carlo} [8]  or {\it molecular
    dynamics } [9,10,11,12,14] techniques. The    Monte-Carlo
    calculations  are  usually  faster  but  the  moves  of   the
    particles  are  artificial  rather  than dynamical.  For this
    reason only the equilibrium properties can be calculated.

         In this section we  present an  ideal gas  model based on
    molecular dynamics.  The modelled objects are mass points and
    a  container  made  up  from  various  types  of  walls.  The
    particles  playing   the  role   of  gas   molecules  collide
    elastically  with  each  other.   The  different interactions
    modellized  by  the  different  walls  are represented by the
    corresponding rules for collision against the walls.

         This program
    is useful as a visual aid and experimental tool
    in  various  levels  of  education  including  the  following
    topics:
    molecular motion, temperature, first and
    second laws of thermodynamics, speed distribution, fluctuations [13],
    molecule formation, relaxation phenomena.

         The program displays on the screen the moving particles
         confined to a rectangular vessel and
    three real time diagrams
    (G1, G2, G3) simultaneously. See Figure 5.
         There are  three types  of particles: red, green
    and invisible ones. The invisible particles are useful when
    you want to concentrate just to certain molecules. (E.g.
    in simulation of gas mixing, see later). Just make those
    molecules visible and the others invisible!
    The particle number and the mass of each particle type can
    be chosen
    independently.  The maximal number of the particles is at
    most 999.

        The types of the walls are:

        - adiabatic: i.e. the particles collide elastically against
    it,

        - diathermic: it can exchange energy both with the
    particles and with the heat reservoir.

        The walls can move along the normal or the tangential
    direction to their surface.
    Moreover, a {\it separation wall} can be placed into the
    container. There is a {\it slit} on this wall, the user can
    open or close this slit or put a {\it Maxwell's demon} into the
    slit.

         The  on-line  diagrams  can  show  a  variety of distribution
    functions (e.g. pressure, density, temperature versus position,
    histograms of velocity and energy),  time-evolution
    functions (e.g.
    entropy).  Time averaged  distribution functions
    may be displayed, as well.
         Other possible useful features are the {\it shot noise}
    (a clack in every wall-particle collision) and
    {\it particle tracing}
    which shows the path of one particular mass point.
\vskip 0.2truein
\noindent{\it 5.1. Free Expansion of the Gas Into Vacuum}
\vskip 0.1truein
         The particles start from the upper left corner into
    different directions with uniform speed distribution.
         The container is quickly filled up by the particles.
        The velocities will be changed by the collisions and one
    arrives to the Maxwell velocity distribution. Its form in two
    dimensions is given by: $ F(v) = A v e ^ { - v^2 / v_0^2 } $.
        If the collisions of particles with each other are
    switched off
    the velocity
    distribution is kept in its initial form in case of
    adiabatic boundaries.
        When at least one
    boundary is diathermic, however, the energy distribution of
    the gas converges to the energy distribution of the
    reservoir.
\vskip 0.2truein
\noindent{\it 5.2. Mixing of Two Types of Particles}
\vskip 0.1truein
         The green particles start from the upper left corner, the red
    ones start from the opposite  corner.  The initial
    velocity and the mass  of both types may be chosen
    independently. A separation wall with a hole may be
    placed into the container with selectable location and hole size.
         By displaying the velocity distribution in G1 and the
    time-temperature
    functions of the two types of particles in G2  and G3
    the thermalization of the system [16]
    may be observed, i.e. the
    average energy for both types of particles will be the same.
    If the
    greens' mass  is greater  than the  reds' mass,  the greens' average
    velocity will be smaller yielding a two peaked
    velocity distribution (see Figure 5).

         The  role of a Maxwell's demon may be played when
    opening or closing the hole on the wall, i.e. the particles may be
    separated to the left and right side of  the container according
    to some attribute (e.q. color, speed).
\vskip 0.2truein
\noindent{\it 5.3. Boltzmann Distribution}
\vskip 0.1truein
         Let's switch the gravitational field on!

         If collisions are allowed only against the walls the particles
    move  along  independent  parabolic  paths.   If they start
    with the same speed each particle reaches the same maximal
    height, above that height one gets vacuum.  This strange picture
    is dramatically changed when mutual collisions are permitted.
    Energetic particles reach  high elevations while slow particles
    bounce just on the floor.  The particle density versus height
    may be checked and it proves to be similar to the theoretical
    exponential curve. When for the two types of particles different
    masses are chosen, the lighter ones reach higher elevations due to
    their larger average velocity.
\vskip 0.2truein
\noindent{\it 5.4. Maxwell's Demon}
\vskip 0.1truein
         The demon is sitting in the hole of the separation wall
    dividing the container into two parts.  She lets pass
    particles coming from the left through the hole only if their
    energy is greater than  a  preselected  threshold.
    Particles coming  from the  right are  let through  only with
    energy smaller  than the  threshold.
         When the threshold energy is zero the demon
    operates as a pump. The particles starting from  the upper left
    corner first uniformly fill both sides of the vessel then the
    demon slowly pumps them to the right side.  This strange process
    is shown in Figure 6.
\vskip 0.2truein
\noindent{\it 5.5. Motion of a Piston}
\vskip 0.1truein
         The right wall is pushed into the container then  pulled
    out  again.   The  compression  rate  and the velocity of the
    piston may be selected. Particles colliding against the moving
    piston change their energy leading to heating or cooling  the
    gas. The process is adiabatic, because there is no heat transfer,
    just mechanical work is done.

         How much work is done by a given compression?  That  depends
    on the velocity of the piston as it may be demonstrated by
    plotting the average temperature versus  time.  Maximal work  is
    done by  slow, quasistatic piston motion -  one arrives to the
    law [15] $pV^{\kappa} = const$. The work is almost zero  if you
    pull out quickly the piston, because just only a  few  particles
    hit the piston during its motion. In this  case the process  is
    isothermic, i.e. $pV  = const.$
\vfill
\eject
\noindent {\bf References}
\vskip .1truein
\parindent -0.20truecm
\baselineskip=0.80truecm plus 0.1truecm
    1\ R. Gore, National Geographic {\bf 160} (1984) 3

    2\ D. Greenspan, SIAM J. Appl. Math. {\bf 20} (1971) 67

    3\ M. L. Cartwright, J. E. Littlewood and J. London, Mat. Soc.
    {\bf 20} (1945) 180

    4\ N. Levinson, Ann. Math. {\bf 50} (1949) 127; S. Smale, Bull.
    Am. Math. Soc. {\bf 73} (1967) 747; G. Guckenheimer, Physica
    {\bf 1D} (1980) 227

    5\ A. S. Pikovsky and M. I. Robinovich, Physica {\bf 2D} (1981) 8

    6\ B. Mandelbrot, {\it Fractals: Form, Chance and Dimension} (W.
    H. Freedman, San Francisco, 1977)

    7\ E. H. Kennard,
       {\it Kinetic Theory of Gases
       With an Introduction to Statistical Mechanics
       }
       (McGraw-Hill, New York, 1938)

    8\ J. Novak and A. B. Bortz,
       Am. J. Phys. {\bf 38} (1970) 1402

    9\ P. Empedocles,
       J. Chem. Educ. {\bf 51} (1974) 593

   10\ B. J. Adler and T. E. Wainwright,
       J. Chem. Phys. {\bf 31} (1959) 459

   11\ {\it Kinetic Theory by Computer Animation}, a film
       produced by J. T. Fitch, J. L. Kinsey and S. F. Martin.
       (Kaima Co., Dept. P 2. Concord, MA 01742.)

\noindent
       Reviewed in Am. J. Phys. {\bf 44} (1976) 810

   12\ E. T. Lane, {\it Simulated Waves and Paricles}, a program for
       APPLE II.C. (CONDUIT RM 4557, Oakdale Hall,
       The University of Iowa, Iowa City 52242)

\noindent
       This program, however, don't model the mutual collisions of
       particles.

   13\  J. R. Ray,
        Am. J. Phys. {\bf 50} (1982) 1035

   14\  T. Tajima, A. Clark, G. G. Craddock, D. L. Gilden,
        W. K. Leung, Y. M. Li, J. A. Robertson and
        B. J. Saltzman,
        Am. J. Phys. {\bf 53} (1985) 365

   15\  M. I. Sobel,
        Am. J. Phys. {\bf 48} (1980) 877

   16\  J. Berger,
        Am. J. Phys. {\bf 56} (1988) 923
\vfill
\eject
\noindent {\bf Figure Captions}
\vskip .1truein
\parindent 0truecm
    Figure 1. Distance of the moons versus time. $R$ is the average
    distance of the moons measured from Saturn. $T$ is the average
    time of the orbits of the moons around Saturn. Model parameters
    were used [2]. (Mass of Saturn $M=10$, masses of the moons $m=0.01$,
    $dt=0.01$, initial distance of the moons from Saturn $R_1=101$
    and $R=99$, $T\approx 132.5$.)
\vskip .2truein
    Figure 2. Scheme of the oscillator and the V--I model characteristics
of the tunnel diode.
\vskip .2truein
    Figure 3. Attractor of the Van der Pol oscillator.
    $X = {{I}\over {I_m}}$, $Y = {{U}\over
    {I_m}}\sqrt{{{C}\over{L}}}$,$Z = {{V}\over {V_m}}$ \hfill \break
\vskip .2truein
    Figure 4. Recursive graphs: {\it a.} Sierpinski carpet, {\it b.}
    Sierpinski graph, {\it c.} Hilbert graph, {\it d.} Bethe lattice.
\vskip .2truein
    Figure 5. Demonstration of thermal equilibrium and
    relaxation. The vessel contains 150 red and 75 green
    particles (colors not shown in this figure).
         Initially they had different temperatures.
         The mass of the green particles is 36 times
    greater than the reds' mass yielding a two peaked velocity
    distribution shown on G1.
         Curves G2 and G3 show the average temperature of the
    green and red particles versus time, respectively.
\vskip .2truein
    Figure 6. Maxwell's demon in action.
         The demon's threshold energy is zero.
         G1 shows the entropy of the system, G2 shows the horizontal
    density distribution.
\bye
