This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

Make an animation showing how an initially localized perturbation evolves on a torus under the wave equation

(((∂^2)ψ(t;(Subscript[φ, 1]),(Subscript[φ, 2])))/(∂(t^2)))=Δ ψ(t;(Subscript[φ, 1]),(Subscript[φ, 2]))
Δ=(1/(Sqrt[g((Subscript[φ, 1]),(Subscript[φ, 2]))]))(∑_(i,j=1)^2)(∂/(∂(Subscript[φ, i])))((Sqrt[g((Subscript[φ, 1]),(Subscript[φ, 2]))])(g^(i j))((Subscript[φ, 1]),(Subscript[φ, 2]))(∂/(∂(Subscript[φ, j])))).

Here (Subscript[g, i j])((Subscript[φ, 1]),(Subscript[φ, 2])) is the metric tensor and g((Subscript[φ, 1]),(Subscript[φ, 2])) its determinant.

Use the torus parametrization x=((Subscript[r, 1])+(Subscript[r, 2])cos((Subscript[φ, 2])))cos((Subscript[φ, 1])),  y=((Subscript[r, 1])+(Subscript[r, 2])cos((Subscript[φ, 2])))sin((Subscript[φ, 1])),  z=(Subscript[r, 2])sin((Subscript[φ, 2])).

(* calculate g_ij *)
𝓍𝓎𝓏 = {(r1 + r2 Cos[φ[2]]) Cos[φ[1]], (r1 + r2 Cos[φ[2]]) Sin[φ[1]], r2 Sin[φ[2]]};
{𝓍𝓎𝓏φ1, 𝓍𝓎𝓏φ2} = D[𝓍𝓎𝓏, #]& /@ {φ[1], φ[2]};
gll = {{𝓍𝓎𝓏φ1.𝓍𝓎𝓏φ1, 𝓍𝓎𝓏φ1.𝓍𝓎𝓏φ2}, {𝓍𝓎𝓏φ1.𝓍𝓎𝓏φ2, 𝓍𝓎𝓏φ2.𝓍𝓎𝓏φ2}}

The line element has a very simple form.

g = Det[gll] // Simplify

And here is the Laplace operator in (Subscript[φ, 1]), (Subscript[φ, 2])-coordinates [✶], [✶].

(* calculate g^ij *)
guu = Factor //@ Together //@ Inverse[gll];
Sum[1/Sqrt[g] D[Sqrt[g] guu[[i, j]] D[Ψ[φ[1], φ[2]], φ[i]], φ[j]], 
    {i, 2}, {j, 2}] // FullSimplify

We will, of course, solve the wave equation numerically for the torus. In the following, we choose the parameters (Subscript[r, 1])=2 and (Subscript[r, 2])=1.

r1 = 1; r2 = 1/2;

For an animation, we have to calculate many frames. To carry out the numerical PDE solution in an memory-efficient manner, we will use the functions NDSolve`ProcessEquations and NDSolve`Iterate. The following input sets up the problem. We use an initially localized, smooth perturbation. Because of the periodic boundary conditions in (Subscript[φ, 1]), (Subscript[φ, 2]), we use the "Pseudospectral" method setting for the method of lines.

{state} = 
NDSolve`ProcessEquations[{D[u[t, φ1, φ2], t, t] == 
   1/r2^2 D[u[t, φ1, φ2], φ2, φ2] - 
   Sin[φ2]/(r1 + r2 Cos[φ2]) D[u[t, φ1, φ2], φ2] +
   1/(r1 + r2 Cos[φ2])^2 D[u[t, φ1, φ2], φ1, φ1],
  (* periodic boundary conditions in φ1, φ2 *)
  u[t, -Pi, φ2] == u[t, +Pi, φ2],
  Derivative[0, 1, 0][u][t, -Pi, φ2] == Derivative[0, 1, 0][u][t, +Pi, φ2],
  u[t, φ1, -Pi] == u[t, φ1, +Pi],
  Derivative[0, 0, 1][u][t, φ1, -Pi] == Derivative[0, 0, 1][u][t, φ1, +Pi],
  (* initial perturbation *)
  u[0, φ1, φ2] == If[Sqrt[φ1^2 + φ2^2] < 1/3 Pi, 
                 (Sign[#] #^2)&[Cos[3/2 Pi Sqrt[φ1^2 + φ2^2]/(1/3 Pi)]], 0], 
  Derivative[1, 0, 0][u][0, φ1, φ2] == 0},
  {}, t, {φ1, -Pi, Pi}, {φ2, -Pi, Pi},
  (* option settings for the numerical solution *)
  MaxSteps -> 10^5, PrecisionGoal -> 2, AccuracyGoal -> 2, 
  Method -> {"MethodOfLines", "SpatialDiscretization" -> 
              {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral", 
               "MaxPoints" -> {120, 60}, "MinPoints" -> {120, 60}}}]

To visualize the propagation of the perturbation on the torus surface, we implement a function coloredTorusPlot3D. Unperturbed regions ψ((Subscript[φ, 1]),(Subscript[φ, 2]))≈0 are shown in green, positive perturbations ψ((Subscript[φ, 1]),(Subscript[φ, 2]))>0 in red, and negative perturbations ψ((Subscript[φ, 1]),(Subscript[φ, 2]))<0 in blue. In addition, we represent the size of the perturbations as deformations of the torus in its local normal direction.

coloredTorusPlot3D[sols_, {ppφ1_, ppφ2_}, opts___] := 
Module[{α = 1/4, torusPoint, torusNormal, data, points, toColor, polys},
 (* point and normal on the torus *)
 torusPoint[{φ1_, φ2_}, r1_, r2_] = {r1 Cos[φ1] + r2 Cos[φ1] Cos[φ2],
                        r1 Sin[φ1] + r2 Sin[φ1] Cos[φ2], r2 Sin[φ2]};
 torusNormal[{φ1_, φ2_}, r1_, r2_] =  
             {Cos[φ1] Cos[φ2], Cos[φ2] Sin[φ1], Sin[φ2]};
 (* scaled perturbations *)
 data = Table[Evaluate[{{φ1, φ2}, ArcTan[25 sols[[1, 2]]]/(Pi/2)}], 
              {φ1, -Pi, Pi, 2Pi/ppφ1}, {φ2, -Pi, Pi, 2Pi/ppφ2}]; 
 (* add perturbations in normal direction *)
 points = Apply[(torusPoint[#1, r1, r2] + 
                α r2 #2 torusNormal[#1, r1, r2])&, data, {2}];
 (* color according to perturbations *)
 toColor[x_] := SurfaceColor[#, #, 2.2]&[
                RGBColor[UnitStep[x] x, 1 - Abs[x], -UnitStep[-x] x]];
 (* polygons *)
 polys = Table[{toColor[(data[[i, j, 2]] + data[[i + 1, j, 2]] + 
                         data[[i + 1, j + 1, 2]] + data[[i, j + 1, 2]])/4], 
                Polygon[{points[[i, j]], points[[i + 1, j]], 
                         points[[i + 1,j + 1]], points[[i, j + 1]]}]}, 
               {i, Length[data] - 1}, {j, Length[data[[1]]] - 1}];
 (* graphics of the torus with perturbations *)
 Graphics3D[{EdgeForm[], polys}, opts, Boxed -> False, Axes -> False,
            PlotRange -> {{-2, 2}, {-2, 2}, {-3/4, 3/4}},
            (* view toward initial perturbation *) 
            ViewPoint -> {3, 0, 2.4}]]

Now, it is straightforward to propagate the solution forward in time and to visualize the resulting wave propagation. The function waveEquationTorusAnimation generates the actual animation.

waveEquationTorusAnimation[ΔT_, frames_, DoTable_, showIt_, {ppφ1_, ppφ2_}] := 
Module[{uCurrent},
DoTable[(* visualize solution *) 
        (* solve wave equation in time steps *)
        NDSolve`Iterate[state, j ΔT];
        (* form current solution *)
        uCurrent = NDSolve`ProcessSolutions[state, "Forward"];
        showIt @ coloredTorusPlot3D[uCurrent, {ppφ1, ppφ2}],
        {j, 0, frames}]]

And here is the resulting animation. Due to the superposition of the wave contributions along different rays, we see complicated time-dependent patterns emerge. Interestingly, in the beginning the initial pulse splits into four single pulse, two traveling along the “small” circumference and two traveling along the “large” circumference. From the 2D nature of the torus surface one might expect a approximately radial propagation of the initial disturbance. So, in some sense, the wave propagation on this torus is more in (ℝ^1)×(ℝ^1) than in (ℝ^2).

Show[GraphicsArray[#]]& /@ 
     Partition[waveEquationTorusAnimation[2/3, 8, Table, Identity, {240, 120}], 3]

waveEquationTorusAnimation[1/10, 90, Do, Show, {240, 120}]