Solving Advanced Physics Problems 
with Mathematica

Michael Trott, Wolfram Research, Inc.

Introduction

● Mathematica is a very valuable tool to solve theoretical physics problems

● Each version of Mathematica makes it more useful and opens new areas 
   (for instance, Version 5: sparse array capabilities)

● Present here a selection of hopefully interesting examples (with lots of visualizations and few formulas)

● All code, complete explanations, detailed references are given in the Mathematica GuideBook for Programming and 
   the Mathematica GuideBook for Graphics (Springer-Verlag 2004; redistributed with permission)
   and additions to these books on http://www.MathematicaGuideBooks.org

What makes Mathematica useful to solve physics problems?

What Mathematica features are important?

At one time or another every mathematica functions comes in handy

Most frequently and virtually irreplacable functions:

• Significance arithmetic

• Differentiation, integration of simple functions (potentially milions of times)

• Resultant, GroebnerBasis, CylindricalDecomposition

• Together (to keep things quick and tidy)

• Fourier, FindMinimum, NDSolve

• Sparse linear algebra

• Special functions of mathematical physics

And, of course,

• Programming language and smooth reuse of results in numerics, symbolics, graphics

What Mathematica features are not important?

• Integration of complicated functions

• Functions with unpredictable runtime complexity (like Simplify and FullSimplify)

• Buttons and palettes of any form

What is different from paper-and-pencil calculations?

• More larger steps (⟹programming) and fewer smaller steps required in Mathematica

• Use more series expansions (for instance for the eigenvalues of the harmonic oscillator)

• Different views on the problem are more easily investigated

• Manipulate and carry around larger expressions

• Adapt to uniquely determined branch cuts for contour integrals

• Try to rephrase as an algebraic problems, rather than an analysis problem

What to keep in mind?

• Write version-independent code (do not rely on the form of integrals or the order of solutions from Solve)

• Do all steps programmatically (never edit output!)

• Compare theoretical complexity with achieved runtime complexity

• Comment your code

Areas not to tackle—yet

If you do not want to write larger amounts of packages

• Integration-related parts of differential geometry (Stokes theorem)

• Functional integration-related problems

• Local quantum physics

Ten amuse bouche (of recreational taste)

Hofstadter butterflies

Problem

Energy spectrum of a quantum-mechanical particle under the simultaneous influence of a 2D periodic potential and a perpendicular magnetic field.

(((("p")^2)/(2m))+V("r"))(Subscript[Ψ, E])("r")=E(Subscript[Ψ, E])("r")

"A"(ϵ,m,α,ν)=(({{ϵ-2cos(2π m α-ν), -1}, {1, 0}}))

Q(ϵ,α,q)=(∏_(k=0)^(q-1))"A"(ϵ,k,α,(π/(2q))).

Core Mathematica code

Function[α, {Line[{{#[[1]], α}, {#[[2]], α}}], 
              Line[{{#[[1]], 1 - α}, {#[[2]], 1 - α}}]}& /@
  Partition[Sort[Re[#]& /@ Select[Level[
  {NRoots[# == -4, C], NRoots[# == +4, C]}, {-1}], NumberQ]], 2]&[
  (#[[1, 1]] + #[[2, 2]])&[Fold[Expand[#1.#2]&, First[#], Rest[#]]&[
  Function[x, N[{{C - 2 Cos[2Pi x α - Pi/(2 Denominator[α])], -1},
                {1, 0}}], {Listable}][Range[0, Denominator[α] - 1]]]]]] /@
Union[Flatten[Range[#/2]/#& /@ Range[2, maxDenominator]]]}]

Results

Organic shapes from simple rules

Problem

The time evolution of a a system of coupled cells on a 2D lattice.

(Subscript[c, i,j]^((n+1)))=(1-(Subscript[ε, 1])-(Subscript[ε, 2]))f((Subscript[c, i,j]^((n))))+((Subscript[ε, 2])/(d^2))(∑_(i,j=1)^d)f((Subscript[c, i,j]^((n))))+(Subscript[ε, 1])((1/6)(f((Subscript[c, i+1,j]^((n))))+f((Subscript[c, i-1,j]^((n))))+f((Subscript[c, i,j+1]^((n))))+f((Subscript[c, i,j+1]^((n)))))+(1/12)(f((Subscript[c, i+1,j+1]^((n))))+f((Subscript[c, i+1,j-1]^((n))))+f((Subscript[c, i-1,j+1]^((n))))+f((Subscript[c, i-1,j-1]^((n)))))).

Core Mathematica code

Compile[{{L, _Real, 2}, ε1, ε2, a},
Module[{l = 1. - a L^2, n = Length[L]},
(* nonlinear local term *)
(1 - ε1 - ε2) l +
(* nearest neighbor interactions *)
ε1 (ListConvolve[{{0., 1., 0.}, {1., 0., 1.}, {0., 1., 0.}}, l, {2, 2}]/6 + 
    ListConvolve[{{1., 0., 1.}, {0., 0., 0.}, {1., 0., 1.}}, l, {2, 2}]/12) +
(* global average *)
ε2/n^2 Sum[l[[i, j]], {i, n}, {j, n}]]];

Results

Closed orbits of the three-body problem

Problem

Find periodic solutions of the restricted three-body problem

ξ''(t)=ξ(t)+2 η'(t)-(1-μ)(((ξ(t)+μ))/( (Subscript[ρ, 1])))-μ (((ξ(t)-(1-μ)))/(Subscript[ρ, 2]))
η''(t)=η(t)-2 ξ'(t)-(1-μ) ((η(t))/(Subscript[ρ, 1]))-μ ((η(t))/(Subscript[ρ, 2]))
(Subscript[ρ, 1])=((((((ξ(t)+μ))^2)+((η(t))^2)))^(3/2))
(Subscript[ρ, 2])=((((((ξ(t)-(1-μ)))^2)+((η(t))^2)))^(3/2)).

Core Mathematica code

(* solve the differential equations for given initial conditions *)
solveODEs[{ξ0_, Η0_}, T_, opts___] :=
Block[{(* moon-earth value *) μ = 0.01227747},
eqs = Module[{ν = 1 - μ, ρ1, ρ2},
{ρ1, ρ2} = ((ξ[t] + {μ, -ν})^2 + η[t]^2)^(3/2);
(* use first order system *)
{ξ'[t] == Ξ[t], η'[t] == Η[t],
 Ξ'[t] == ξ[t] + 2 Η[t] - ν(ξ[t] + μ)/ρ1 - μ(ξ[t] - ν)/ρ2,
 Η'[t] ==  η[t] - 2 Ξ[t] - ν η[t]/ρ1 - μ η[t]/ρ2}];
Check[NDSolve[Flatten[{eqs, 
     ξ[0] == ξ0, Ξ[0] == 0, η[0] == 0, Η[0] == Η0}],
     {ξ, η, Ξ, Η}, {t, 0, T}, opts, MaxSteps -> 10^4,
     PrecisionGoal -> 10, AccuracyGoal -> 10], $Failed]]

(* determine minimum distance in initial condition space *)
minDiff[{ξ0_, Η0_}, T_, nsol_] := 
Module[{δ, pl, δs, fm},
 (* summed squared difference to initial conditions *)
 δ[t_?NumberQ] = (ξ[t] - ξ0)^2 + η[t]^2 +
                  Ξ[t]^2 + (Η[t] - Η0)^2 /. nsol[[1]];
 (* find global minimum using Plot's adaptive refinement *)
 pl = Plot[Evaluate[δ[t]], {t, 5, T},
        AxesOrigin -> {0, 0}, PlotRange -> All,
        PlotStyle -> {Hue[0]}, DisplayFunction -> Identity];
 points = Cases[pl, _Line, Infinity][[1, 1]];     
 δs = points[[Position[points, Min[Last /@ points]][[1, 1]], 1]];
 (* refine minimum if possible *)
 Off[InterpolatingFunction::dmval];
 If[δs < 6., $Failed,
    fm = FindMinimum[δ[t], {t, δs, δs + 10^-3}];
    {fm[[2, 1, 2]], fm[[1]]}]]

minDiff[{ξ0_, Η0_}, T_, $Failed] = $Failed;

(* estimate length of the orbit *)
orbitLengthEstimate[nsol_, T_] := 
(* some length of discretized orbit *)
Plus @@ (Sqrt[#.#]& /@ (Subtract @@@ 
      Partition[Table[Evaluate[{ξ[t], η[t]} /. nsol[[1]]],
              {t, 0, T, 10^-2}], 2, 1]));

Results

Quantum carpets

Problem

Time-evolution of the probability density of a particle in a square well with homogeneous initial conditions

i ((∂ψ(ξ,τ))/(∂τ))=-(1/(4π))(((∂^2)ψ(ξ,τ))/(∂(ξ^2)))
ψ(0,τ)=ψ(1,τ)=0
ψ(ξ,0)=1    for  0≤τ≤1

Core Mathematica code

ψcf = Compile[{{ξ, _Real}, {τ, _Real}, {n, _Integer}, {m, _Integer}}, 
Module[{ψ0}, 
       ψ0 = Sum[1/(2k + 1) Sin[(2k + 1) Pi ξ] Exp[-I Pi/4 (2k + 1)^2 τ], 
               {k, 0, n}];
       ψ0 + FoldList[Plus, 0, Table[1/(2k + 1)*
                     Sin[(2k + 1) Pi ξ] Exp[-I Pi/4 (2k + 1)^2 τ], 
                            {k, n + 1, n + m}]]]];

data = Table[Table[SequenceLimit[ψcf[ξ, τ, 1000, 50]],
              {ξ, 0, 1/2, 1/2/ppξ}], {τ, 0, 1, 1/ppτ}];

Results

High-order perturbation theory

Problem

Carry out the Rayleigh–Schrödinger perturbation theory up to order 200 and resum the resulting diverging series for the quartic oscillator.

((Subscript[(OverHat[L]), 0])+g (OverHat[V]))(Subscript[ψ, n]^((g)))(z)=(Subscript[ε, n]^((g)))(Subscript[ψ, n]^((g)))(z)

(Subscript[(OverHat[L]), 0])(Subscript[ψ, n]^((0)))(z)=(Subscript[ε, n]^((0)))(Subscript[ψ, n]^((0)))(z)

(Subscript[ε, n]^((g)))=(∑_(j=0)^∞)(g^j)  (^((n)))(Subscript[ε, j])

(Subscript[ψ, n]^((g)))(x)=(∑_(j=0)^∞)(g^j) ((∑_(k=0)^∞)(^((n)))(Subscript[c, j,k])(Subscript[ψ, k]^((0)))(z)).

(^((n)))(Subscript[ε, j])=(∑_(l=0)^∞)(Subscript[V, m,l])((Subscript[𝒩, l])/(Subscript[𝒩, m]))(^((n)))(Subscript[c, j-1,l])

(^((n)))(Subscript[c, j,k])=(1/((Subscript[ε, k]^((0)))-(Subscript[ε, n]^((0)))))((∑_(l=1)^(j-1))(^((n)))(Subscript[ε, l])(^((n)))(Subscript[c, j-l,k])+(∑_(l=0)^∞)(Subscript[V, k,l])((Subscript[𝒩, l])/(Subscript[𝒩, k]))(^((n)))(Subscript[c, j-1,l])).=

-(Subscript[ψ, n]^″)(z)+((z^2)+(z^4))(Subscript[ψ, n])(z)=(Subscript[ε, n])(Subscript[ψ, n])(z)

Core Mathematica code

((* remember values of V *)
 𝒱[n_, m_] := 𝒱[n, m] = V[n, m];
 (* the recursive definitions for ϵ and c *)
 ϵ[0] = N[1, prec];
 ϵ[i_] := ϵ[i] = SetPrecision[Sum[𝒱[0, k] c[i - 1, k], {k, -4, 4, 2}], prec];
 c[0, 0] = N[1, prec];
 c[i_, 0] = 0;
 c[0, k_] = 0;
 c[i_, k_] := (c[i, k] = 0) /; k > 4i;
 c[i_, k_] := c[i, k] = (* use prescribed precision at each step *) 
   SetPrecision[
   (Sum[ϵ[j] c[i - j, k], {j, i - 1}] -
    Sum[𝒱[k, n] c[i - 1, n], {n, k - 4, k + 4, 2}])/(2k), prec])

formalPowerSeriesToDiagonalPadé[as_, M_, x_] := 
Det[Table[If[j != M, as[[i + j + 2]], 
             Sum[as[[k - M + i + 1]] x^k, {k, M - i, M}]],
          {j, 0, M}, {i, 0, M}]]/
Det[Table[If[j != M, as[[i + j + 2]], x^(M - i)], {j, 0, M}, {i, 0, M}]]

Results

A sliding spinning coin

Problem

The equations of motion for a spinning sliding coin on a table are:

(Overscript[𝓍, ..])(τ)={({{((4 )/(3 π ε(τ)))((((ε(τ))^2)+1) E(((ε(τ))^2))+((ε^2)-1) K(((ε(τ))^2))), if |ε(τ)|≤1}, {((4 sgn(ε(τ)))/(3 π))  ((((ε(τ))^2)+1) E(((ε(τ))^(-2)))-(((ε(τ))^2)-1) K(((ε(τ))^(-2)))), if |ε(τ)|≥1}})

(Overscript[φ, ..])(τ)={({{(4/(9 π)) ((4-2 ((ε(τ))^2)) E(((ε(τ))^2))+(((ε(τ))^2)-1) K(((ε(τ))^2))), if |ε(τ)|≤1}, {((4 sgn(ε(τ)))/(9π ε(τ))) ((4 ((ε(τ))^2)-2 ((ε(τ))^4)) E(((ε(τ))^(-2)))+(2 ((ε(τ))^4)-5 ((ε(τ))^2)+3) K(((ε(τ))^(-2)))), if |ε(τ)|≥1}})

ε(τ)=(OverDot[𝓍])(τ)/(OverDot[φ])(τ)

Core Mathematica code

ℱ[ε_] = 
Which[ε == 1, 8/(3 Pi), (* no rotation *) ε === ComplexInfinity, 1, 
Abs[ε] < 1,  
4/3 ((ε^2 + 1) EllipticE[ε^2] + (ε^2 - 1) EllipticK[ε^2])/(Pi ε),
Abs[ε] > 1,
4/3 Sign[ε] ((ε^2 + 1) EllipticE[1/ε^2] - (ε^2 - 1) EllipticK[1/ε^2])/(Pi)];

𝒯[ε_] = 
Which[ε == 1, 8/(9 Pi), (* no rotation *) ε === ComplexInfinity, 0, 
Abs[ε] < 1, 
4/9 ((4 - 2 ε^2) EllipticE[ε^2] + (ε^2 - 1) EllipticK[ε^2])/Pi,
Abs[ε] > 1,
4/9 Sign[ε] ((4 ε^2 - 2 ε^4) EllipticE[1/ε^2] + 
     (2ε^4 - 5ε^2 + 3) EllipticK[1/ε^2])/(Pi ε)];

Results

Operator splitting

Problem

Find numbers (Subscript[ω, k]), such that

(e^(t((OverHat[A])+(OverHat[B]))))= (e^(t (Subscript[ω, 1])(OverHat[A])))(e^(t (Subscript[ω, 2])(OverHat[B])))(e^(t (Subscript[ω, 3])(OverHat[A])))(e^(t (Subscript[ω, 4])(OverHat[B])))(e^(t (Subscript[ω, 5])(OverHat[A])))(e^(t (Subscript[ω, 6])(OverHat[B])))(e^(t (Subscript[ω, 7])(OverHat[A])))+O((t^5))

Core Mathematica code

a___ ⊗ (⊗b___) ⊗ c___ := a ⊗ b ⊗ c; 
a___ ⊗ (b_ + c_) ⊗ d___ := a ⊗ b ⊗ d + a ⊗ c ⊗ d; 
⊗ x_?notOperatorQ := x; 
a___ ⊗ x_?notOperatorQ ⊗ b___ := x (a ⊗ b); 
a___ ⊗ (x__?notOperatorQ y_) ⊗ b___ := x(a ⊗ y ⊗ b); 
a___ ⊗ b_Plus^n_ ⊗ c___ := (a ⊗ ## ⊗ c)& @@ Table[b, {n}]; 
a___ ⊗ OverHat[A_]^n_. ⊗ OverHat[A_]^(m_.) ⊗ c___ := 
       a ⊗ OverHat[A]^(n + m) ⊗ c;

Results

Time-evolution in a smooth three-well potential

Problem

Visualize the time-development of a an initially localized wave packet in a cos-triple well potential.

i((∂Ψ(x,t))/(∂t))=-(((∂^2)Ψ(x,t))/(∂(x^2)))-(Subscript[V, 0]) cos(x) Ψ(x,t)
Ψ(-3π,t)=Ψ(3π,t)
Ψ(x,0)=cos(𝒦 x)(cos^2)((( π )/2)(x/(Subscript[x, M]))) θ((Subscript[x, M]^2)-(x^2))

Core Mathematica code

o = 160;
εData = Table[MathieuCharacteristicA[k/3, -4 V0]/8 // N[#, 40]&,
{k, 0, 2o, 2}];

pp = 600;
ψData = Table[Table[{x, ψ[x, N[εData[[k]]]]}, 
                    {x, 0, 3Pi, 3Pi/pp}], {k, o}];

Do[ψIpo[k] = Interpolation[ψData[[k]], InterpolationOrder -> 16], {k, o}]

Results

Brain growth model

Problem

Solve the following PDE system for random initial conditions:

((∂ψ(t,x,y))/(∂t))==(((∂^2)ψ(t,x,y))/(∂x ∂x))+(((∂^2)ψ(t,x,y))/(∂y ∂y))+γ (-(1/3) ((ψ(t,x,y))^3)+ψ(t,x,y)+η(t,x,y))

((∂η(t,x,y))/(∂t))==δ ((((∂^2)η(t,x,y))/(∂x ∂x))+(((∂^2)η(t,x,y))/(∂y ∂y)))-((ν+β η(t,x,y)+ψ(t,x,y))/γ)

Core Mathematica code

Module[{ndsol},
(* solve PDE system *)
ndsol = NDSolve[Flatten[{eqs1[{γ, β, ν, δ}], 
         ψ[0, x, y] == ric[nψ, mψ, seedψ][x/L, y/L], 
         η[0, x, y] == ric[nη, mη, seedη][x/L, y/L],
         ψ[t, -L, y] == ψ[t, +L, y], ψ[t, x, -L] == ψ[t, x, +L],
         η[t, -L, y] == η[t, +L, y], η[t, x, -L] == η[t, x, +L]}],
         {ψ, η}, {t, 0, T}, {x, -L, L}, {y, -L, L}, 
         Sequence @@ {ndsolveOpts}, MaxSteps -> 10^5,
         Method -> {"MethodOfLines", "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> do, 
                    "MaxPoints" -> {ppL, ppL}, "MinPoints" -> {ppL, ppL}}}]]

Results

Eigenmodes of a vibrating tetrahedron

Problem

Visualize the eigenmodes of the Helmholtz equation in a tetrahedron with Dirichlet boundary conditions.

(Subscript[ψ, l,m,n])(x,y,z)=|(({{1, 1, 1, 1}, {(Subscript[g, l])(x,z,3), (Subscript[g, l])(y,-z,1), (Subscript[g, l])(-x,z,-1), (Subscript[g, l])(-y,-z,-3)}, {(Subscript[g, m])(x,z,3), (Subscript[g, m])(y,-z,1), (Subscript[g, m])(-x,z,-1), (Subscript[g, m])(-y,-z,-3)}, {(Subscript[g, n])(x,z,3), (Subscript[g, n])(y,-z,1), (Subscript[g, n])(-x,z,-1), (Subscript[g, n])(-y,-z,-3)}}))|

(Subscript[g, α])(ξ,η,j)=exp(i α (ξ /(2^(1/2))+η/2+j π/4))

Core Mathematica code

xyzTo𝓈𝓉𝓊[{𝓈_, 𝓉_, 𝓊_}] = 
({x, y, z} /. Solve[{x + y == 𝓈 (1 - 𝓊)/Sqrt[2] Pi , 
                     x - y == 𝓉 (𝓊 + 1)/Sqrt[2] Pi , 
                     z     == 𝓊 Pi/2}, {x, y, z}][[1]])

interiorNodalSurface[{l_, m_, n_}, reIm_, opts___] :=
Module[{ε = 10^-8, tetrahedronEdges, tetrahedronOutline, addColor},
 (* contour plot in slice of 𝓈,𝓉,𝓊-space *)
 cp𝓈𝓉𝓊 = ContourPlot3D[Evaluate[reIm @ ψ[{l, m, n}, xyzTo𝓈𝓉𝓊[{𝓈, 𝓉, 𝓊}]]],
            {𝓈, -1 + ε, 1 - ε}, {𝓉, -1 + ε, 1 - ε}, {𝓊, -1 + ε, 1 - ε},
            Contours -> {0}, PlotPoints -> {14, 3},
            MaxRecursion -> 1, DisplayFunction -> Identity];
 (* triangulated faces of the tetrahedron *)
 tetrahedronOutline = {Thickness[0.003], GrayLevel[0], 
 Map[Line[Append[#, First[#]]]&, Nest[
     Flatten[Apply[{{#1, (#1 + #2)/2, (#1 + #3)/2}, 
                    {(#1 + #2)/2, #2, (#2 + #3)/2},
                    {(#2 + #3)/2, #3, (#1 + #3)/2}, 
                    {(#1 + #2)/2, (#2 + #3)/2, (#1 + #3)/2}}&, 
                                               #, {2}], 1]&,
            N[{{{p1, p2, p3}, {p1, p2, p4},
                {p1, p3, p4}, {p2, p3, p4}}}], 2], {-3}]};
 (* color nodal surfaces according to distance from origin *)
 addColor[Polygon[ℓ_]] := SurfaceColor[#, #, 2.4]&[
                          Hue[0.23 Sqrt[#.#]&[Plus @@ ℓ/Length[ℓ]]]]]

Show[GraphicsArray[#]]& /@
 With[{ins = interiorNodalSurface[##, DisplayFunction -> Identity]&},
     (* six examples *)
     {{ins[{2, 4, 7}, Re], ins[{1, 7, 9}, Re], ins[{3, 7,  8}, Im]},
      {ins[{1, 5, 7}, Im], ins[{3, 7, 8}, Re], ins[{1, 7, 15}, Im]}}]

Results

Ten sample case studies

Better check more terms: Generalized WKBK approximation

Problem description

The Schrödinger equation can be rewritten as a Riccati equation, which in turn can be differentiated N times and backsubstitution yields a higher-order differential equation. A generalized WKBK quantization condition can be naturally stated (T. Hyouguchi, R. Seto, M. Ueda, S. Adachi. Ann. Phys. 312, 177 (2004))

-ψ''(x)+V(x)ψ(x)=ε ψ(x)

φ''(x)+((φ'(x))^2)=(Subscript[𝒦, ε])(x)   with  (Subscript[𝒦, ε])(x)=ε-V(x)

(1/(2π i))(∮_𝒞)y((Subscript[𝒦, (Subscript[ε, n])]^((k)))(x))dx=n+1-(1/2)(Subscript[δ, N,2])

(φ^((N)))(x)+(Subscript[𝓅, N])(φ'(x),(Subscript[𝒦, ε]^((k)))(x))

(Subscript[𝓅, N])(y((Subscript[𝒦, (Subscript[ε, n])]^((k)))(x)),(Subscript[𝒦, (Subscript[ε, n])]^((k)))(x))=0

The conjecture that all higher order quantization conditions are correct for the harmonic oscillator appears naturally after checking that the lowest orders yield the exact result.

Why is this problem interesting?

• Any quantum mechanics-related problem is interesting
• Simple quantization conditions are useful for many applications

Solution strategy

• Derive coupled nonlinear system of ODEs equivalent to the contour integral
• Integrate ODE system with high precision

Crucial Mathematica functions and features

• High-precision arithmetic

• Resultant, NDSolve

Results

Differentiated Riccati ODEs:

{((φ'(x))^2)+𝒦(x)+φ''(x),
-2 ((φ'(x))^3)-2 𝒦(x) φ'(x)+𝒦'(x)+Derivative[3][φ](x),
6 ((φ'(x))^4)+8 𝒦(x) ((φ'(x))^2)-2 𝒦'(x) φ'(x)+2 ((𝒦(x))^2)+𝒦''(x)+Derivative[4][φ](x),-24 ((φ'(x))^5)-40 𝒦(x) ((φ'(x))^3)+10 𝒦'(x) ((φ'(x))^2)-16 ((𝒦(x))^2) φ'(x)-2 𝒦''(x) φ'(x)+6 𝒦(x) 𝒦'(x)+Derivative[3][𝒦](x)+Derivative[5][φ](x),120 ((φ'(x))^6)+240 𝒦(x) ((φ'(x))^4)-60 𝒦'(x) ((φ'(x))^3)+136 ((𝒦(x))^2) ((φ'(x))^2)+12 𝒦''(x) ((φ'(x))^2)-52 𝒦(x) 𝒦'(x) φ'(x)-2 Derivative[3][𝒦](x) φ'(x)+16 ((𝒦(x))^3)+6 ((𝒦'(x))^2)+8 𝒦(x) 𝒦''(x)+Derivative[4][𝒦](x)+Derivative[6][φ](x)}

Contour integral integrands:

{-(x^2)+(y^2)+ε,(y^3)+(ε-(x^2)) y+x,
(y^4)+(4/3) (ε-(x^2)) (y^2)+((2 x y)/3)+(1/6) (2 (((ε-(x^2)))^2)-2),
(y^5)+(5/3) (ε-(x^2)) (y^3)+((5 x (y^2))/6)+(1/12) (8 (((ε-(x^2)))^2)-2) y+(1/2) x (ε-(x^2)),(y^6)+2 (ε-(x^2)) (y^4)+x (y^3)+(1/30) (34 (((ε-(x^2)))^2)-6) (y^2)+(13/15) x (ε-(x^2)) y+(1/120) (16 (((ε-(x^2)))^3)-16 (ε-(x^2))+24 (x^2))}

⟹ Quantization condition is not exact for N=7,8,…!

Implementation details

Approximate number of lines of core Mathematica code: 100
Running time (on a 2 GHz computer): ≈(1/3)hour

Many edges and corners: Higher order Brillouin zones of 3D cubic lattices

Problem description

The first Brillouin zone of a lattice Λ with lattice points (Subscript[(OverVector[g]), i]) is the closure of the set of all points (OverVector[x]) such that |(OverVector[x])-(OverVector[0])|≤|(OverVector[x])-(Subscript[(OverVector[g]), i])| for all (Subscript[(OverVector[g]), i])≠(OverVector[0]). (In other words, the first Brillouin zone is the Voronoi cell of the lattice around the origin.) The second Brillouin zone is formed by all points (OverVector[x]) such that |(OverVector[x])-(OverVector[0])|≤|(OverVector[x])-(Subscript[(OverVector[g]), i])| for all (Subscript[(OverVector[g]), i]) not already used in the inequalities of the first Brillouin zone. The third Brillouin zone is formed by all points (OverVector[x]) such that |(OverVector[x])-(OverVector[0])|≤|(OverVector[x])-(Subscript[(OverVector[g]), i])| for all (Subscript[(OverVector[g]), i]) not already used in the inequalities of the first and second Brillouin zone, the fourth is and so on. In other words, the (n+1)th Brillouin zones is the set of points that a line to them crosses exactly n bisector planes. (In computational geometry, a nth-order Brillouin zones would be called an nth-degree Voronoi region or the nth nearest point Voronoi diagram.)

Why is this problem interesting?

• Only the first few Brillouin zones are typically shown in textbooks
• Little is known on how number of faces, surface area, … vary with order

Solution strategy

• Use symmetry
• Calculate bisector planes
• Calculate intersection of bisector planes
• Form bounding polygons
• Assemble bounding polygons to polyhedra

Crucial Mathematica functions and features

• Dynamic caching

• LinearSolve

Results

SC:

BCC:

FCC:

Implementation details

Approximate number of lines of core Mathematica code: 500
Running time (on a 2 GHz computer): ≈(1/2)hour

The end was fate: Long range correlations in Hamlet

Problem description

We start by choosing a (largely arbitrary) binary encoding (1->1,0->-1) for the letters and symbols of the text. We replace all letters and symbols by their encoded version and concatenate all strings. As a result, we get a long sequence (Subscript[s, i]) of -1’s and 1’s. (Typically, the absolute number of 0’s and 1’s will not be random, and on average, -1’s or 1’s will dominate.) Now we look for correlations in the sequence (Subscript[s, i]). The position (Subscript[x, i]) of our random walker is given by (Subscript[x, i])=(∑_(j=1)^i)(Subscript[s, j]), and its walked distance between step k and k+l is d(k,l)=(Subscript[x, k+l])-(Subscript[x, k]). For a true random walk with no correlation between the steps, we expect the mean square fluctuation of (Subscript[((f(l))^2)=⟨((d(k,l))^2)⟩, k])-(Subscript[⟨d(k,l)⟩, k]^2) to have the form f(l)∼(Sqrt[l]). Here, (Subscript[⟨…⟩, k]) denotes the averaging over k. For human-written texts, one finds f(l)∼(l^α) with α>0.5. This indicates long-range correlations are present in the text. (An ideal random walk would have α=0.5.)

Why is this problem interesting?

• Universal for all longer texts and all languages
• Correlations in natural systems

Solution strategy

• Use encoding to convert to a random 1D walk
• Calculate walk explicitly
• Calculate correlations

Crucial Mathematica functions and features

• Packed arrays

• ListConvolve

Results

Hamlet

Resulting dependency between the logarithm of the square root of the mean square fluctuations and the window length.

Mathematica notebook of Chapter 1 of the Programming GuideBooks

Implementation details

Approximate number of lines of core Mathematica code: 100
Running time (on a 2 GHz computer): ≈(1/4)hour

Things don't come together: div(OverVector[B])=0 does not mean closed field lines

Problem description

The magnetic field lines of a current carrying curved wire are defined through:

(OverVector[B])((OverVector[x]))=(∫_(-∞)^∞)(((OverVector[γ])'(s)×((OverVector[x])-(OverVector[γ])(s)))/(|(OverVector[x])-(OverVector[γ])(s)((|)^3)))ds=(∫_(-∞)^∞)(((OverVector[γ])'(t)×((OverVector[x])-(OverVector[γ])(t)))/(|(OverVector[x])-(OverVector[γ])(t)((|)^3))) d t

(OverVector[c])'(t)=(((OverVector[B])((OverVector[c])(t)))/(|(OverVector[B])((OverVector[c])(t))|))

Frequently one reads that from div(OverVector[B])=0 it follows that the field lines close. But from a differential formulation of the density of field lines nothing follows for the geometry in the large of the integral curves.

Why is this problem interesting?

• Under which conditions are field lines closed?
• What are "typical" forms of nonclosed field lines?

Solution strategy

• Calculate magnetic field
• Calculate field lines

Crucial Mathematica functions and features

• NDSolve, NIntegrate

Results

Implementation details

Approximate number of lines of core Mathematica code: 300
Running time (on a 2 GHz computer): ≈(1/2)hour

Fractal tunes: Vibrating square Koch drum

Problem description

The Helmholtz equation for a 2D domain with Dirichlet boundary conditions is:

-Δ u(x,y)=(ω^2) u(x,y)

u(x,y)(Subscript[|, ∂(Subscript[Ω, Koch])])=0

Why is this problem interesting?

• How do the eigenfunctions look for fractal domains?
• Berry corrections to Weyl's law:

n(ω)~(A/(4π)) (ω^2)-{({{c (ω^d), for fractal boundaries}, {(a/(4π))  ω, for smooth boundaries.}})

Solution strategy

• Discretize 2D domain
• Incorporate boundary conditions
• Solve resulting sparse system

Crucial Mathematica functions and features

• Dynamic caching

• D, InterpolatingPolynomial, Eigenvalues

Results

Implementation details

Approximate number of lines of core Mathematica code: 200
Running time (on a 2 GHz computer): ≈3hour

Stairway to heaven: Pendulum in complex time

Problem description

The solution for the elongation of a pendulum are:

φ''(t)=-sin(φ(t)) ∧ φ'(0)=2sin((Subscript[φ, 0])/2)

φ(t)=2 arcsin(sin((Subscript[φ, 0])/2) sn(t|(sin^2)((Subscript[φ, 0])/2)))

How does this solution behave as a function of complex time t?

Why is this problem interesting?

• Nontrivial complex minimum action paths
• Interesting solution curves

Solution strategy

• Find singularities and branch cut
• Carry out correct analytic continuation

Crucial Mathematica functions and features

• NDSolve

Results

Implementation details

Approximate number of lines of core Mathematica code: 100
Running time (on a 2 GHz computer): ≈(1/3)hour

Start to lase: Scattering on a dielectric cylinder

Problem description

The scattering of a vertically polarized electromagnetic wave (Subscript[ℰ, z])(r,φ)(e^(i ω t)) impinging from the left on an infinite dielectric cylinder of radius R and refractive index n can be described by the following equations:

(Subscript[ℰ, z])(r,φ)={({{(Subscript[ℰ, 0])(∑_(m=-∞)^∞)(Subscript[a, m])(Subscript[J, m])(k r)  (e^(i m φ)) , for   r≤R}, {(Subscript[ℰ, 0])(e^(i k x))+(Subscript[ℰ, 0])(∑_(m=-∞)^∞)(Subscript[b, m])(Subscript[H, m]^((1)))(k r)  (e^(i m φ)), for   r≥R}})

(Subscript[a, m])=(i^m) (((Subscript[J, m]^′)(k R)(Subscript[H, 0]^((1)))(k R)-(Subscript[J, m])(k R)(Subscript[H, 0]^((1)′))(k R))/(n  (Subscript[J, m]^′)(n k R)(Subscript[H, 0]^((1)))(k R)-(Subscript[J, m])(n k R)(Subscript[H, 0]^((1)′))(k R)
))
(Subscript[b, m])=(i^m) (((Subscript[J, m]^′)(k R)(Subscript[J, m])(n k R)-(Subscript[J, m])(k R)(Subscript[J, m]^′)(n k R))/(n  (Subscript[J, m]^′)(n k R)(Subscript[H, 0]^((1)))(k R)-(Subscript[J, m])(n k R)(Subscript[H, 0]^((1)′))(k R))).

Here (Subscript[ℰ, 0]) is the amplitude of the incoming wave, k is the wave vector (k=2π/λ), n is the index of refraction, and (Subscript[H, m]^((1)))(z) is again the Hankel function of the first kind.

Why is this problem interesting?

• Lasing conditions for dielectric nanocylinders
• Resonance phenomena
• Wave versions of caustics

Solution strategy

• Efficiently calculate large number of Bessel functions
• Carry out summation

Crucial Mathematica functions and features

• Compilation, dynamic caching

• BesselJ, BesselY, Compile

Results

Animation

Implementation details

Approximate number of lines of core Mathematica code: 100
Running time (on a 2 GHz computer): ≈1hour

Drums that sound alike: Isospectral polygons

Problem description

In 1966 Marc Katz asked the question: Can one hear the shape of a drum?
This question was answered negative in 1991.
Two examples of different–looking and identically-sounding 2D drums are "bilby" and "hawk":

Δψ(x,y)=-(ω^2)ψ(x,y)

Why is this problem interesting?

• Correlation between nodal line length, amplitude autocorrelation function
• Dependence of boundary conditions

Solution strategy

• High-order finite element discretization
• Numerical solution of the resulting sparse eigenvalue system

Crucial Mathematica functions and features

• Compilation, dynamic caching

• Integrate, Eigenvalues

Results

Implementation details

Approximate number of lines of core Mathematica code: 200
Running time (on a 2 GHz computer): ≈2hour

Well, well, well: Hénon–Heiles eigenfunctions

Problem description

Calculate and visualize eigenfunctions in the Hénon–Heiles-like potential:

(-((ℏ^2)/(2M))(Subscript[Δ, x,y])+V(x, y))(Subscript[Ψ, μ])(x, y)=(Subscript[ε, μ])(Subscript[Ψ, μ])(x, y)

V(x,y)=α ((((x^2)+(y^2)))^2)+β ((x^2)+(y^2))+γ ((x^3)-3x (y^2))=α (r^4)+β (r^2)+γ cos(3φ)(r^3).

Why is this problem interesting?

• Nontrivial 2D eigenvalue problem corresponding to a classically chaotic system
• Interesting solution curves

Solution strategy

• Expand in eigenfunctions of the 2D harmonic oscillator with "appropriate" ω
• Diagonalize resulting sparse matrices

Crucial Mathematica functions and features

• NDSolve

Results

Implementation details

Approximate number of lines of core Mathematica code: 200
Running time (on a 2 GHz computer): ≈2hour

Where positive and negative mix: Wigner function of a Duffing oscillator

Problem description

The scattering of a vertically polarized electromagnetic wave (Subscript[ℰ, z])(r,φ)(e^(i ω t)) impinging from the left on an infinite dielectric cylinder of radius R and refractive index n can be described by the following equations:

i ℏ ((∂ψ(x,t))/(∂t))=-((ℏ^2)/(2m))(((∂^2)ψ(x,t))/(∂(t^2)))+(α (x^4)-β (x^2)-γ x cos(ω t))ψ(x,t)

ψ(x,0)=(1/(Sqrt[(Sqrt[2 π]) (Subscript[σ, x])]))exp(-((((x-(Subscript[x, 0])))^2)/(4 (Subscript[σ, x]^2)))+((i (Subscript[p, 0]) x)/ℏ)).

W( x,p, t)=(1/(2π ℏ))(∫_(-∞)^∞)(e^(i (p/ℏ) y))ψ(x-(y/2),t)(OverBar[ψ(x+(y/2),t)])dy

((∂ W( x,p, t))/(∂t))=-(p/m)((∂ W( x,p, t))/(∂x))-(2β x-4 α (x^3)+γ cos(ω t))((∂ W( x,p, t))/(∂p))+(ℏ^3)(α/8) x(((∂^3) W( x,p, t))/(∂(p^3))).

Why is this problem interesting?

• Classical–quantum transition
• Quantitative and qualitative features of quantized chaotic systems
• Appearance of sub-ℏ-structures in the Wigner function

Solution strategy

• Solve time-dependent Schrödinger equation to high precision (6…8 digits)
• Calculate Wigner function through a fractional FFT

Crucial Mathematica functions and features

• Compilation

• NDSolve, Fourier (especially FourierParameters)

Results

Animation

Implementation details

Approximate number of lines of core Mathematica code: 200
Running time (on a 2 GHz computer): ≈6hours

Ten (I think doable) physics problems for Mathematica

1) Develop a regularization procedure for the divergent integrals of the Rayleigh–Schrödinger perturbation theory for V(x)=λ exp((x^4)) in the harmonic oscillator basis

2) Try to find functional equations for p(t;q)=(∏_(k=1)^∞)(1-t (q^k)) for t=exp(i p/q π) for small integer p, q

3) Find the simplest finite energy, propagating, rational function (in x, y, z, and t) solution of the free field Maxwell equations

4) Find proper spatial discretization for the numerical solution of the Schrödinger equation for potentials where wave packets return in finite time from reflection at ∞ (like V(x)=-exp(x))

5) Find a closed knotted field line of a current-carrying trefoil knot

6) Do the nth approximation of the half-space free particle path integral in closed form

7) Calculate the implicit polynomial form of the Poynting vector magnitude for the Lienard-Wiechert potential of a accelerated point charge

8) Investigate the vibration patterns of n-dimensional Hilbert curves (n=2,3,4,5,6,…) and their spectral properties

9) Investigate qualitative and quantitative properties of the geodesics on triple periodic minimal surfaces

10) Numerically calculate number of mutually unbiased bases in small nonprime dimensions (n=6 for a start)