Message Boards Message Boards

[WSS19] Solving the wave equation on a torus

In this post, I will share the code and methods I used to solve the wave equation on a torus. To construct a solution to the wave equation on a torus, we need the Laplace-Beltrami operator from differential geometry. For more information on this operator, I recommend looking at this paper: http://www.math.mcgill.ca/toth/spectral%20geometry.pdf. Define coordinates, dimension, surface, and metric:

n = 2; coords = {u, v};
surf = {Cos[u] (2 + Cos[v]), (2 + Cos[v]) Sin[u], Sin[v]};
g = 
FullSimplify[
Table[D[surf, coords[[i]]].D[surf, coords[[j]]], {i, 1, n}, {j, 1, 
n}]];
ParametricPlot3D[surf, {u, 0, 2 \[Pi]}, {v, 0, 2 \[Pi]}]

enter image description here

Find the inverse metric:

ginv = Inverse[g];

Find the determinant of the metric:

detg = Det[g];

For some extra fun, we can find the Christoffel symbols of the second kind, the Riemann curvature tensor, the Ricci curvature tensor, and the Riemann scalar curvature:

\[CapitalGamma] = 
  FullSimplify[ginv.Table[
     1/2 (D[g, coords[[j]]][[k, i]] + 
        D[g, coords[[i]]][[k, j]] - 
        D[g, coords[[k]]][[i, j]]), {k, 1, n}, {i, 1, n}, {j,
       1, n}]];
riemann = 
  FullSimplify[
   Table[Sum[
     D[\[CapitalGamma][[l, i, k]], coords[[j]]] - 
      D[\[CapitalGamma][[l, i, j]], 
       coords[[k]]] + \[CapitalGamma][[l, j, s]] \[CapitalGamma][[s, 
        i, k]] - \[CapitalGamma][[l, k, s]] \[CapitalGamma][[s, i, 
        j]], {s, 1, n}], {l, 1, n}, {i, 1, n}, {j, 1, n}, {k, 1, n}]];
ricci = FullSimplify[
   Sum[Table[
     D[\[CapitalGamma][[l, i, j]], coords[[l]]] - 
      D[\[CapitalGamma][[l, i, l]], 
       coords[[j]]] + \[CapitalGamma][[m, i, j]] \[CapitalGamma][[l, 
        l, m]] - \[CapitalGamma][[m, i, l]] \[CapitalGamma][[l, j, 
        m]], {i, 1, n}, {j, 1, n}], {m, 1, n}, {l, 1, n}]];
scalar = FullSimplify[
   Sum[ginv[[i, j]] ricci[[i, j]], {i, 1, n}, {j, 1, n}]];

However, what really interests us is the Laplace-Beltrami operator.

\[ScriptCapitalL][f_, coord_] := 
  Sum[1/Sqrt[
      det[g] D[ginv[[\[Mu], \[Nu]]] Sqrt[
       detg] D[f, coord[[\[Nu]]]], coord[[\[Mu]]]], {\[Mu], 
    1 n}, {\[Nu], 1, n}];

To solve a differential equation on a torus, we first need to create a gluing diagram to represent how our coordinates u and v map on to the torus:

enter image description here

We can define our quotient space:

quotientspace = Rectangle[{0, 0}, {2 \[Pi], 2 \[Pi]}];

We then need to specify appropriate boundary conditions. I will use the PeriodicBoundaryConditions symbol allowing us to connect the two edges of the gluing diagram. This can also be thought of as creating a large grid of squares where the solution to the partial differential equation repeats itself in each square.

bcs = {PeriodicBoundaryCondition[\[Psi][t, u, v], u == 2 \[Pi], 
   TranslationTransform[{-2 \[Pi], 0}]], 
  PeriodicBoundaryCondition[\[Psi][t, u, v], v == 2 \[Pi], 
   TranslationTransform[{0, -2 \[Pi]}]];

We then should specify some options:

opts = Method -> {"MethodOfLines", "TemporalVariable" -> t, 
    "SpatialDiscretization" -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.001}}};

Let's solve with the initial conditions being a small bump at u = Pi and v = Pi.

sol = NDSolveValue[{\[ScriptCapitalL][ \[Psi][t, u, v], coords] - 
      D[\[Psi][t, u, v], {t, 2}] == 0, 
    bcs, \[Psi][0, u, v] == 
     Exp[-(8 (u - \[Pi]))^2 - (8 (v - \[Pi]))^2], 
    Derivative[1, 0, 0][\[Psi]][0, u, v] == 0}, \[Psi], {t, 0, 
    10}, {u, v} \[Element] quotientspace, opts]; 

To visualize the solution, we will use a color function. Let's define the color function first.

colorf=Function[t, Function[{x,y,z,u,v},ColorData["FuchsiaTones"][Rescale[sol[t,u,v],{-.1/3,.1/3}]]]];

To animate, we need to generate a set of frames:

frames = ParallelTable[
   ParametricPlot3D[surf, {u, 0, 2 \[Pi]}, {v, 0, 2 \[Pi]}, 
    MeshFunctions -> Automatic, ColorFunction -> colorf[t], 
    ColorFunctionScaling -> False, Boxed -> False, Axes -> False, 
    ViewPoint -> {0, 15, 10}], {t, 0, 10, 1/30}];

Putting this all together in combination with some heavy computing power (thanks Boston University SCC), we get an animation that can be found at this link: https://drive.google.com/file/d/1gf0fMxZeMSInLbLemtFkvg1cGDb5ss/view?usp=sharing. Here are some screenshots from that animation: enter image description here enter image description here enter image description here enter image description here enter image description here enter image description here enter image description here

As seen in the animation the PeriodicBoundaryCondition allows the wave to interact with itself. The interaction of the wave with itself creates an increasingly chaotic solution. From my experience adjusting the size of the torus and the size of the initial bump. There is one additional representation I would like to leave you with. We can represent the solution by multiplying the solution by the tangent vector to create a torus that 'wiggles.'

frames = Table[
   ParametricPlot3D[{Cos[
       u] (2 + (1 + 1.2 \[Psi][t, u, v]) Cos[
          v]), (2 + (1 + 1.2 \[Psi][t, u, v]) Cos[v]) Sin[u], 
     Sin[v]}, {u, 0, 2 \[Pi]}, {v, 0, 2 \[Pi]}, 
    ViewPoint -> {0, 15, 10}, PlotRange -> 3.5, Boxed -> False, 
    Axes -> False], {t, 0, 20, .1}];

enter image description here

Attachments:
POSTED BY: Emmy Blumenthal
2 Replies

enter image description here - Congratulations! This post is now featured in our Staff Pick column as distinguished by a badge on your profile of a Featured Contributor! Thank you, keep it coming, and consider contributing your work to the The Notebook Archive!

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract