Procedural approach to classical mechanics

Table of Contents

Note. For more details see Structure and Interpretation of Classical Mechanics (SICM) by G. J. Sussman and J. Wisdom (2001). The book revisits the fundamentals of classical mechanics, integrating theoretical concepts with formalized computational procedures.

1 Traditional notation

Physics notation is often ambiguous. Unstated notational conventions are effective tools if the context is clear. However, a formal definition of procedures about how to automatically compute (as by a computer) a given quantity requires unequivocal syntax.

As an example, let's consider Euler-Lagrange equations. In the traditional notation they read:

\begin{equation} \frac{d}{dt} \frac{\partial L}{\partial \dot{q}^i} - \frac{\partial L}{\partial q^i} = 0 \;. \end{equation}

The use of ambiguous notation (the Lagrangian \(L\) must be interpreted as a function of position \(q^i\) and velocity \(\dot{q}^i\), but for the time derivative to make sense, solution paths must be inserted in the partial derivatives) is convenient in this simple case. Still, implicit and equivocal mathematical notation is an obstacle to interpret the equation above programmatically.

2 Functional notation

M. Spivak in Calculus on Manifolds argues in favor of a functional notation. For instance, the chain rule would read \(D(f \circ g) = ((Df) \circ g) \cdot Dg\), so at \(a\):

\begin{equation} D(f \circ g)(a) = Df(g(a)) \cdot Dg(a) \end{equation}

Given the proper basic symbols definitions, this expression is unequivocal and it does not require the introduction of irrelevant letters (as is the case of the traditional notation).

The main advantage of functional notation over traditional notation is that it is well adapted to be interpreted automatically, independently from the context. Programming languages cannot deal in a robust, general way with apparent contradicting statements that a human would know how to interpret given the context. (There are exceptions, though. See for instance the prototypical language discussed in G. J. Sussman and A. Radul, The Art of the Propagator, 2009.)

2.1 Euler-Lagrange equations

Using a functional notation, Euler-Lagrange equations read as follows:

\begin{equation} D(\partial_2L \circ \Gamma[q]) - \partial_1L \circ \Gamma[q] = 0 \;. \end{equation}

The Lagrangian is a function of the tuple composed by the time, coordinates and velocities \(L(t, x, v)\). The partial derivative suffix indicates the argument position in the \((t, x, v)\) tuple. The function \(\Gamma\) maps time to the local tuple \(\Gamma[q](t) = (t, q(t), Dq(t))\), given the coordinate path \(q\). This notation can be easily written in terms of a robust procedural syntax:

(define ((Lagrange-equations Lagrangian) q)
  (- (D (compose ((partial 2) Lagrangian) (Gamma q)))
     (compose ((partial 1) Lagrangian) (Gamma q))))

This defines a procedure called Lagrange-equations that takes a Lagrangian and a coordinate path q as inputs to compute the Euler-Lagrange equations. Without going into the details of the syntax definition, let's stress that the first argument of each parenthesis names a procedure. For instance, (+ 1 2) indicates the sum of the integers 1 and 2, which evaluates to 3. Internal lists are evaluated first: given the expression (+ 1 (* 2 3)), first the product between 2 and 3 is evaluated to 6, which is then summed to 1.

Note the compose procedure to compose the Lagrangian procedure with the one resulting from evaluating (Gamma q). As mentioned above, given a coordinate path (here up indicates an up-tuple, think about it as a vector):

(define q
  (up (literal-function 'x)
      (literal-function 'y)
      (literal-function 'z)))

the procedure Gamma returns the local tuple composed by the time, the coordinate tuple and the velocity tuple:

(print-expression ((Gamma q) 't))

 (up (x t) (y t) (z t))
 (up ((D x) t) ((D y) t) ((D z) t)))

(The result of evaluating the procedure is shown enclosed by #|...|#.)

The procedure defined for the Euler-Lagrange equation can be applied to any Lagrangian. Let's consider a couple of examples.

2.1.1 The free particle

Let's apply the procedure defined above to the Lagrangian \(L=\frac{1}{2} m v^2\) for a free particle of mass \(m\). The corresponding procedure is:

(define ((L-free-particle mass) local)
  (let ((v (velocity local)))
    (* 1/2 mass (dot-product v v))))

Let's compose the Euler-Lagrange equations procedure with the free particle Lagrangian

(show-expression (((Lagrange-equations (L-free-particle 'm))
 (literal-function 'x)) 't))

(* m (((expt D 2) x) t))

This provides the expected equation of motion for a free particle \(mD^2x(t)=0\).

2.1.2 The harmonic oscillator

The Lagrangian \(L(t, q, v)=\frac{1}{2}mv^2 - \frac{1}{2}kq^2\) of a particle of mass \(m\) constrained at the end of a spring (characterized by a spring constant \(k\)) is described by the following procedure :

(define ((L-harmonic m k) local)
  (let ((q (coordinate local))
        (v (velocity local)))
    (- (* 1/2 m (square v)) (* 1/2 k (square q)))))

Let's check that the motion of a harmonic oscillator \(x(t)=a cos(\omega t)+\phi\):

(define (proposed-solution t)
  (* 'a (cos (+ (* 'omega t) 'phi))))

satisfies the Euler-Lagrange equations given the Lagrangian above:

 (((Lagrange-equations (L-harmonic 'm 'k))

This gives

\begin{equation} a \cos\left( \omega t + \phi \right) (k - m {\omega}^{2}) \end{equation}

As expected, the motion is harmonic and the solution allowed are those satisfying \(\omega^2=k/m\).

3 About the procedural notation

The procedures formalized in the examples above are robust, i.e., not bound to a specific context. The Euler-Lagrange equations procedure generalizes automatically to arbitrary Lagrangians or coordinate paths. It can be composed with other arbitrary functions, too (as long as their domains overlap).

The procedural syntax used in the examples is the Scheme dialect of the Lisp family of programming languages. Lisp was designed by John McCarthy in 1960 (which makes it one of the older programming languages still in use, second only to Fortran). The name stands for List Processing, because one of the key ideas is to use a simple data structure called a list for both code and data. The Roots of Lisp by P. Graham (2001) revisits McCarthy's original article illustrating the code/data evaluation procedure.

Here we skipped several basic definitions and only aimed to give a flavor about the Scheme notation used in SICM. For more details, refer to:

  • SICM, appendix Scheme: an introduction to the Scheme programming language.
  • SICM, appendix Our notation: details about the functional mathematical notation adopted by the book and how generic arithmetic extends basic operations, such as addition and multiplication, to a variety of mathematical types. The discussion about differentiation and of tuple (a generalization of tensors) arithmetic is also interesting.

4 Appendix: Procedures composition

Note. This section is rather technical and aimed to illustrate a subtlety used above that shows the flexibility of the Lisp syntax (despite its simplicity).

The procedure L-free-particle was defined as a function of mass and of the local tuple as:

(define ((L-free-particle mass) local)
  (let ((v (velocity local)))
    (* 1/2 mass (dot-product v v))))

It may seem simpler to define it as (define (L-free-particle mass local) ...). However, we want to compose (L-free-particle mass) with other functions by itself.

Let's consider a minimal example. The following definitions are equivalent and bound the procedure (lambda (x y) (+ x y)) (where lambda denotes an anonymous function) to the symbol foo:

(define foo (lambda (x y) (+ x y)))

(define (foo x y) (+ x y)) ; Same as before, using syntactic sugar.

Let's consider also the following equivalent definitions:

(define bar (lambda (x) (lambda (y) (+ x y))))

(define (bar x) (lambda (y) (+ x y)))   ; Syntactic sugar.

(define ((bar x) y) (+ x y))            ; More syntactic sugar.

For the simple purpose of evaluating \(x+y\) (e.g., for \(x=2\) and \(y=3\)), both cases work:

(foo 2 3)
;Value: 5

((bar 2) 3)
;Value: 5

However, the second case allows to compose bar with other functions of one argument. For instance, let's compute \(x^3+y\) (for \(x=2\) and \(y=3\)) as:

(((compose bar cube) 2) 3)
;Value: 11

A similar composition of foo with cube is not possible because foo requires two arguments, while cube one.

It is remarkable how a simple idea such as using lists for both code and data naturally leads to (a priori non-trivial) procedure composition.

5 Appendix: scmutils

The examples shown above are based on the powerful scmutils library for the Scheme programming language. The library is libre software, developed for education and research at MIT.

Once installed, all the procedures illustrated above can be easily evaluated through the interactive Scheme REPL (read-eval-print loop), similar to (but way more powerful than) an interactive shell.

5.1 Installation

Official installation instructions are provided here. However, a local installation may be preferable:

  1. Download the system from the official website.
  2. Extract the tar ball in a local folder (e.g., ~/scmutils):
    tar -xvzf [your-tarball.tar.gz]

    This will extract two folders: bin and scmutils.

  3. Now you have two options: access the Scheme REPL with Edwin (an Emacs-like editor written in Scheme) or Emacs.
    1. Edwin. Set the SCMUTILS_ROOT environment variable to override the default location (/usr/local) for the scmutils root directory:
      export SCMUTILS_ROOT

      Launch bin/mechanics script.

    2. Emacs. Copy the following mechanics Elisp command (e.g., in ~/.emacs or in the *scratch* buffer) and evaluate it (change the root path below to the full one corresponding to the scmutils folder extracted from the tar ball):
      (defun mechanics ()
        (let ((root "/path/to/scmutils/"))
           (concat root "/mit-scheme/bin/scheme --library " root "/mit-scheme/lib"))))

      Launch the command as M-x mechanics from any Emacs buffer. This will open the Scheme REPL with the scmutils library already available.

  4. If everything seems to work you can now remove the tarball.

5.2 Usage

After having launched the Scheme REPL in Ediwn/Emacs, open a new buffer C-x 2. Open a source code file C-x C-f demo.scm. After adding some code, evaluate and send sexps to the REPL by typing C-x C-e at the end of the last parenthesis, or C-M-x within a given sexp. (This assumes that only default modes are installed in Emacs. If, for instance, Geiser mode is also installed, see its usage instructions or toggle geiser-mode to disable it.)

To quit the REPL, type (exit).

Useful documents are included in the scmutils/manual/ folder. Scheme sources for the entire system are included in the scmutils/src/ folder.

Tip: when editing Scheme source code in Emacs, install the useful paredit extension.

Back to blog

Date: <2018-02-28 Fri>

Author: Francesco Montanari

Created: 2018-02-28 Wed 23:18

Emacs 24.5.1 (Org mode 8.2.10)