[Mathematia]微分方程式とグレブナー基底で複素多価関数による曲線の像を求める

複素多価関数

複素関数のうち、$\sqrt{z}$、$\log{z}$、$\arccos{z}$といったものは多価関数です。これらは一つの点に対して２つ以上の点が対応するような関数です。

これらの関数は連続関数なので、連続曲線を連続曲線に写します。一価関数の場合、それ以上言うことはないのですが、
連続な多価関数の場合、面白いことがおきます。それは、「点がある位置から出発して再びもとの位置に戻ってくるような軌道を描くとき、関数により写されたその点の像がもとに戻ってこない場合がある」ということです。

もとに戻ってこない場合

もとに戻ってくる場合

プログラム上の関数として多価関数のこのような振る舞いを自然に再現するのは無理なので、どのプログラム言語も多価関数はとりうる複数の値の内、主値と呼ばれる代表値を返す一価関数として扱われています。しかしその結果、連続曲線の像が不連続点をもつこようになってしまいます。

いかにして多価関数の自然な振る舞いを再現するかというのが本記事の主題です。ただし、ここでは冪根だけを含む多価関数を扱います。

多価関数のとりうるすべての値を表示する

まず冪根を含む関数が取りうるすべての値をリストとして返す関数を作成してみます。

BranchRootsMod[exp_] := Module[
   {p, f, e, k, t},
   If[
    Head[exp] === List,
    t = BranchRootsMod /@ exp;
    ,
    p = Position[exp, (a_^_Rational | a_^-_Rational) /; ! NumericQ[a]];
    If[
     Length[p] > 0,
     f = First[p];
     If[
      Length[f] == 0,
      e = exp;
      k = e /. _^rr_ :> 1/rr;
      t = 
       Table[root[First[e], 1/k]*Exp[I (2 Pi i )/k], {i, 0, 
         Abs@k - 1}];
      ,
      e = Extract[exp, f];
      k = e /. _^rr_ :> 1/rr;
       t = 
       Table[ReplacePart[exp, 
         f -> root[First[e], 1/k]*Exp[I ( 2 Pi i)/k]], {i, 0, 
         Abs@k - 1}];
      ];
     ,
     t = exp;
     ];
    ];
   t
   ];
BranchRoots[exp_] := 
  Module[{e = (FixedPoint[BranchRootsMod, exp] /. root -> Power)}, 
   If[Head[e] === List, Flatten[e], e]];

BranchRootsModは式の中に冪根が含まれていかを探索し、冪根が見つかったら、冪根がとりうる複数の値それぞれに対し、その値を使った元の式の複製を作り、リストとして返します。BranchRootはBranchRootsModを式に繰り返し作用させ、冪根が見つかるたびにネストするリストを作っていきます。最終的に形が変わらなくなった式をFixedPointでとり、Flattenでリストのネストを取り除きます。

この関数を用いると、以下のムービーのように、冪根を含む多価関数がとりうるすべての値を表示することができます。しかし各点に色をつけてみると、色が途中で瞬間的に別の色に変わっているのが分かります。つまり、一価関数にされたことによる多価関数の不自然な振る舞いはまだ残っています。

Manipulate[
 Module[{g1, g2, cl = ColorData[97, "ColorList"], 
   list = (Through[{Re, Im}@#] & /@ (Flatten[{BranchRoots[exp]}]) /. 
      x :> Exp[I \[Theta]]), glist},
  g1 = Labeled[
    Graphics[
     Table[{Blue, PointSize[0.02], Point[list[[k]]]}, {k, 
       Length[list]}],
     PlotRange -> {{-3, 3}, {-3, 3}}, ImageSize -> 500
     ],
    Style[exp, 20], Top
    ];
  g2 = Labeled[
    Graphics[
     Table[{cl[[Mod[k, Length[cl], 1]]], PointSize[0.02], 
       Point[list[[k]]]}, {k, Length[list]}],
     PlotRange -> {{-3, 3}, {-3, 3}}, ImageSize -> 500
     ],
    Style[exp, 20], Top
    ];
  Row[{g1, g2}]

  ]
 ,
 {\[Theta], -12 \[Pi], 12 \[Pi]},
 {exp, (x^(1/3) + x)^(1/3)}
 ]

グレブナー基底と微分方程式を用いたテクニック

この不自然な振る舞いを解決するために、多価関数を多価関数を使わないで表現するという方法を使います。
例えば $w(z)=\sqrt{z}$ だったら、これを変形して $f(w,z)= w^2-z = 0$のような陰関数表示にします。$w(z)=\sqrt{z}$の場合、陰関数表示は簡単に得られましたが、$w(z)=\sqrt{z}+z^{1/3}$のような場合、これは容易いことではありません。それをやってくれるのがグレブナー基底です。

はいできました。グレブナー基底で陰関数表示$f(w,z)$を得たら、元の問題は、$z=z(t)$という軌道を描くときに、$f(w,z)=0$を満たす$w$はどのような軌道を描くかという問題になります。これは微分方程式の問題です。$w$は以下のような微分方程式を満たします。
$$
\frac{dw}{dt}=-\frac{f_z}{f_w}\frac{dz}{dt}
$$
この微分方程式をある初期値$w(t_0)=w_0$に対して解けば求める軌道w(t)が得られます。これらのことをすべてやってくる関数ProjectedTrajectoryを作ります。

ProjectedTrajectory[exp_, var_, traj_, range_] := 
 Module[{w, t, trajd, wd, func, sol},
  func = Simplify[First[GroebnerBasis[{w == exp}, {var, w}]]];
  t = range[[1]];
  trajd = D[traj, t];
  wd = w'[
     t] == ((-(D[func, var]/D[func, w]) trajd) /. {w -> w[t], 
       var -> traj});
  sol = First[
    NDSolve[{wd, 
      w[range[[2]]] == exp /. (var :> (traj /. t -> range[[2]]))}, 
     w, {t, range[[2]], range[[3]]}]];
  w /. sol
  ]

可視化の例：

tmax = 12 \[Pi];
exp = (x^(1/3) + x^(1/2) - 0.5)^(1/3);
traj = ProjectedTrajectory[exp, x, Exp[I t], {t, 0, 2*tmax}];
list = Table[ReIm@traj[t], {t, 0, tmax, tmax/ 256}];

Manipulate[
 Module[{pts},
  pts = ReIm@# & /@ BranchRoots[exp] /. x -> Exp[I t];
  Graphics[
   {
    {Blue, Circle[#, 0.1]} & /@ pts,
    {Orange, PointSize[0.02], Point[ReIm@traj[t]]},
    {Orange, Line@list}
    },
   PlotRange -> {{-3, 3}, {-3, 3}}
   ]],
 {t, 0, tmax}
 ]

ありがとうございました。
この記事は自分がMathematica.stackexchangeで質問したときに、解答が自分が全く思いつかなかったテクニックで感動したので、皆さんにも伝えようと思い書きました。
https://mathematica.stackexchange.com/questions/207792/better-way-to-treat-branch-cut-problem-in-mathematica

すずめ変形コード：

img=Import["すずめ.png"];
Manipulate[
 Module[{gimg, iimg, g, file, \[Theta]},
  \[Theta] = (2 \[Pi])/100 (k - 1);
  gimg = ParametricPlot[
    {u, v}, {u, -Pi, Pi}, {v, -Pi, Pi}, 
    MeshStyle -> {Directive[Thick, Green], Directive[Thick, Orange]}, 
    PlotStyle -> LightYellow, Frame -> None,
    PlotRange -> Pi,
    Mesh -> Full,
    ImageSize -> 1000,
    Epilog -> {
      Rotate[
       {{
         Texture[img],
         Polygon[{{-1, -1}, {1, -1}, {1, 1}, {-1, 1}},

          VertexTextureCoordinates -> {{0, 0}, {1, 0}, {1, 1}, {0, 
             1}}
          ]
         },
        {Blue, PointSize[0.02], Point[{1/Sqrt[2], 1/Sqrt[2]}]}
        },
       \[Theta], {1.5, 1.5}
       ],
      {Blue, Thickness[0.01], 
       Circle[{1.5, 1.5}, 
        Norm[{1.5, 1.5} - {1/Sqrt[2], 1/Sqrt[2]}], {-\[Pi] 3/
           4, -\[Pi] 3 /4 + \[Theta]}]}
      }
    ];
  iimg = Image[
    ImageTransformation[
     Image[gimg],
     {Re@cf[#[[1]] + I #[[2]]], Im@cf[#[[1]] + I #[[2]]]} &,
     DataRange -> {{-Pi, Pi}, {-Pi, Pi}},
     PlotRange -> 0.7*{{-Pi, Pi}, {-Pi, Pi}}
     ], ImageSize -> 1000];
  g =Grid[
     {
      {Style[TraditionalForm[z], 50], Null, 
       Style[TraditionalForm[Sqrt[z]], 50]},
      {
       Image[gimg, ImageSize -> 800],
       Graphics[{Black, Opacity[0.5], Thickness[0.1], Arrowheads[0.4],
          Arrow[{{0, 0}, {1, 0}}]}, 
        PlotRange -> {{0, 1}, {-0.5, 0.5}}, ImageSize -> 100],
       Image[iimg, ImageSize -> 800]
       }
      }
     ]
   ],
 {k, 101}
 ]