Asymptoteでレムニスケートを描く

レムニスケートは$(x^2 +y^2)^2 =a^2 (x^2 -y^2)$で定義される曲線($(\pm a, 0)$は$x$軸との交点を表す)です。例えば$a=1$として$(x^2+y^2)^2=x^2- y^2$を媒介変数表示すると
\[
x(t)=\frac{t(t^2+1)}{t^4+1},\qquad y(t)=-\frac{t(t^2-1)}{t^4+1}
\]
と表せることが知られていますが、媒介変数表示でレムニスケートのグラフを描こうとすると$t$は実数全体を動かす必要があるのでどのように描いたらよいのかわかりませんでした。

Asymptoteで用意されているcontourモジュールを使えば、$x, y$について陽に解かれていない陰関数$f(x, y)=0$の形のグラフを描くことが出来ます。

size(500,0);
import contour;

real a=1; //レムニスケートのx軸との交点座標
real[] c={0};
real f(real x,real y){return (x^2+y^2)^2-(a^2)*(x^2-y^2);};
draw(contour(f,(-1,-1),(1,1),c));

ここでは$f(x, y)=(x^2+y^2)^2 -a^2(x^2-y^2)$とし、$f(x, y)=0$のグラフを描かせています。

 contourでは2変数関数$f(x, y)$と$x$-$y$平面の矩形を指定するための2点(ここでは$(-1,-1)$, $(1,1)$を指定しているので$(x,y)$は$\{(x, y)\in \mathbb{R}^2 \,|\, -1\leq x\leq 1,\quad -1\leq y\leq 1\}$で表される矩形上を動く)、それに実数値の配列を指定する必要があります。配列を指定することで、等高線を描くことが出来るようになっています。

$x$-$y$平面上の曲線が媒介変数表示で$\ell(t)=(x(t),y(t))$の形で表されるならば、曲面$z=g(x, y)$上の曲線で$x$-$y$平面上に射影した像が$\ell(t)$になるものが$\varphi(t)=(x(t),y(t),g(x(t),y(t))$のように簡単に定義出来ます。

しかし$f(x, y)=0$を満たす曲線$\ell$に対して、曲面$z=g(x, y)$上にのっている曲線で$x$-$y$平面に射影すると$\ell$になるようなものの描き方がよくわかりません。

contourやcontour3モジュールでは解決出来なさそうに思えるので、別の方法を考える必要がありそうです。

何故こんなことをしたいかというと、次のような図を完成させたいからです。

この図では曲面上に媒介変数表示のレムニスケートを描いています。ぱっと見はうまくいっているように見えますが、視点を変更して上から見ると、曲線が始点に戻る前に途切れていることがわかります。

これでも十分なように思えますが、曲線の定義域やノードの数を非常に多くしないと形が崩れてしまい、コンパイルにも時間がかかるので、定義式のimplicitな形から描かせたいと思案しています。

Asymptote シルクハット型$C^\infty$級関数

前回$t\leq 0$で$g(t)=0$, $0\leq t\leq 1$で$0\leq g(t)\leq 1$, $t\geq 1$で$g(t)=1$を満たす
$C^\infty$級関数を使いました。今回はこの関数を使って半径$1$の閉円板上で$\varphi(x, y)=1$, 半径2の開円板の外部で$\varphi(x, y)=0$となる関数$\varphi : \mathbb{R}^2 \to \mathbb{R}$の図を書いてみました。
初めて挑戦したことと言えば、関数を条件分岐で定義したくらいです。

シルクハット型というほど立ち上がっていません。次は真横から見た状態です。

かなり綺麗な曲線になっています。$e^{-1/t}$という関数の$t=0$付近での挙動はあまり問題にならないのでしょうか。

size(500,0);

import graph3;
currentprojection=orthographic(2,-1,0.5);

real f(real t){
  if(t>0)
    {return exp(-1/t);}
  else {return 0;}
};
real g(real t){return f(t)/(f(t)+f(1-t));}
real h(real t){return g(t+2)*g(2-t);}

triple F(pair z){return (z.x,z.y,h(abs(z)));};
surface s=surface(F,(-2.5,-2.5),(2.5,2.5),50,50,Spline);
draw(s,lightgray,meshpen=0.01+gray,nolight);

$\sqrt{z}$や$\log z$のRiemann面についての例がこちらで見つけられます。 バネの作図例もあるようでした。さすがにどれも手馴れている感じだなあ。

Asymptoteで$\sqrt{z(z-1)}$のRiemann面を作図する2

引き続き$f(z)=\sqrt{z(z-1)}$のRiemann面について。
今度は双曲線関数の代わりに、多様体論で1の分割を作る際などに用いられる関数を使ってみました。
つまり
\[
f(t)
=
\begin{cases}
e^{-1/t} &(t>0)\\ 0 & (t\leq 0)
\end{cases}
\]
と置き、$g(t)$を
\[
g(t)=\frac{f(t)}{f(t)+f(1-t)}
\]
と定義します。$g(t)$は$g(0)=0$, $g(1)=1$を満たしているので、2つの関数を滑らかにつなぐときに便利です。さらに
\[
\varphi(x)=
\begin{cases}
-\tanh x &(x<0)\\ 0 &(0\leq x\leq 1)\\ \tanh (x-1) &(x>1)
\end{cases}
\]
と定義、$F(x,0)=\varphi(x)$, $F(x,1)=1$を満たす$F(x, y)$を
\[
F(x, y) = (1-g(y))\varphi(x) +g(y)
\]
と定義します。この関数を使い、今回作図したものは次のようになりました。

今回もパーツを寄せ集めて描画しているので、滑らかであって欲しい部分が滑らかでないという問題が残っていますが、前回よりはマシになったのではないでしょうか。曲面をもう少しうまく定義すれば、より滑らかに出来そうです。

Asymptoteで$\sqrt{z}$のRiemann面を作図する。

$f(z)=\sqrt{z}$のRiemann面を作図してみました。
作図しているのは、$(r, \theta)$を独立変数として
\[
\begin{pmatrix}
x(r, \theta)\\ y(r, \theta)\\ z(r, \theta)
\end{pmatrix}
=
\begin{pmatrix}
r\cos \theta \\ r\sin \theta \\ \tfrac{r}{2}\cos \tfrac{\theta}{2}
\end{pmatrix}
\]
という曲面です。$\theta =\pi, 3\pi$が負の虚軸上に対応していて、$\theta =4\pi$で元に戻ります。同様にして、$z(r, \theta)=\tfrac{r}{2}\cos \tfrac{\theta}{n}$などとすれば、$f(z)=\sqrt[n]{z}$に対応するRiemann面が作図出来ます。

size(7cm);
import graph3;
currentprojection=orthographic(2,-3,3);

triple R(pair z) {return ((z.x)*cos(z.y),(z.x)*sin(z.y),(z.x/2)*cos((z.y)/2));};

real a=10;
real b=4;
pen p=black;
pen q=lightgray;
surface s=surface(R,(0,0),(a,b*pi),8,24,Spline);
draw(s,q,meshpen=p);

Asymptoteでバネを描く

Asymptoteで2次元のバネを描いてみます。
まず螺旋$\ell$を$\ell(t)=(a\cos t, a\sin t, t)$と定義し、$\ell$を$y$軸周りに回転したものを$z$軸の正の側から眺める($x$-$y$平面に射影する)と、バネに見えるはずです。
$y$軸周りを$z$軸を$x$軸に向ける向きに角度$\theta$だけ回転する回転行列は
\[
R_y(\theta)
=\begin{pmatrix}\cos \theta & 0 &\sin \theta \\0 & 1 & 0\\-\sin \theta & 0 & \cos \theta\end{pmatrix}
\]
で与えられるので、$\ell(t)$を$R_y(\theta)$で回転させると
\[
\begin{pmatrix}
x \\ y\\ z
\end{pmatrix}
= \begin{pmatrix}\cos \theta & 0 &\sin \theta \\0 & 1 & 0\\-\sin \theta & 0 & \cos \theta\end{pmatrix}
\begin{pmatrix}
 a\cos t \\ b\sin t \\ t
\end{pmatrix}
=
\begin{pmatrix}
a\cos \theta \cos  t + t\sin \theta \\
a\sin t \\ -a\sin\theta \cos t + t \cos \theta
\end{pmatrix}
\]
となります。したがってこの螺旋を$x$-$y$平面に射影したものは
\[
\begin{pmatrix}
x \\ y
\end{pmatrix}
=
\begin{pmatrix}
a\cos \theta \cos  t + t\sin \theta \\
a\sin t
\end{pmatrix}\]
です。これをAsymptoteで描画してみたのが次の図です。

size(10cm);

import graph;
real a=10; // 螺旋の半径
real b=70pi/180; // y軸周りの回転角度
real c=7;

//螺旋 l(t)=(a\cos t, a\sin t, t)をy軸周りに角b回転させたものをx-y平面に射影した図形
pair l(real t) {return (a*(cos(b))*cos(t)+t*sin(b),a*sin(t));};
draw(graph(l,-(2c-1)*pi,(2c)*pi,400,operator ..));

ここでは、400個の点をベジエ曲線で補間したグラフを描いています。200個よりも少ないと歪みが顕著に見えるようになりました。点の数を指定しないと、定義域を広げるほど歪みが激しくなっていきます。
以前にinkscapeでのバネの描き方を紹介しましたが、inkscapeでも数学関数が使えるので同様の方法でバネを描くことが出来ます。バネを伸ばしたり縮めたりするのはinkscapeのほうが楽ですね。

Asymptoteの練習 Log zのRiemann面的な何か

graph3モジュールを使って$\log z$のRiemann面っぽいものを作図してみました。
アイディアはMathematicaのチュートリアルを参考にしました。
配色はColor Scheme Designerを参考に色々試してみましたが、結局グッとくるものができませんでした。Tikz and PGF exampleには配色のセンスのよい図表が沢山あり、あんなふうにしてみたいのですがセンスはどうにもならなさそうです。TikzはTeXをゴリゴリ使える人が使っているという印象があります。以前にTikzを使って花火を打ち上げていた人もいたような。
Asymptoteの用例集も配色にこだわれば使用人口が増えると思うのですが…。
見た目って大切ですよね。

メッシュの数を減らしすぎると意図している曲面にならないようでした。
影を付けるとドス黒い色になったので却下しました。
欲を言えば、Mathematicaのチュートリアルのように矩形に切り取ったものを表示させたいのですが、今のところ方法がわからないのでとりあえず。

size(10cm);
import graph3;
currentprojection=orthographic(3,1,1);

triple R(pair z) {return ((z.x)*cos(z.y),(z.x)*sin(z.y),z.y);};

real a=12;
real b=6;
pen p=rgb("29467F");
pen q=rgb("4575D4");
surface s=surface(R,(0,0),(a,b*pi),6,30,Spline);

draw(s,q,meshpen=p,nolight);

graph3モジュールでは、$x$-$y$平面上の曲面$z=f(x, y)$の他に、$u, v$を独立変数とした$(x, y, z)=(x(u, v), y(u, v), z(u, v))$タイプの曲面を扱うことが出来ます。
今回は後者の方法で$z=(x,y)$を独立２変数として$(x\cos y, x\sin y, y)$を$(x,y)\in [0,12]\times [0,6\pi]$の領域で描画しています。独立変数$(x,y)$はsurfaceを定義する際に定義した、2点$(0,0)$, $(12,6\pi)$を頂点とする矩形上を動きます。

今度は$\sqrt{z}$のリーマン面に挑戦だ。

Asymptoteの練習 無理関数の分岐

無理関数 $f(z)=\sqrt{z(z-1)}$ の分岐(branch cut)についての説明図。
以前にAblowitz & FokasのCOMPLEX VARIABLESという本で読んだ時にわかりやすいと思ったので載せておきます。

$z=r_1 e^{i\theta_1}$, $z-1= r_2e^{i\theta_2}$と置けば、$f(z)=\sqrt{r_1 r_2}\exp(\frac{i}{2}(\theta_1+\theta_2))$と表せる。$\Theta =\frac{1}{2}(\theta_1 + \theta_2)$と置けば、$f(z)=\sqrt{r_1 r_2}e^{i\Theta}$であるが、これを$0\leq \theta_i < 2\pi$ ($i=1,2$)の範囲で考えると、下図のように、矢印方向から$z$の実軸に近づくときに$\Theta$の値が図に書き込んだ値に近づく。

 $\mathrm{Re}\, z>1$や$\mathrm{Re}\, z<0$では実軸の上から近づいても下から近づいても$\Theta$の値が$\bmod 2\pi$で変わらないので、$f(z)$の値は連続であるが、$0<\mathrm{Re}\, z<1$においては実軸の上から近づけた$\Theta$の値と下から近づけた値とが$\bmod 2\pi$で一致しないため、$f(z)$はここで不連続となる。

さらに$2\pi \leq \theta_i < 4\pi$ $(i=1,2)$で同様のことを考えてみると、これら2枚のz平面をつなぎ合わせることで$f(z)$が1価正則となるような面(Riemann面)が構成出来る。

さらに、$-\pi \leq \theta_1 <\pi$, $0\leq \theta_2 <2\pi$とすれば、branch cutは$0$から負の実軸上を通り、$\infty$を通って$1$に至るような形になることが示せるという話。多くの本ではここまで細かく書いていないので、頭の悪い私には有りがたかったです。 さらに$z=u^2$なる変数変換を行なって$g(u)=f(u^2)=\sqrt{u^2(u^2 -1)}$とすれば、$u$平面で考えたときに$z=0$や$z=1$に対応する$u=0$や$u=1$における分岐点が解消されており、例えば$g$は$u=0$を除いた$u=0$の近傍で正則なので、$u=0$を中心とするLaurent展開が可能(puiseux展開)。

size(6cm);
defaultpen(fontsize(10pt));

dot((0,0));
dot((1,0));
label("$0$",(0,0),SW);
label("$1$",(1,0),SE);

real b=0.6;
draw((-b,0)--(1.05+b,0),Arrow(SimpleHead,7));
label("$\mathrm{Re} z$",(1.55,0),S);

pen dashed=linetype("4 4");
real r1 =0.2;
draw(arc((0,0),r1,0,360),dotted,ArcArrow(SimpleHead));
draw(arc((1,0),r1,0,360),dotted,ArcArrow(SimpleHead));

real a=0.6;
draw((0,a)--(0,-a),dashed);
draw((1,a)--(1,-a),dashed);
draw((0,0)--(1,0),1.2+black);
draw((0.5,a/2)--(0.5,0.02),0.3+black,Arrow(SimpleHead,5));
label("$\Theta =\frac{1}{2}\pi$",(0.5,a/2),N);
draw((0.5,-a/2)--(0.5,-0.02),0.3+black,Arrow(SimpleHead,5));
label("$\Theta =\frac{3}{2}\pi$",(0.5,-a/2),S);
draw((-0.35,a/2)--(-0.35,0.02),0.3+black,Arrow(SimpleHead,5));
label("$\Theta =\pi$",(-0.35,a/2),N);
draw((-0.35,-a/2)--(-0.35,-0.02),0.3+black,Arrow(SimpleHead,5));
label("$\Theta =\pi$",(-0.35,-a/2),S);
draw((1.35,a/2)--(1.35,0.02),0.3+black,Arrow(SimpleHead,5));
label("$\Theta =0$",(1.35,a/2),N);
draw((1.35,-a/2)--(1.35,-0.02),0.3+black,Arrow(SimpleHead,5));
label("$\Theta \equiv 0$",(1.35,-a/2),S);
label("$(\mathrm{mod}\,\, 2\pi)$",(1.35,-a/2)+(0.1,-a/2.5));