1999 Computational Physics Exam

Instructions

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

The exam consists of 6 questions of approximately equal difficulty and several hints as to how to proceed. You should attempt to answer  whole questions whose total value is 100 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.

. Convergence Acceleration [20]

It is very common in mathematical physics and engineering to expand functions or solve differential equations by expansion into a Fourier basis or into a system of orthogonal polynomials.  However, because these basis functions are well-behaved, convergence near singular points or discontinuities is slow — a famous example being the Gibbs' phenomena for Fourier series. If the form of the singularity is known then one can explicitly take this into account using series re-arrangement to dramatically improve the convergence of series expansions.

As a specific example, consider the generating function for the Legendre polynomials, (Subscript[P, n])(x),

(1/(Sqrt[1 - 2  x t + (t^2)]))= (∑_(n = 0)^∞)(Subscript[P, n])(x)(t^n), -1 < x < 1,|t|<1.

x (Subscript[P, n])(x)=((n (Subscript[P, n - 1])(x)+(n + 1) (Subscript[P, n + 1])(x))/(2 n + 1))·

(1 - x)(∑_(n = 0)^∞)(Subscript[a, n]) (Subscript[P, n])(x)≡  (∑_(n = 0)^∞)(-((n (Subscript[a, n - 1]))/(2 n - 1)) + (Subscript[a, n]) - 
                (((n + 1) (Subscript[a, n + 1]))/(2 n + 3))) (Subscript[P, n])(x).

(1/(Sqrt[1 - x]))=(Sqrt[2]) (∑_(n = 0)^∞)(Subscript[P, n])(x).

. Maxwell Construction [20]

The equation of state of an ideal gas of 𝒩 noninteracting classical particles reads

p V == 𝒩 (Subscript[k, B]) T,

where (Subscript[k, B]) is the Boltzmann constant. A simple approximation, taking into account the interaction between particles, is the van der Waals equation of state,

\(TraditionalForm
    \`\(\(\((p + a\/V\^2)\) \((V - b)\) == 
      𝒩\ \(k\_B\) T\)\(,\)\)\)

with parameters a and b.  As you will see, there exists a critical temperature (Subscript[T, c]) below which loops appear in the p - V diagram.  These loops are forbidden for thermodynamical reasons. For T < (Subscript[T, c]), a phase transition between gas and liquid takes place; at T = (Subscript[T, c]) the difference between the liquid and gaseous phases vanishes.

(Subscript[T, c]) and (Subscript[V, c]) are determined by the condition that the first and second derivatives of p(V, T) with respect to V vanish, i.e., ((∂p(V, T))/(∂V))==0 and (((∂^2) p(V, T))/(∂(V^2)))==0.

\(TraditionalForm
    \`∫\_v\_1\%v\_3\(\( p\&~\)(v)\) dv = 
      \((v\_3 - v\_1)\) \(\(\(p\_t\)\&~\)\(·\)\)\)

({{\(\(p\&~\) \((v\_1, t)\)\), ==, \(\(\(p\&~\)(v\_3, t)\)\(,
            \)\)}, {\(∫\_v\_1\%v\_3\(\( p\&~\)(v, t)\) dv\), ==, \(\((v\_3 - v\_1)\) \(\(\(p\&~\)(v\_1, t)\)\(.\)\)\)}})

. Normal Modes of a Circular Membrane [20]

If there is no external force, the (wave) equation for free vibrations of a uniform isotropic membrane reads

(((∂^2) u)/(∂(t^2)))==(c^2) Δ u,

where Δ is the two-dimensional Laplacian. For a circular membrane of radius 1 fixed at its edge, polar coordinates (r, θ) are convenient; the boundary condition is u(1, θ, t) == 0 and the Laplacian reads

\(TraditionalForm
    \`Δ[f_] := 
      \(1\/r\) ∂\((r ∂f\/∂r)\)\/∂r + 
        \(1\/r\^2\) ∂\^2 f\/∂θ\^2\)

. Ground State Energy of Helium [20]

In atomic units, the Hamiltonian for a two-electron atom with nuclear charge (SuperPlus[𝒵]) reads

\(TraditionalForm
    \`\(\(ℋ = 
      \(-\(1\/2\)\) \((Δ\_1 + Δ\_2)\) - 
        𝒵\/r\_1 - 𝒵\/r\_2 + 1\/r\_12\)\(,
    \)\)\)

where the interparticle coordinates —  (Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 12])— are defined in Figure 1.

$FormatType = TraditionalForm;

$TextStyle = {FontSize -> 10};

\(TraditionalForm
    \`\(Show[Graphics[
          {Arrow[{0, 0}, \ {1, 0.1}]\ , Arrow[{0, 0}, \ {0.7, 1}], 
            Arrow[{0.7, 1}, {1, 0.1}], AbsolutePointSize[40], \ 
            Point[{0, 0}], AbsolutePointSize[20], Point[{1, 0.1}], 
            Point[{0.7, 1}], Text[r\_1, \ {0.3, 0.55}], 
            Text[r\_2, \ {0.5, 0}], Text[r\_12, \ {0.9, 0.55}], GrayLevel[1], 
            Text[1, \ {0.7, 1}], Text[2, \ {1, 0.1}], 
            Text[\(𝒵\^+\), \ {0, 0}]}], 
        AspectRatio -> Automatic, 
        PlotRange -> {{\(-0.1\), 1.1}, {\(-0.1\), 1.1}}];\)\)

In interparticle coordinates, the matrix element of a general operator, ℒ, is defined by the following integral:

\(TraditionalForm
    \`⟨ϕ 
          |ℒ| 
          ϕ⟩ = 
      8 \( π\^2\) 
        \(∫\_0\%∞\( r\_1\) dr\_1 
            \(∫\_0\%∞\( r\_2\) dr\_2 
                \(∫\_|r\_1 - r\_2|\%\(r\_1 + r\_2\)\(r\_12\) \(ϕ\^*\) 
                    ℒ\ ϕ 
                    \(\(dr\_12\)\(·\)\)\)\)\)\)

In assignment 6 you showed that for

ϕ((Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 12]))=(e^(-(Subscript[λ, 1]) (Subscript[r, 1]) - (Subscript[λ, 2]) (Subscript[r, 2]) 
                - (Subscript[λ, 12]) (Subscript[r, 12]))),

the matrix element of an arbitrary monomial, (Subscript[r, 1]^i) (Subscript[r, 2]^j) (Subscript[r, 12]^k) with i, j, k ≥  -1, can be computed by differentiating a particular matrix element, viz,

\(TraditionalForm
    \`\(\(\(∂\^\(i + j + k + 3\)\/\(∂λ\_1\%\(i + 1\) ∂λ\_2\%\(j + 1\) 
              ∂λ\_12\%\(k + 1\)\)\) 
        ⟨ϕ 
            |1\/\(\(r\_1\) \(r\_2\) r\_12\)| ϕ⟩ = 
      \(\((\(-2\))\)\^\(i + j + k + 3\)\) 
        ⟨ϕ 
            |\(r\_1\%i\) \(r\_2\%j\) 
                r\_12\%k| ϕ⟩\)\(,
    \)\)\)

where

\(TraditionalForm
    \`⟨ϕ 
          |1\/\(\(r\_1\) \(r\_2\) r\_12\)| ϕ⟩ = 
      \(\(\(2  π\^2\)\/\(\((λ\_1 + λ\_2)\) 
              \((λ\_1 + λ\_12)\) 
              \((λ\_2 + λ\_12)\)\)\)\(·\)\)\)

The importance of this result is that matrix elements of any combination of integral powers (≥-1) of the interparticle coordinates can be obtained by differentiation with respect to the parameters, (Subscript[λ, i]).

In interparticle coordinates the (s-state) kinetic energy operator, \(TraditionalForm
      \`𝒯 = 
        \(-\(1\/2\)\) \((Δ\_1 + Δ\_2)\)\), reads

𝒯(ψ_):=-(1/2) ((((∂^2) ψ)/(∂((Subscript[r, 1])^2)))+(((∂^2) ψ)/(∂((Subscript[r, 2])^2)))+2(((∂^2) ψ)/(∂((Subscript[r, 12])^2)))+(2/(Subscript[r, 1]))((∂ψ)/(∂(Subscript[r, 1])))+(2/(Subscript[r, 2]))((∂ψ)/(∂(Subscript[r, 2])))+(4/(Subscript[r, 12]))((∂ψ)/(∂(Subscript[r, 12])))+((((Subscript[r, 1]^2) - (Subscript[r, 2]^2) + (Subscript[r, 12]^2)))/((Subscript[r, 1]) (Subscript[r, 12])))(((∂^2) ψ)/(∂(Subscript[r, 1]) ∂(Subscript[r, 12])))+((((Subscript[r, 2]^2) - (Subscript[r, 1]^2) + (Subscript[r, 12]^2)))/((Subscript[r, 2]) (Subscript[r, 12])))(((∂^2) ψ)/(∂(Subscript[r, 2]) ∂(Subscript[r, 12]))))

the potential energy operator, 𝒱, is

\(TraditionalForm
    \`\(𝒱 = 
        \(-\(𝒵\/r\_1\)\) - 𝒵\/r\_2 + 1\/r\_12;
      \)\)

and the Hamiltonian, ℋ, is

\(TraditionalForm
    \`ℋ(ψ_) := 
      𝒯(ψ) + 𝒱\ ψ\)

In Mathematica notation, the matrix element \(TraditionalForm
      \`⟨i, j, k⟩ ≡ 
        ⟨ϕ 
            |\(r\_1\%i\) \(r\_2\%j\) 
                r\_12\%k| ϕ⟩\) can be computed and saved using the following code:

\(TraditionalForm
    \`⟨i_, j_, k_⟩ := 
      \(⟨i, j, k⟩ = 
        \(1\/\((\(-2\))\)\^\(i + j + k + 3\)\) 
          ∂\_\({λ\_1, i + 1}, {λ\_2, j + 1}, {λ\_12, k + 1}\)\(\(2  π\^2\)\/\(\((λ\_1 + λ\_2)\) 
                  \((λ\_1 + λ\_12)\) 
                  \((λ\_2 + λ\_12)\)\)\)\)\)

Using this result we can compute matrix elements of any zeroth order approximation consisting of a linear combination of terms of the form  \(TraditionalForm
      \`\(P(r\_1, r\_2, r\_12)\) 
        e\^\(\(-λ\_1\) r\_1 - \(λ\_2\) r\_2 - \(λ\_12\) r\_12\)\), where P((Subscript[r, 1]), (Subscript[r, 2]), (Subscript[r, 12])) is a polynomial in the interparticle coordinates.

Consider the trial wavefunction \(TraditionalForm
      \`ψ\_0 = 
        e\^\(\(-𝒵\)\ r\_1 - 𝒵\ r\_2 + r\_12/2\)\).

\(TraditionalForm
    \`\(p = 1 + 2\/\(\(a\)\(\ \)\) + 3\/b + 4  a + 5  a\ b;\)\)

Expand[(a^i) (b^j) p]

\(TraditionalForm
    \`% /. \(a\^\(i + m_. \)\) b\^\(j + n_. \) -> 
        ⟨m, n⟩\)

. Pauli spin matrices and Dirac gamma matrices [20]

The Pauli spin matrices,

σ/:(Subscript[σ, 1])=(({{0, 1}, {1, 0}}));

σ/:(Subscript[σ, 2])=(({{0, -I}, {I, 0}}));

σ/:(Subscript[σ, 3])=(({{1, 0}, {0, -1}}));

along with the 2×2 identity matrix:

(Subscript[ℐ, n])_ := IdentityMatrix[n]

(Subscript[ℐ, 2])

are used in a matrix representation of spin operators in quantum mechanics.

((Subscript[a, 1]) (Subscript[σ, 1]) + (Subscript[a, 2]) (Subscript[σ, 2]) + 
              (Subscript[a, 3]) (Subscript[σ, 3])) . 
          ((Subscript[b, 1]) (Subscript[σ, 1]) + (Subscript[b, 2]) (Subscript[σ, 2]) + 
              (Subscript[b, 3]) (Subscript[σ, 3]))=={(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])} . {(Subscript[b, 1]), (Subscript[b, 2]), (Subscript[b, 3])} (Subscript[ℐ, 2])+i ({(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3])}×{(Subscript[b, 1]), (Subscript[b, 2]), (Subscript[b, 3])}).{(Subscript[σ, 1]), (Subscript[σ, 2]), (Subscript[σ, 3])}. [3]

(Subscript[γ, k])=(({{𝒪, -i (Subscript[σ, k])}, {i (Subscript[σ, k]), 𝒪}})), k = 1,2,3,
(Subscript[γ, 4])=(({{(Subscript[ℐ, 2]), 𝒪}, {𝒪, -(Subscript[ℐ, 2])}})),
(Subscript[γ, 5]) =  
          (Subscript[γ, 1]) . (Subscript[γ, 2]) . (Subscript[γ, 3]) . (Subscript[γ, 4]),

. Series [20]

.. Discrete nonlinear transform [12]

Consider the formal Taylor series

\(TraditionalForm
    \`\(\(∑\+\(n = 0\)\%∞\(\((\(-1\))\)\^n\) \(a\_n\) x\^n\)\(,
    \)\)\)

where (Subscript[a, 0]) = 1and the continued fraction

\(TraditionalForm
    \`\(\(1\/\(1 + 
            \(\(b\_1\) x\)\/\(1 + \(\(b\_2\) x\)\/\(1 + \(\(b\_3\) x\)\/\(1 + …\)\)\)\)\)\(·\)\)\)

.. Kepler equation [8]

The Kepler equation, s = u + e sin s, arises in celestial mechanics. Defining s = u + a reduces the equation to a = e sin(u + a) which can be formally solved by introducing the recurrence relation (Subscript[a, k]) = e sin(u + (Subscript[a, k - 1])). Since e is a small parameter, as k -> ∞, (Subscript[a, k - 1]) -> (Subscript[a, k]) -> a, and hence the Kepler equation is formally satisfied.