Exam 2000

Instructions

This exam, worth 40% of the total mark for Computational Physics, is due in by 5 pm on Monday May 22. Each student is to do the exam by themselves and submit their solutions individually.

The exam consists of 6 questions of varying difficulty with numerous hints as to how to proceed. You should attempt to answer whole questions whose total value is at least 110 marks. The marks for each question is clearly indicated in red brackets (i.e., [n]) and, obviously, the marks provide a good guide to the question's difficulty. Bonus marks will be given for elegant or unconventional solutions.

I quote from the Departmental Guidelines on Plagiarism:

Like the scientific process, assessment of students relies on the integrity of the participants. The Department expects that any piece of work submitted by a student for assessment will be essentially their own work, and that contributions of others to that work will be appropriately acknowledged. In essence, students are expected to write their own essays and assignments, just as they are expected to sit the exams themselves.

Students should generally avoid verbatim copying of published work, even when the work is cited. They are strongly encouraged to express their knowledge and thoughts in their own words, rather than those of others, as this greatly assists the learning process. It is acknowledged however that in some circumstances it is desirable to quote directly from a published work.

Submitting your Solution

Please name your exam solution using your initials, e.g., PCA.nb, and submit to ftp://ftp.physics.uwa.edu.au/incoming/Exam.

. Thermodynamics [40]

.. Equation of State [1]

Equations of state (EOS) are widely used — both in textbooks and in industry — for representing the pressure (P), volume (V), and temperature (T) relationships of real fluids. Many EOS are explicit in P and cubic in V:

P = ((R T)/(V - b)) -  
        (a/(w (b^2) + u V b + (V^2))),

where R is the gas constant, and a and b are parameters to be determined. In some EOS the parameter a is a function of T, i.e., a(T). The quantities u and v are usually given integer values.

.. Polynomial Form [3]

The general EOS can be expressed as a cubic polynomial in V.

.. Special Cases [2]

.. Compressibility [2]

An important thermodynamical quantity is the compressibility, defined by z = ((P V)/(R T)).

.. Critical Point [5]

The parameters a and b (or A and B) may be evaluated by making use of the fact that the P-V relationship has a point of inflection at the critical point. That is:

((∂P)/(∂V))=0 and (((∂^2) P)/(∂(V^2)))=0

when P = (Subscript[P, c]), T = (Subscript[T, c]), V = (Subscript[V, c]), and (Subscript[z, c]) = (((Subscript[P, c]) (Subscript[V, c]))/(R (Subscript[T, c]))), where (Subscript[P, c]), (Subscript[T, c]), (Subscript[V, c]), and (Subscript[z, c]) are the critical pressure, temperature, volume, and compressibility respectively. These properties have been determined experimentally for a large number of fluids and have been estimated for many others.

Equation (.) is equivalent to requiring that the volume or compressibility polynomials have three equal real roots:

(((V - (Subscript[V, c])))^3)=0=(((z - (Subscript[z, c])))^3).

.. Redlich-Kwong EOS [7]

The Redlich-Kwong (RK) EOS reads

P = ((R T)/(V - b)) - (a/((T^(1/2)) 
                V (V + b))).

..  Maxwell Construction [20]

As you will see, below the critical temperature, t = 1, loops appear in the P-V diagram. These loops are forbidden for thermodynamical reasons. For t < 1, a phase transition between gas and liquid takes place; at t = 1 the difference between the liquid and gaseous phases vanishes.

Over a certain temperature range, three volumes (Subscript[v, 1]) ≤ (Subscript[v, 2]) ≤ (Subscript[v, 3]) correspond to each pressure p. For thermodynamic reasons, the transition from the gas with large volume (Subscript[v, 3]), to the liquid, with small volume (Subscript[v, 1]), takes place at the pressure \(TraditionalForm\`p\&_\) for which

\(TraditionalForm\`∫\_\(v\_1\)\%\(v\_3\)\(p(
            v)\) dv = \((v\_3 - 
            v\_1)\) \(\(p\&_\)\( \)\(.\)\)\)

At \(TraditionalForm\`p\&_\), gas and liquid are simultaneously present for all volumes v with (Subscript[v, 1]) ≤ v ≤ (Subscript[v, 3]); one observes a two-phase mixture — the coexistence region. Geometrically, equation (.) means that the area between the curve p(v) and the straight line \(TraditionalForm\`p\&_ = constant\) for (Subscript[v, 1]) ≤ v ≤ (Subscript[v, 2]) is the same as the corresponding area for (Subscript[v, 2]) ≤ v ≤ (Subscript[v, 3]). This is known as the Maxwell construction. We want to construct the line \(TraditionalForm\`p\&_ = constant\); the two-phase line. For the two-phase line we need two equations in order to determine the two unknown volumes (Subscript[v, 1]) and (Subscript[v, 3]):

({{p((Subscript[v, 1]), t) = p((Subscript[v, 3]), t),}, {(∫_((Subscript[v, 1]))^((Subscript[v, 3])))p(v, 
                    t) dv = ((Subscript[v, 3]) - (Subscript[v, 1])) p((Subscript[v, 3]), 
                    t).}})

There is no simple analytical solution to equations (.).

. Solitons [40]

The Korteweg de Vries (KdV) equation, describing water waves in shallow channels, is given by

((∂u)/(∂t))-6 u ((∂u)/(∂x))+(((∂^3) u)/(∂(x^3)))=0,

where u is the wave amplitude.

One-soliton Solution [9]

Two-soliton Solution [8]

Observe that the smaller slower soliton of the two-soliton solution appears identical to the one-soliton solution, but with a small phase shift.

Plot[{(Subscript[u, 1])(x, 2), (Subscript[u, 2])(x, 2)}, {x, 0, 15}, 
        PlotRange -> All];

From your visualization of the two-soliton solution, you should find that the velocity of the smaller slower soliton is

(Subscript[v, 1]) = +4;

In a frame of reference moving with this velocity (x -> x + (Subscript[v, 1]) t), the shape of the slow soliton as t -> ∞ is given as follows.

(Subscript[u, 2, slow])(x_, ∞) = 
      (lim_(t -> ∞))(Subscript[u, 2])(x + (Subscript[v, 1]) t , t)

Conservation Laws [10]

A continuity equation, relating a general density ρ and its corresponding current J, of the form

(∂_t) ρ + (∂_x) J = 
      0,

leads to a related conservation law. If both ρ and J depend on t, u, ∂u/∂x, (∂^2) u/∂(x^2), … but not on ∂u/∂t, and J vanishes as x -> ±∞, integrating equation (.) over \(TraditionalForm\`x ∈ \((\(-∞\), ∞)\)\), leads to

\(TraditionalForm\`∂\_t\(∫\_\(-∞\)\%∞ ρ dx\) = \(0 => ∫\_\(-∞\)\%∞ ρ dx = \(\(constant\)\(.\)\)\)\)

For example, the KdV equation (.) is already in the form of the continuity equation (.) with the density

ρ = u(x, t);

if the corresponding current J, is

J=∫((((∂^3) u(x, 
                      t))/(∂(x^3)))-6 u(x, t) ((∂u(x, t))/(∂x)))\[DifferentialD]x

because the continuity equation (.) yields the KdV equation (.).

(∂_t) ρ + (∂_x) J == 
      0

Hence, we conclude from the conservation equation (.) that (∫_(-∞)^∞) u(x, 
              t) dx = constant for solutions to the KdV equation (.) satisfying u(±∞, t) = 0. Physically, this conservation law expresses conservation of mass.

n-soliton solution [13]

There is an interesting relationship between the soliton solutions of the KdV equation (.) and the Schrödinger equation. Using this relationship, the n-soliton solution can be generated in an exact way.

Initial values (t = 0) of the one- and two-soliton solutions read.

(Subscript[u, 1])(x, 0)

(Subscript[u, 2])(x, 0) // Simplify

It turns out that the initial value of the n-soliton solution is (Subscript[u, n])(x, 
          0) = -n (n + 1) (sech^2)(x) = 
          2  (Subscript[𝒱, n])(x). If (Subscript[𝒱, n])(x) is interpreted as the potential of the Schrödinger equation, that is

\(TraditionalForm\`\(\(\(\(-1\)\/2\) ∂\^2\( ψ(x)\)\/∂x\^2 - \(1\/2\) \(n (n + 1)\) \(\(sech\^2\)(
            x)\) \(ψ(x)\) = ℰ\ \(ψ(
          x)\)\)\(,\)\)\)

the (unnormalized) solution is

(Subscript[ψ, n_, m_])(x_)=(Subscript[P, n]^m)(tanh(x));

where (Subscript[P, n]^m)(z) is the associated Legendre function, and the energy is (Subscript[ℰ, m]) = -(m^2)/2. For each n, there exists n bound states with 1 ≤ m ≤ n.

The next step in determining the n-soliton solution (much theory omitted here) is to form the n×n matrix, (Subscript[𝒜, n]), whose coefficients are

(Subscript[𝒜, n_, m_, k_]) := (Subscript[δ, m, k]) + ((Subscript[c, n, m]^2)/(m + 
                k)) (E^(8  (m^3) t - (m + k) x))

where

(Subscript[c, n_, m_]) := (Subscript[c, n, m]) = 
        (lim_(x -> ∞))Expand[(Subscript[𝒩, n, 
                  m]) (Subscript[ψ, n, m])(
                x) (E^(m x))]

and

(Subscript[δ, m_, k_]) := 
      KroneckerDelta[m, k]

Finally, the n-soliton solution is given by

(Subscript[u, n])(x, t)=-2 (((∂^2)log(\(|𝒜\_n|\)))/(∂(x^2))),

where \(TraditionalForm\`|𝒜\_n|\) is the determinant of (Subscript[𝒜, n]). For example, the two-soliton matrix is

(Subscript[𝒜, 2]) = 
      Table[(Subscript[𝒜, 2, m, k]), {m, 2}, {k, 2}]

and the n-soliton solution is

-2 (((∂^2)log(\[LeftBracketingBar](Subscript[𝒜, 2])\[RightBracketingBar]))/(∂(x^2)))

We compare this result with (Subscript[u, 2])(x, t) given above in question ..

Simplify[% == (Subscript[u, 2])(x, t)]

. Hydrogen Molecular Ion [15]

In this problem you will examine the bound states of a system comprising a single electron moving in the field of two fixed nuclei of arbitrary charge (commonly known as the generalized hydrogen molecular ion problem or reduced two-centre problem). This configuration is the simplest conceivable molecule, and yet is the most complex exactly solvable system.

An important special case is the hydrogen molecular ion, (Subscript[H, 2]^(+)), which is composed of two (fixed) protons and one electron. In atomic units ((Subscript[m, e]) = e=ℏ =(Subscript[a, 0]) =1), the Schrödinger equation for the motion of the electron reads

ℋ Φ(r)={(Underscript[-(1/2) (Subscript[Δ, r]), (Underscript[⏟, 𝒯])])+(Underscript[(-(1/(Subscript[r, A])) - (1/(Subscript[r, B])) + (1/R)), (Underscript[⏟, 𝒱])])}Φ(r)=ℰ Φ(r),

where R is the internuclear separation, and (Subscript[r, A]), (Subscript[r, B]), and r are the position vectors of the electron with respect to the protons A and B, and to the midpoint of the internuclear line, respectively.

The vectors (Subscript[r, A]), (Subscript[r, B]), and r are not independent since (Subscript[r, A])=r+(1/2)R and (Subscript[r, B])=r-(1/2)R.

LCAO Method

Although it is possible to solve the Schrödinger equation for (Subscript[H, 2]^(+)) exactly, it is more convenient to solve the equation approximately. When R is large, the electron will be attached to one of the protons. If the system is in the ground state, and the electron is attached to proton A, the solution will be of the form

Φ(R;r)=(Subscript[ψ, 1,s])((Subscript[r, A])),

where (Subscript[ψ, 1,s])(
          r) = (π^(-1/2)) (e^(-r)) is the normalized ground state wave function for atomic hydrogen. Clearly such a wave function does not have the required symmetry about the midpoint of the internuclear line. The linear combination of atomic orbitals (LCAO) method constructs an approximate solution by forming the even (+) and odd (-) linear combinations,

(Subscript[Φ, ±])(R;r)=(1/(Sqrt[2])) ((Subscript[ψ, 1,s])(
                  (Subscript[r, A])) ± (Subscript[ψ, 1,s])(
                  (Subscript[r, B]))).

We will use these combinations as trial functions in a variational calculation of the energy,

(Subscript[ℰ, ±])(R)=((∫(SuperStar[(Subscript[Φ, ±])(R;r)])ℋ (Subscript[Φ, ±])(R;r)dr)/(∫(SuperStar[(Subscript[Φ, ±])(R;r)]) (Subscript[Φ, ±])(R;r)dr)).

Probability Densities [4]

The electron probability density in the states \(TraditionalForm\`\(Φ\_+\)\) and \(TraditionalForm\`\(Φ\_-\)\) is given by \(TraditionalForm\`|\(Φ\_+\)|\^2\) and \(TraditionalForm\`|\(Φ\_-\)|\^2\), respectively, so that the corresponding charge densities are \(TraditionalForm\`\(ρ\_+\) = \(-e\) |\(Φ\_+\)|\^2\) and \(TraditionalForm\`\(ρ\_-\) = \(-e\) |\(Φ\_-\)|\^2\).

Along the internuclear line (arbitrarily chosen to be the z-axis) the states \(TraditionalForm\`\(Φ\_+\)(R, z)\) and \(TraditionalForm\`\(Φ\_-\)(R, z)\) can be defined by

(Subscript[ψ, 1,s])(
          r_) = (π^(-1/2)) (E^(-r));

(Subscript[Φ, σ])_(R_, 
        z_) = (((Subscript[ψ, 1, s])((Subscript[r, A])) + σ (Subscript[ψ, 1, s])((Subscript[r, B])))/(Sqrt[2])) /. {(Subscript[r, A]) -> \[LeftBracketingBar]z + (R/2)\[RightBracketingBar], 
          (Subscript[r, B]) -> \[LeftBracketingBar]z - 
                (R/2)\[RightBracketingBar]}

where σ = +1 corresponds to \(TraditionalForm\`\(Φ\_+\)\) and σ = -1 corresponds to \(TraditionalForm\`\(Φ\_-\)\).

Confocal Elliptic Coordinates

Integrals which arise in two-centre problems are called two-centre integrals and are most easily evaluated by using confocal elliptic coordinates, defined by

({{ξ = 
              (1/R) ((Subscript[r, A]) + 
                    (Subscript[r, B])),, \(1 ≤ ξ < ∞\)}, {η = 
              (1/R) ((Subscript[r, A]) - 
                    (Subscript[r, B])),, \(\(-1\) ≤ η ≤ 1\)}, {ϕ,, 0 ≤ ϕ < 2  π}})

where ϕ is the usual azimuthal angle about the z-axis. The volume element dr, expressed in confocal elliptic coordinates (ξ, η, ϕ) is given by

\(TraditionalForm\`dr = \(R\^3\/8\) \((ξ\^2 - η\^2)\)\ dξ dη \(\(dϕ\)\(.\)\)\)

Denominator of (Subscript[ℰ, ±])(R) — Normalization [3]

The normalization integral, for the even and odd cases, reads

∫(SuperStar[(Subscript[Φ, ±])(R;r)]) (Subscript[Φ, ±])(R;r)dr=(1/2)∫\((\(\(ψ\_\(1,s\)\)(r\_A)\)\^2 + \(\(\(ψ\_\(1,s\)\)(r\_B)\)\^2 ± 
                    2 \(\( ψ\_\(1,s\)\)(
                        r\_A)\) \(\(ψ\_\(1,s\)\)(
                        r\_B)\)\))\)dr.

Introducing the change of variables,

change = {(Subscript[r, A]) -> (R/2) (ξ + η),  (Subscript[r, B]) -> (R/2) (ξ - η)};

and defining

⟨a_ \[VerticalSeparator] 
          b_⟩ := (( (R^3))/8) (∫_(-1)^1)(∫_1^∞)(∫_0^(2  π))(a b /. 
                      change) ((ξ^2) - (η^2)) \[DifferentialD]ϕ \[DifferentialD]ξ \[DifferentialD]η // 
        Simplify

we can evaluate ∫(((Subscript[ψ, 1 s])((Subscript[r, A])))^2)dr as follows.

⟨(Subscript[ψ, 1,s])(
          (Subscript[r, A])) \[VerticalSeparator] (Subscript[ψ, 1,s])(
          (Subscript[r, A]))⟩

Numerator of (Subscript[ℰ, ±])(R) [2]

From the Schrödinger equation for hydrogen (noting that \(TraditionalForm\`Δ\_r ≡ Δ\_\(r\_A\) ≡ Δ\_\(r\_B\)\)),

{-(1/2) (Subscript[Δ, r]) - (1/r)}(Subscript[ψ, 1,s])(r)=\(ℰ\_\(1,s\)\ \(\(ψ\_\(1,s\)\)(r)\)\), (Subscript[ℰ, 1,s]) = -1/2,

the Schrödinger equation (.) for (Subscript[H, 2]^(+)) yields

\(TraditionalForm\`\(\(ℋ\ \(\(ψ\_\(1 s\)\)(
          r\_A)\) = \((ℰ\_\(1,s\) - 1\/r\_B + 
            1\/R\ )\) \(\(ψ\_\(1 s\)\)(
          r\_A)\)\)\(,\)\)\)

so

(Subscript[ℋ, A,A])=∫(Subscript[ψ, 1,s])((Subscript[r, A]))ℋ (Subscript[ψ, 1,s])((Subscript[r, A]))dr=∫(Subscript[ψ, 1,s])(
                  (Subscript[r, A]))((Subscript[ℰ, 1,s]) - 
                    (1/(Subscript[r, B])) + 
                    (1/R) )(Subscript[ψ, 1,s])((Subscript[r, A]))dr,

and

(Subscript[ℋ, A,B])=∫(Subscript[ψ, 1,s])((Subscript[r, A]))ℋ (Subscript[ψ, 1,s])((Subscript[r, B]))dr=∫(Subscript[ψ, 1,s])(
                (Subscript[r, A]))((Subscript[ℰ, 1,s]) - 
                  (1/(Subscript[r, A])) + (1/R) )(Subscript[ψ, 1,s])(
                (Subscript[r, B]))dr.

(Subscript[ℰ, ±])(R) [6]

. Series [15]

.. Ideal Fermi Gas [9]

The equation of state is a functional relationship between the thermodynamic parameters, say P, V, and T, of the form f(P, V, T) = 0. The parametric form of the equation of state for a spinless ideal Fermi gas is (see Huang, K, Statistical Mechanics, pp187-8, 1987)

\(TraditionalForm\`P\/\(k\ T\) = \(1\/λ\^3\) \(\(f\_\(5/2\)\)(
          z)\), \ 

    1\/v = \(1\/λ\^3\) \(\(f\_\(3/2\)\)(z)\), \)

where v = V/N, \(TraditionalForm\`λ = ℏ \@\(\( 2  π\)\/\(m\ k\ T\)\)\) is the thermal wavelength, and \(TraditionalForm\`\(f\_n\)(
          z) = ∑\+\(l = 1\)\%∞\(\(\((\(-1\))\)\^\(l + 1\)\) z\^l\)\/l\^n\).

.. Cumulant Expansion [6]

The cumulant expansion plays a central role in statistical physics. Let λ be a real parameter and (Subscript[a, n]), (Subscript[b, n]) ∈ ℝ. Assume that

\(TraditionalForm\`\(\(e\^\(∑\+\(n = 1\)\%∞\( b\_n\ λ\^n\)\/\(n!\)\) = ∑\+\(n = 0\)\%∞\(\( a\_n\) λ\^n\)\/\(n!\)\)\(,\)\)\)

where (Subscript[a, 0]) = 1.

. Convergence Acceleration [20]

.. Shank's Transform [14]

One can improve the convergence of the summation of a series by removing its transient behaviour. Suppose that the sum to n terms of a given series satisfies (Subscript[s, n]) = S + α (q^n), (|q| < 1) so that \(TraditionalForm\`s\_n\(⟶\+\(n -> ∞\)\)S\).  To determine S, α and q one requires three equations. These can be provided by computing, e.g., (Subscript[s, n - 1]), (Subscript[s, n]), (Subscript[s, n + 1]).

.. Richardson Extrapolation [6]

One can improve the convergence of the summation of a series by removing its asymptotic behaviour. Suppose that the sum to n terms of a given series satisfies (Subscript[s, n]) = S + ((Subscript[q, 1])/n) + ((Subscript[q, 2])/(n^2)) + … as n -> ∞. Richardson extrapolation consists of truncating at (Subscript[q, p])/(n^p) and solving the p + 1 equations

({{(Subscript[s, n]), =, S + ((Subscript[q, 1])/n) + ((Subscript[q, 2])/(n^2)) + … + 
              ((Subscript[q, p])/(n^p)),}, {(Subscript[s, n + 1]), =, S + ((Subscript[q, 1])/(n + 1)) + 
              ((Subscript[q, 2])/(((n + 1))^2)) + … + 
              ((Subscript[q, p])/(((n + 1))^p)),}, {(Subscript[s, n + 2]), =, S + ((Subscript[q, 1])/(n + 2)) + 
              ((Subscript[q, 2])/(((n + 2))^2)) + … + 
              ((Subscript[q, p])/(((n + 2))^p)),}, {⋮,  , ⋮}, {(Subscript[s, n + p]), =, \(S + q\_1\/\(n + p\) + 
              q\_2\/\((n + p)\)\^2 + … + \(\(q\_p\/\((n + p)\)\^p\)\(.\)\)\)}})

The solution is

\(TraditionalForm\`S\^\([n, p]\) = ∑\+\(k = 0\)\%p\(\(\(\(\(\( 
                          s\_\(n + k\)\)(
                      n + k)\)\^p\) \((\(-1\))\)\^\(k + p\)\)\/\(\(k!\)\ \(\((p - k)\)!\)\)\)\(·\)\)\)

. Green Function Solution to an Inverse Problem [20]

A powerful and general method for solving inhomogenous differential equations such as L[f(x)] = g(x), where L is a differential operator with homogenous boundary conditions, f(a) = 0,  f(b) = 0, is provided by the Green function, G(x, y), which formally satisfies L[G(x, y)] = δ(x - y), where δ is the Dirac delta function (distribution), and the (homogenous) boundary conditions

G(a, y) = 0,  G(b, y) = 0,  G(x, a) = 0,  
    G(x, b) = 0.

Green functions are due to George Green (1793-1841), an English mathematician whose work was not appreciated during his lifetime — largely due to his unusual methodology.

The strength of the method is that G(x, y) is independent of the inhomogenous part of the differential equation, g(x). The solution to the differential equation is

f(
        x) = (∫_a^b) G(x, y) g(
            y) dy.

Consider the differential equation describing a forced simple harmonic oscillator:

(x^″)(t) + x(t) = 
      g(t),

along with the boundary conditions

x(0) = x((π/2)) = 0.

The (region-dependent) Green function, independent of g(t), is

G(t, τ)={({{-cos(τ) sin(t),  
                    0 ≤ 
                      t ≤ τ ≤ (π/2)}, {-cos(t) sin(τ),  
                    0 ≤ τ ≤ 
                      t ≤ (π/2)}})·

Consider the inverse problem of determining the inhomogenous part, g(t), of the differential equation by measuring x(t) at a finite number of times. Discretizing the equation, we obtain the following matrix equation:

\(TraditionalForm\`\(\(x(
        t) = \(∫\_0\%\(π/2\)\ \(G(
              t, τ)\) \(g(τ)\) dτ⟶\(x(t\_i)\) = \(∑\+\(j = 0\)\%m\ \(G(
                t\_i, τ\_j)\) \(g(τ\_j)\) δτ\_j⟶x\_i = ∑\+j\( G\_\(i,j\)\) 
              g\_j\)\)\)\(,\)\)\)

where (Subscript[G, i,j]) = G(
            (Subscript[t, i]), (Subscript[τ, j])) δ(Subscript[τ, j]). The discretized forcing term (Subscript[g, j]) ≡ g((Subscript[τ, j])) can be computed from the inverse of the (discretized) Green function,

\(TraditionalForm\`x\_i = \(∑\+j\( G\_\(i,j\)\) 
              g\_j => 
          g\_j = ∑\+i\( G\_\(i,j\)\^\(-1\)\) \(\(x\_i\)\(·\)\)\)\)

Note that, in general, the matrix (Subscript[G, i,j]) is rectangular because there can be a different number of (Subscript[t, i]) and (Subscript[τ, j]) samples. Hence no inverse matrix exists in the usual sense. However, values of (Subscript[g, j]) = g((Subscript[τ, j])) can be computed using the pseudo or Moore-Penrose inverse (PseudoInverse).