ヨーキョクデイ

100% pure impurities, which may imply some value. （j は虚数単位）

調和数列の積の部分分数分解とラプラス変換に関する雑多

数学

序（昔話）

さかのぼること高 3 時代、とあるクラス仲間が熱心に部分分数分解していた。ちょっと難しめの数列や積分の問題で部分分数分解は便利ツールとして使用されていたわけが、高々 3 つくらいの積しか登場しなかったはずではある。そこで彼は何を思ったのか、6 つの積バージョンの一例を計算していたのだが、次の式がターゲットだった。

$$\frac{1}{n(n+1)(n+2)(n+3)(n+4)(n+5)}$$

これは興味深いとその結果を持ち帰り、検算してみたところ、彼の計算ミスを発見するわけだが、結果として面白い規則性を発見することになる。

$$\frac{1}{120} \left( \frac{1}{n}-\frac{5}{n+1}+\frac{10}{n+2}-\frac{10}{n+3}+\frac{5}{n+4}-\frac{1}{n+5} \right)$$

なんということだろうか、二項係数の交代和のような何かが出てくるのだ。そして 120 とは何か。そこで、全貌を見るがためにこれの一般化を試みた。すなわち、調和数列の項の積の部分分数分解を試みようという話である。そして証明まで一気にやってしまったのは覚えている。

結果を先に書くと、より一般的な形として次のようになることが証明できた。

$$\prod_{k=0}^{n-1} \frac{1}{a+kd} = \frac{1}{(n-1)! \, d^{n-1}} \sum_{k=0}^{n-1} \frac{(-1)^k \binom{n-1}{k}}{a+kd} \quad (n = 1, 2, 3, \ldots; d \neq 0)$$

調和数列の積の部分分数分解

当時は未熟だったゆえに上記の式は多少 overkill で、これを変形することで一般性を失うことなく簡潔な表記ができることが今回わかった。

公式

非負整数 $n$ と複素数 $z$ について、次の等式が成り立つ。

$$\prod_{k=0}^{n} \frac{1}{z+k} = \frac{1}{n!} \sum_{k=0}^n (-1)^k \binom{n}{k}\frac{1}{z+k}$$

ただし、$z \neq 0, -1, -2, -3, \ldots, -n$。また、$\binom{n}{k}$ は二項係数。

証明

往時は前者のパターンを証明したが、後者の公式を改めて数学的帰納法で証明するが、先にいくつかの道具を用意しておくことにしよう。

道具 1：通分と部分分数分解

$$\frac{1}{x} - \frac{1}{y} = \frac{y-x}{xy}$$ゆえに$$\frac{1}{xy} = \frac{1}{y-x} \qty(\frac{1}{x} - \frac{1}{y} )$$

道具 2：二項係数の性質

$$\begin{align}
\binom{n}{k} &= \frac{n!}{(n-k)! \, k!} \\
&= \frac{n}{n-k} \cdot \frac{(n-1)!}{(n-k-1)! \,k!} \\
&= \frac{n}{n-k} \binom{n-1}{k}
\end{align}$$ゆえに$$\binom{n-1}{k} = \frac{n-k}{n} \binom{n}{k}$$

道具 3：二項定理の利用

二項定理$$(x+y)^n = \sum_{k=0}^n \binom{n}{k} x^{n-k} y^k$$において、$x=1, y=-1$ として、$$(1-1)^n = \sum_{k=0}^n \binom{n}{k} 1^{n-k} (-1)^k$$すなわち$$\sum_{k=0}^n (-1)^k \binom{n}{k} = 0$$

数学的帰納法による証明

まず $n=0$ のとき、左辺について、$$\prod_{k=0}^0 \frac{1}{z+k} = \frac{1}{z}$$

また、右辺について、$$\frac{1}{0!} \sum_{k=0}^0 (-1)^k \binom{0}{k} \frac{1}{z+k} = 1 \cdot (-1)^0 \binom{0}{0} \frac{1}{z} = \frac{1}{z}$$であるから、成り立つ。

次に、$n=m$ のときに成り立つと仮定すると、$$\prod_{k=0}^m \frac{1}{z+k} = \frac{1}{m!} \sum_{k=0}^m (-1)^k \binom{m}{k} \frac{1}{z+k}$$である。この両辺に $\frac{1}{z+m+1}$ を掛けることを考える。

すると、左辺について、$$\qty( \prod_{k=0}^m \frac{1}{z+k} ) \cdot \frac{1}{z+m+1} = \prod_{k=0}^{m+1} \frac{1}{z+k}$$である。右辺について順を追って変形していく。

$$\begin{align}
&\phantom{{}={}} \frac{1}{m!} \qty[ \sum_{k=0}^m (-1)^k \binom{m}{k} \frac{1}{z+k} ] \cdot \frac{1}{z+m+1} \\
&= \frac{1}{m!} \sum_{k=0}^m \qty[ (-1)^k \binom{m}{k} \frac{1}{z+k} \cdot \frac{1}{z+m+1} ] \\
&= \frac{1}{m!} \sum_{k=0}^m \qty[ (-1)^k \binom{m}{k} \frac{1}{m+1-k} \qty( \frac{1}{z+k} - \frac{1}{z+m+1} ) ] \\
&= \frac{1}{m!} \sum_{k=0}^m \qty[ (-1)^k \frac{m+1-k}{m+1} \binom{m+1}{k} \frac{1}{m+1-k} \qty( \frac{1}{z+k} - \frac{1}{z+m+1} ) ] \\
&= \frac{1}{(m+1)!} \sum_{k=0}^m \qty[ (-1)^k \binom{m+1}{k} \qty( \frac{1}{z+k} - \frac{1}{z+m+1} ) ] \\
&= \frac{1}{(m+1)!} \qty[ \sum_{k=0}^m (-1)^k \binom{m+1}{k} \frac{1}{z+k} - \frac{1}{z+m+1} \sum_{k=0}^m (-1)^k \binom{m+1}{k} ] \\
&= \frac{1}{(m+1)!} \qty[ \sum_{k=0}^m (-1)^k \binom{m+1}{k} \frac{1}{z+k} - \frac{1}{z+m+1} \qty( \sum_{k=0}^{m+1} (-1)^k \binom{m+1}{k} - (-1)^{m+1} \binom{m+1}{m+1} ) ] \\
&= \frac{1}{(m+1)!} \qty[ \sum_{k=0}^m (-1)^k \binom{m+1}{k} \frac{1}{z+k} - \frac{1}{z+m+1} \qty(0 - (-1)^{m+1} \binom{m+1}{m+1} ) ] \\
&= \frac{1}{(m+1)!} \qty[ \sum_{k=0}^m (-1)^k \binom{m+1}{k} \frac{1}{z+k} + (-1)^{m+1} \binom{m+1}{m+1} \frac{1}{z+m+1} ] \\
&= \frac{1}{(m+1)!} \sum_{k=0}^{m+1} (-1)^k \binom{m+1}{k} \frac{1}{z+k}
\end{align}$$

よって、$n=m+1$ のときも成立する。ゆえに任意の $n$ について成立することが示された。

ここまでは当時の記述を書き直したものだ。今思えばこれが数学趣味の原点となった数式いじりなのだ。

より簡潔な証明

往時は Heaviside の目隠し法なぞと呼ばれる方法など知らず、それから 2 年ほどしてから Laplace 変換や複素積分の場面で利用することになるのだが、それを使って証明する。

$$\prod_{k=0}^{n} \frac{1}{z+k} = \sum_{k=0}^n \frac{a_k}{z+k}$$と部分分数分解できるから、$k=m$ のときの係数 $a_m$ を求めればよい。この両辺に $z+m$ を掛けて、$z = -m$ とすることを考える。

右辺について、$$
\begin{align}
&\phantom{{}={}} \eval{ (z+m) \sum_{k=0}^n \frac{a_k}{z+k} }_{z=-m} \\
&= \eval{ \qty[ \sum_{k=0}^{m-1} \frac{a_k (z+m)}{z+k} + a_m + \sum_{k=m+1}^{n} \frac{a_k (z+m)}{z+k} ] }_{z=-m} \\
&= \sum_{k=0}^{m-1} \frac{a_k (-m+m)}{-m+k} + a_m + \sum_{k=m+1}^{n} \frac{a_k (-m+m)}{-m+k} \\
&= a_m
\end{align}$$である。左辺について、$$
\begin{align}
&\phantom{{}={}} \eval{ (z+m) \prod_{k=0}^n \frac{1}{z+k} }_{z=-m} \\
&= \eval{ \qty[ \prod_{k=0}^{m-1} \frac{1}{z+k} ] \qty[ \prod_{k=m+1}^{n} \frac{1}{z+k} ] }_{z=-m} \\
&= \qty[ \prod_{k=0}^{m-1} \frac{1}{-m+k} ] \qty[ \prod_{k=m+1}^{n} \frac{1}{-m+k} ] \\
&= \qty[ (-1)^{m} \prod_{k=0}^{m-1} \frac{1}{m-k} ] \qty[ \prod_{k=1}^{n-m} \frac{1}{k} ] \\
&= (-1)^{m} \frac{1}{m!} \cdot \frac{1}{(n-m)!} \\
&= \frac{1}{n!} (-1)^{m} \binom{n}{m}
\end{align}$$であるから、辺々を比較して、$$a_m = \frac{1}{n!} (-1)^{m} \binom{n}{m}$$ゆえ、示された。

補足

2006 阪大（理・前）の大問 3 がまさにこれを別アプローチで証明させる問題らしい。

今回示した公式の左辺は上昇階乗冪の逆数として書けて、$$\frac{1}{z^\overline{n+1}} = \frac{1}{n!} \sum_{k=0}^n (-1)^k \binom{n}{k} \frac{1}{z+k}$$である。

グーグル先生によると、背景として次のものが存在するらしい。

二項変換

次のようにして数列 $\{a_n\}$ から数列 $\{s_n\}$ をでっち上げることを二項変換と呼ぶらしい。

$$s_n = \sum_{k=0}^n {(-1)}^k \binom{n}{k} a_k$$

Melzak の公式

非負整数 $ n$ と高々 $ n$ 次の多項式 $f(x)$ について、次の等式が成り立つらしい。

$$\sum_{k=0}^n {(-1)}^k \binom{n}{k} \frac{f(x-k)}{y+k} = n!\,f(x+y) \prod_{k=0}^n \frac{1}{y+k}$$

ただし $ y \neq 0,-1,-2,\ldots,-n$。これを Melzak の公式あるいは Melzak の恒等式と呼ぶらしい。二項変換の例である。

arxiv.org

今回証明したものは $f(x) = 1$ の例であって、上記によるとこの等式は well-known らしい。

arxiv.org

調和数列の積の部分分数分解と Laplace 変換

数列といえば z 変換かもしれないが。ここではとにかく形式的に式をいじることにする。

Laplace 変換の定義

$f(t)$ の（片側）Laplace 変換 $F(s) = \mathcal{L} [f]$ を次のように定義する。

$$F(s) = \mathcal{L} [f] = \int_0^{\infty} f(t) e^{-st} \dd{t}$$

指数減衰

ここでは唐突ながらも $s$ の関数として$$F(s) = \frac{1}{s+\alpha} $$を考える。これを逆 Laplace 変換した $f(t)$ は、$$f(t) = \mathcal{L}^{-1} \qty[ F(s) ] = e^{-\alpha t}$$となる。単位ステップ関数を掛けるべきかもしれないが、面倒なので省略した。

部分分数分解の公式

上記を踏まえて、調和数列の積の部分分数分解の公式の逆変換を考えると、$$
\begin{align}
&\phantom{{}={}} \mathcal{L}^{-1} \qty[ \sum_{k=0}^n (-1)^k \binom{n}{k} \frac{1}{s+k\alpha} ] \\
&= \sum_{k=0}^n (-1)^k \binom{n}{k} \mathcal{L}^{-1} \qty[ \frac{1}{s+k\alpha} ] \\
&= \sum_{k=0}^n (-1)^k \binom{n}{k} e^{-k\alpha t} \\
&= \sum_{k=0}^n \binom{n}{k} (-e^{-\alpha t})^k \\
&= (1 - e^{-\alpha t})^n
\end{align}$$となる。これにより、次のように Laplace 変換できるだろう。

$$\mathcal{L} \qty[ (1 - e^{-\alpha t})^n ] = \sum_{k=0}^n (-1)^k \binom{n}{k} \frac{1}{s+k\alpha}$$

序盤に書いた煩雑なバージョン、$$\prod_{k=0}^{n-1} \frac{1}{a+kd} = \frac{1}{(n-1)! \, d^{n-1}} \sum_{k=0}^{n-1} \frac{(-1)^k \binom{n-1}{k}}{a+kd}$$を書き直して、$$\prod_{k=0}^n \frac{1}{a+kd} = \frac{1}{n! \, d^n} \sum_{k=0}^n (-1)^k \binom{n}{k} \frac{1}{a+kd}$$となるので、$d$ を $\alpha$ に、$a$ を $s$ に置き直せば、$$\mathcal{L} \qty[ (1 - e^{-\alpha t})^n ] = n! \, \alpha^n \prod_{k=0}^n \frac{1}{s+k\alpha}$$という具合で積の形で書けることがわかる。特に $n=1$ のときバージョンの
$$\mathcal{L} \qty[ 1 - e^{-\alpha t} ] = \frac{1}{s}-\frac{1}{s+\alpha} = \frac{\alpha}{s(s+\alpha)}$$は、力学的には空気抵抗がある場合の落下速度や、電気回路的には RC 直列回路におけるコンデンサの電圧についての過渡現象なんかを表す微分方程式を、Laplace 変換を使って解こうというときあたりにモロに出てくる形であるわけではある。

そもそも Laplace 変換表を眺めていたら今回の部分分数分解の公式を思い出して、この記事を書くことになったのだが、この節はおまけという名の実は本編なのだ。

$n$ 乗バージョンに何か物理的意味はあるのだろうか。立ち上がりがなまってくる感じだと思うが。確率論方面に何かありそうな感じではある。また、連続なパラメータであるほうが面白そうではある。

Borel 変換と総括的な応用

Borel 総和法というのがあるらしい。

Borel 変換

次のような形式的冪級数 $f(z)$ を考える。

$$f(z) = \sum_{k=0}^\infty a_k z^k$$

$f(z)$ の Borel 変換 $\mathcal{B} [f] (t)$ を、次のように定義する。

$$\mathcal{B} [f] (t) = \sum_{k=0}^\infty a_k \frac{t^k}{k!}$$

Borel 変換 と Laplace 変換

冪乗を階乗で割ったやつ

ところで、唐突ながらも Laplace 変換について、$$\mathcal{L}^{-1} \qty[ \frac{1}{s^{n+1}} ](t) = \frac{t^n}{n!}$$であったから、$$\mathcal{B} [f] (t) = \sum_{k=0}^\infty a_k \mathcal{L}^{-1} \qty[ \frac{1}{s^{n+1}} ](t)$$である。

Borel 変換の再構築

上記の結果から、$f(z)$ について、$z= 1/s$ として両辺に $1/s$ を掛ける。

$$\begin{align} \frac{1}{s} f \qty(\frac{1}{s} )
&= \frac{1}{s} \sum_{k=0}^\infty a_k \frac{1}{s^n} \\
&= \sum_{k=0}^\infty a_k \frac{1}{s^{n+1}}
\end{align}$$

さらに、両辺を（形式的に）逆 Laplace 変換し、和と積分を交換する。

$$\begin{align} \mathcal{L}^{-1} \qty[ \frac{1}{s} f \qty(\frac{1}{s} ) ]
&= \mathcal{L}^{-1} \qty[ \sum_{k=0}^\infty a_k \frac{1}{s^{n+1}} ] \\
&= \sum_{k=0}^\infty a_k \mathcal{L}^{-1} \qty[ \frac{1}{s^{n+1}} ] \\
&= \mathcal{B} [f]
\end{align}$$となるから、Borel 変換の正体が見えた。

Borel 和

上記の冪級数 $f(z)$ とその Borel 変換 $\mathcal{B} [f] (t)$ に対して、$f(z)$ の Borel 和 $S[f](z)$を、$$S[f](z) = \frac{1}{z} \mathcal{L} [\mathcal{B} [f](t)] \qty(\frac{1}{z} )$$で定義する。

逆 Borel 変換

性質のよい $f(z)$ の Borel 和 $S[f](z)$ は再び $f(z)$ に戻り、すなわち、$$S[f](z) = f(z) = \frac{1}{z} \mathcal{L} [\mathcal{B} [f](t)] \qty(\frac{1}{z} )$$である。この場合これは $\mathcal{B} [f] (t)$ を $f(z)$ に戻す変換となっている。ここで、$t$ の関数 $f^*(t)$ について、逆 Borel 変換 $\mathcal{B}^{-1} [f^*](z)$ を次のように定義する。

$$\mathcal{B}^{-1} [f^*](z) = \frac{1}{z} \mathcal{L} [f^*] \qty(\frac{1}{z} )$$

こんなところでも Laplace 変換をフル活用する。

母関数と Borel 変換

数列 $ \{a_n\}$ について、通常型母関数 $ g(x)$ を次のように定義する。

$$g(x) = \sum_{k=0}^\infty a_k x^k$$

同様に、指数型母関数 $ g^*(x)$ を次のように定義する。

$$g^*(x) = \sum_{k=0}^\infty a_k \frac{x^k}{k!}$$

先ほどの Borel 変換の話はここに響いてきて、つまり、通常型母関数を Borel 変換することで指数型母関数を得られるといううれしさがあるのだ。

第 2 種 Stirling 数の母関数

$S(n,k)$ で第 2 種 Stirling 数を表すとする。ここで、$k=m$ と固定したときの母関数を考える。

指数型母関数

指数型母関数 $g^*(x)$ を考える。

$$g^*(x) = \sum_{n=0}^\infty S(n,m) \frac{x^n}{n!}$$

収束半径は無限大であり、閉じた式としては次のようになる。（Ref. 上記 26.8.12）

$$g^*(x) = \frac{(e^x-1)^m}{m!}$$

指数型母関数から通常型母関数へ

通常型母関数 $g(x)$ は、$$g(x) = \sum_{n=0}^\infty S(n,m) x^n$$で表され、この Borel 変換によって指数型母関数が得られているわけで、つまり、$$\mathcal{B}[g(z)] = \frac{(e^t-1)^m}{m!} $$という関係が得られたということだ。この右辺についての逆 Borel 変換を考えよう。少しさかのぼって、$$\mathcal{L} [ (1 - e^{-\alpha t})^n ] = n! \, \alpha^n \prod_{k=0}^n \frac{1}{s+k\alpha}$$であったことを思い出すと、$\alpha = -1$ として、$$\mathcal{L} [ (1 - e^t)^n ] = n! \, (-1)^n \prod_{k=0}^n \frac{1}{s-k}$$であり、すなわち、$$\mathcal{L} \qty[ \frac{(e^t - 1)^n}{n!} ] = \prod_{k=0}^n \frac{1}{s-k}$$である。すると、Laplace 変換を経由して所望の逆 Borel 変換、すなわち $g(z)$ の Borel 和 $S[g](z)$ を機械的に算出できて、$$
\begin{align} S[g](z)
&= \mathcal{B}^{-1} \qty[ \frac{(e^t-1)^m}{m!} ] \\
&= \frac{1}{z} \eval{ \mathcal{L} \qty[ \frac{(e^t-1)^m}{m!} ] }_{s=\frac{1}{z}} \\
&= \frac{1}{z} \prod_{k=0}^m \frac{1}{\frac{1}{z}-k} \\
&= \frac{1}{z} \prod_{k=0}^m \frac{z}{1-kz} \\
&= \prod_{k=1}^m \frac{z}{1-kz}
\end{align}
$$となり、これは通常型母関数の閉じた形だ。すなわち閉じた形のまま指数型母関数を通常型母関数に写すことができた。まさに総決算である。

通常型母関数

通常型母関数 $g(x)$ を考える。

$$g(x) = \sum_{n=0}^\infty S(n,m) x^n$$

収束半径は $1/m$ であり、その下で閉じた式は次のようになるという。（Ref. 上記 26.8.11）

$$\begin{align} g(x)
&= \frac{x^m}{(1-x)⁢(1-2⁢x)⁢ \cdots ⁢(1-m⁢x)} \\
&= \prod_{k=1}^m \frac{x}{1-kx}
\end{align}$$

雑感

等差数列の項の積は Stirling 数と関連があるわけで、その性質を眺めていていたら母関数に見覚えのある形が現れていて、そこから深掘りする中でこの変換の存在を知り、この応用例を書いたわけだが。

〆

遠くなってしまった昔に見つけたものを再考し、さらに Laplace 変換というスイスアーミーナイフでけちょんけちょんにぶん殴っては切り刻むことにより、新たな知見が得られて、それもう脳内麻薬がドバドバなのだ。

数学的な厳密性は考えずに、形式的な式変形を行ったのだが、ヘヴィサイドさんの名言を言い訳として使わせてほしい。今回のテーマにちなんで。

ja.wikipedia.org

楕円の中心に対する極座標表示と楕円の扇形の面積と

数学

日常を数式を通して見るシリーズ。今回は円と楕円を行ったり来たりしながら、楕円の面積に関する性質を見ていく。シリーズの前回は 観覧車のゴンドラをぐわんぐわん揺らしたい - ヨーキョクデイ ということにしておく。

円と楕円と線型変換

まず、$a >0, b >0$ として、次の方程式で表される楕円を考える。

$$\frac{x^2}{a^2}+\frac{y^2}{b^2} = 1$$

この楕円は原点を中心とする単位円を x 軸方向に $a$ 倍、y 軸方向に $b$ 倍に拡大したものである。ここで、この単位円上の点 $(x_C, y_C)$ を考える。単位円を変形して楕円を作ったとき、この点が対応する楕円上の点を $(x_E, y_E)$ とする。すると、これらの点には次のような関係がある。

$$\begin{align}
\mqty(x_E \\ y_E)
&= \mqty(a x_C \\ b y_C) \\
&= \mqty(a & 0 \\ 0 & b) \mqty(x_C \\ y_C)
\end{align}$$

つまり、$A = \smqty(a & 0 \\ 0 & b)$ なる行列 $A$ で表される線型変換を考えていたことになる。

楕円の面積

それゆえ、この変換で写される図形の面積は元の面積の $\abs{\det A}$ 倍、つまり $ab$ 倍になる。実際に、単位円の面積は $\pi$ であるが、この楕円の面積は $\pi ab$ である。

中心角

単位円と楕円
再び、単位円上の点 $(x_C, y_C)$ を考える。円の中心とこの点を結ぶ線分（半径）が x 軸正の向きとのなす角を $\theta_C \in [0, 2\pi]$ とすると、$$\mqty(x_C \\ y_C) = \mqty(\cos \theta_C \\ \sin \theta_C)$$であるから、この点が $A$ による変換によって写された楕円上の点 $(x_E, y_E)$ は $\theta_C$ を用いて次のよう表せる。

$$\begin{align}
\mqty(x_E \\ y_E)
&= A \mqty(x_C \\ y_C) \\
&= \mqty(a \cos \theta_C \\ b \sin \theta_C)
\end{align}$$

一方で、楕円の中心と点 $(x_E, y_E)$ を結ぶ線分 $r_E$（半径のような何か、正式名称はわからない）が x 軸正の向きとのなす角を $\theta_E \in [0, 2\pi]$ とすると、$$
\begin{align}
\mqty(x_E \\ y_E )
&= \mqty(r_E \cos \theta_E \\ r_E \sin \theta_E) \\
&= r_E \mqty(\cos \theta_E \\ \sin \theta_E )
\end{align}$$となるから、$$\cos \theta_E : \sin \theta_E = a \cos \theta_C : b \sin \theta_C$$あるいは、これを変形して、$$\cos \theta_C : \sin \theta_C = b \cos \theta_E : a \sin \theta_E$$という関係が導かれる。分母が $0$ になりうるので $\tan$ で書かなかった。変換前後で点の象限は変化しないので、普通の逆正接 $\arctan$ を用いずに、$\mathrm{arctan2} (y, x) \in [0, 2\pi)$ を使うことにして、$\theta_C, \theta_E \in [0, 2\pi)$ について、$$\begin{align}
\theta_E(\theta_C) &= \mathrm{arctan2} (b \sin \theta_C, a \cos \theta_C) \\
\theta_C(\theta_E) &= \mathrm{arctan2} (a \sin \theta_E, b \cos \theta_E)
\end{align}$$とでき、また、自然に$\theta_E(\theta_C = 2\pi) = 2\pi$、$\theta_C(\theta_E = 2\pi) = 2\pi$ と対応づけることにする。これによって x 軸に対しての中心角についての相互関係を得る。

en.wikipedia.org

半径と極方程式

さて、単位円の半径は一定して $1$ であるが、先ほどの楕円の半径のような何か $r_E$ を考えると、この長さは楕円上の点の座標に依存する。$r_E=\sqrt{{x_E}^2 + {y_E}^2}$ であるが、$(x_E, y_E) = (a \cos \theta_C, b \sin \theta_C)$ であったから、$\theta_C$ を用いて $$r_E = \sqrt{ a^2 \cos^2 \theta_C + b^2 \sin^2 \theta_C}$$ と書ける。また、$(x_E, y_E) = (r_E \cos \theta_E, r_E \sin \theta_E)$ とも書けたから、これを当初の楕円の方程式に代入して変形することで、$\theta_E$ を用いた式が楽に得られる。

$$r_E = \frac{ab}{\sqrt{ b^2 \cos^2 \theta_E + a^2 \sin^2 \theta_E}}$$

これは $0 \leq \theta_E \leq 2 \pi$ において楕円を表す極方程式となっている。

また、任意の区間 $[\theta_{E\,begin}, \theta_{E\,end}] \subset [0, 2\pi]$ については x 軸に対しての中心角が $\theta_{E\,begin}$ から $\theta_{E\,end}$ であるところの楕円弧を表す。

極方程式と扇形の面積

一般の極座標による関数 $r(\theta)$ について、極方程式 $r = r(\theta), \theta \in [\theta_{begin}, \theta_{end}]$ で表される連続な曲線を考える。ここで、この曲線と、同じく極方程式で表される 2 直線 $\theta = \theta_{begin}$、$ \theta = \theta_{end}$ とで囲まれる領域 $S$ の面積は次の積分で与えられる。

$$\begin{align}
S(\theta_{begin}, \theta_{end})
&= \iint_S \dd{S} \\
&= \int_{\theta_{begin}}^{\theta_{end}} \int_{0}^{r(\theta)} r \dd{r} \dd{\theta} \\
&= \int_{\theta_{begin}}^{\theta_{end}} \frac{1}{2} r(\theta)^2 \dd{\theta}
\end{align}$$

ここで、次を考える。

$$S_E(\theta_E) = \int_{0}^{\theta_E} \frac{1}{2} r_E(\theta)^2 \dd{\theta}$$

前述の公式より、$S_E(\theta_E)$ なる関数は x 軸に対しての中心角 $\theta_E$ であるときの楕円の扇形の面積を表す。ところで、この積分を愚直に計算する必要はない。楕円が行列 $A$ によって単位円を変換したものであり、扇形に関してもそうであるため、対応する部分の単位円の扇形の面積がわかれば、これを $\abs{\det A} = ab$ 倍すれば、楕円の扇型の面積もわかるのだった。伏線回収。改めて、（x 軸に対しての）中心角 $\theta_C$ であるときの単位円の扇形の面積を $S_C(\theta_C)$ とすれば、$$S_C(\theta_C) = \frac{\theta_C}{2}$$である。さらに、 $\theta_E \in [0, 2\pi)$ に対応するものとして $\theta_C(\theta_E) = \mathrm{arctan2} (a \sin \theta_E, b \cos \theta_E)$ を考えていたから、$$
\begin{align}
S_E(\theta_E)
&= ab S_C( \mathrm{arctan2} (a \sin \theta_E, b \cos \theta_E) ) \\
&= \frac{ab}{2} \mathrm{arctan2} (a \sin \theta_E, b \cos \theta_E)
\end{align}$$

また、$\theta_E = 2\pi$ については $\theta_C = 2\pi$ が対応していたから、$$
\begin{align}
S_E(2\pi)
&= ab S_C(2\pi) \\
&= \pi ab
\end{align}$$となるがゆえ、当初に求めた楕円の面積が導出できた。

扇形の面積、再考

2 角で挟まれる楕円の扇形
ここまで考えてきた楕円の扇形は x 軸とある中心角 $\theta_E$ で構成されるものであった。ここでは任意の $\theta_{E\,begin}$ と $\theta_{E\,end}$ ($ 0 \leq \theta_{E\,begin} < \theta_{E\,end} \leq 2\pi$) で挟まれる扇形を考えたい。この面積を $S_E(\theta_{E\,begin}, \theta_{E\,end})$ とすると、先ほどの公式から次のように表せる。

$$\begin{align}
S_E(\theta_{E\,begin}, \theta_{E\,end})
&= \int_{\theta_{E\,begin}}^{\theta_{E\,end}} \frac{1}{2} r_E(\theta)^2 \dd{\theta} \\
&= \int_{0}^{\theta_{E\,end}} \frac{1}{2} r_E(\theta)^2 \dd{\theta} - \int_{0}^{\theta_{E\,begin}} \frac{1}{2} r_E(\theta)^2 \dd{\theta} \\
&= S_E(\theta_{E\,end}) - S_E(\theta_{E\,begin})
\end{align}$$

まあ、積分表記はなくても明らかだが。

例題：楕円形のピザを等分割する

実はこれが本題だった。楕円形のピザを同じ面積のピースに切り分けたい事情があったときに役に立つ話だった。そんなピザがあるかは別として。とにかく、ピザの中心から放射状に切り分けるとして、これを $N (\geq 2)$ 個に分割するにはどういう角度で切っていけばいいのかというレシピを構成するための前置きが長くなったので、そのレシピというかアルゴリズムというかを書く。

  1. 楕円の長軸と短軸を割り出す
  2. そのどちらかを x 軸上に、もう一方を y 軸上に置き、中心を定める
  3. 長軸と短軸の長さから適当な $a$ と $b$ を算出する
  4. $N$ 回行われるのカットのうちのいずれか 1 つにおける、x 軸正の向きに対する角度 $\theta_{E\,0} \in [0, 2\pi)$ を適当に決める（この角度が基準となる）
  5. $\theta_{E\,0}$ より、$\theta_{C\,0} = \theta_{C}(\theta_{E\,0})$ を求める
  6. $k \in \{1, 2, \ldots, N-1\}$ なるすべての $k$ に対して $\theta_{C\,k} = \theta_{C\,0} + 2\pi k/N$ を求める
    1. $\theta_{C\,k} > 2\pi$ の場合は改めて $\theta_{C\,k} = \theta_{C\,k} - 2\pi$ とする
  7. それらの $\theta_{C\,k}$ について、$\theta_{E\,k} = \theta_{E}(\theta_{C\,k})$ を求める
  8. $\{\theta_{E\,0}, \theta_{E\,1}, \ldots, \theta_{E\,N-1}\}$ なる $N$ 個の角度こそが切るべき x 軸正の向きに対する角度である

下に $\frac{x^2}{5^2}+\frac{y^2}{3^2} = 1$ についての例を図示する。

基準 π/4 : 5 分割
基準 0 : 8 分割

今回のための SageMath のコードを晒しておく。

def arctan2_nonneg(y, x):  # ∈ [0, 2*pi)
    t = arctan2(y, x)  # ∈ (-pi, pi]
    return t if t >= 0 else t + 2 * pi


assert arctan2_nonneg(0, 1) == 0
assert arctan2_nonneg(1, 0) == 1 / 2 * pi
assert arctan2_nonneg(0, -1) == pi
assert arctan2_nonneg(-1, -1) == 5 / 4 * pi

a = 5
b = 3


def _normalize(angle):  # ∈ [0, 2*pi]
    if angle < 0:
        return _normalize(angle + 2 * pi)
    elif 2 * pi < angle:
        return _normalize(angle - 2 * pi)
    else:
        return angle


def translate_central_angle_from_c_to_e(t):  # ∈ [0, 2*pi]
    t = _normalize(t)
    assert 0 <= t <= 2 * pi
    return 2 * pi if t == 2 * pi else arctan2_nonneg(b * sin(t), a * cos(t))


def translate_central_angle_from_e_to_c(t):  # ∈ [0, 2*pi]
    t = _normalize(t)
    assert 0 <= t <= 2 * pi
    return 2 * pi if t == 2 * pi else arctan2_nonneg(a * sin(t), b * cos(t))


# (x/a)**2 + (y/b)**2 = 1
def elliptical_sector(central_angle_rel_to_x_axis
                      ):  # central_angle_rel_to_x_axis ∈ [0, 2*pi]
    return a * b * unit_circular_sector(
        translate_central_angle_from_e_to_c(central_angle_rel_to_x_axis))


# x**2 + y**2 = 1
def unit_circular_sector(central_angle):  # central_angle ∈ [0, 2*pi]
    return _normalize(central_angle) / 2


n = 5
t_e_0 = pi / 4

t_c_0 = translate_central_angle_from_e_to_c(t_e_0)
t_e = [
    translate_central_angle_from_c_to_e(t_c_0 + k * 2 * pi / n)
    for k in range(n)
]
r_by_t_e = [a * b / sqrt((b * cos(t))**2 + (a * sin(t))**2) for t in t_e]
seg = [
    line([(0, 0), (r * cos(t), r * sin(t))],
         linestyle="--",
         rgbcolor=(1, 0, 0),
         thickness=2) for (r, t) in zip(r_by_t_e, t_e)
]

E = ellipse((0, 0), a, b)
G = E + reduce(lambda x, y: x + y, seg)
G.show(aspect_ratio=1)

埼玉遠征

日常

現実逃避を目論んで、21 日から 23 日にかけて埼玉県に遊びに行っていたが。主たる目的は川越の観光。観光向けメインストリートと思しきところは人が多すぎて人混み嫌いの俺には向かないかなぁ、と。レトロな感じの建物が並ぶ町並みとか撮りたかったけど、人がね。朝とかならいいのかもしれんが。神社巡りできたのはよし。

野球好きマンなので、念願叶って條辺のうどん屋に行けたのは非常に満足。

観覧車のゴンドラをぐわんぐわん揺らしたい

物理学 数学

さて、観覧車のゴンドラは振り子のようになっていて、理論上は観覧車の回転によってゴンドラが揺れるはずである。この運動を表す運動方程式を求めたい。

いや、これはモチベーションとしては後付けであって、二重振り子に興味を持ち、これは非常にカオスなのでもうちょっと秩序のある振り子の運動を調べたいというきっかけから、支点を周期的に動かせばいいのではということになり、等速円運動くらいがいいよね、という過程のもとであり、これはすなわち観覧車とゴンドラのモデル化だよね、という話なのだ。

それはさておき、さっそく運動方程式を作る。ラグランジュの方程式だ。解析力学。

鉛直面があり、そこに適当な原点を中心とする適当な半径 $r$ の円がある。位置ベクトルを $\vec r$ とする点が円周上を動くとする。その点を支点として長さ $l$ の理想的な棒がぶら下がっており、その先には質量 $m$ の質点がある。この質点の鉛直面上での運動を考える。例によって摩擦や空気抵抗はないものとする。

支点から質点へのベクトルを $\vec l$ とする。すると、質点の位置ベクトル $\vec p$ は$$\vec p = \vec r + \vec l$$で表される。$\vec p$ だけど運動量ではない。

ここで質点の運動エネルギー $T$ を求めたい。縦軸を鉛直下向きを正の向きとして、縦軸と $\vec r$ のなす角を $\phi$、縦軸と $\vec l$ のなす角を $\theta$ とする。これらは時刻 $t$ の関数であり、特に一定の角速度 $\omega$ を用いて $\phi=\omega t$ となる等速円運動を考える。また、$\theta$ こそが最終的に求めたい関数であり、一般化座標として用いるものである。

$$\begin{align}
T &= \frac{m}{2} \dot{\vec p}^2 \\
&= \frac{m}{2} \qty( \dot{\vec r} + \dot{\vec l} )^2 \\
&= \frac{m}{2} \qty( \dot{\vec r}^2 + 2 \dot{\vec r} \cdot \dot{\vec l} + \dot{\vec l}^2 ) \\
&= m \qty( \frac{1}{2}r^2 \dot{\phi}^2 + rl\dot{\phi}\dot{\theta} \cos( \phi-\theta ) + \frac{1}{2}l^2\dot{\theta}^2 )
\end{align}$$

こうだ。$\dot{\vec l}$ は $\vec l$ と垂直だが、$\dot \theta$ の正負によって向きが変わることに注意して内積をとりたいところ。

さらに、ポテンシャル $U$ を求めたい。$\phi=0$ かつ $\theta=0$ のときに最小となるように適当に定める。

$$\begin{align}
U &=-mg \qty( \vec p の鉛直下向き成分 ) \\
&=-mg \qty( r \cos \phi + l \cos \theta )
\end{align}$$

さて、これらを用いてラグランジアン $L=T-U$ を求める。

$$L =m \qty( \frac{1}{2}r^2 \dot{\phi}^2 + rl\dot{\phi}\dot{\theta} \cos( \phi-\theta ) + \frac{1}{2}l^2\dot{\theta}^2 ) + mg \qty(r \cos \phi + l \cos \theta )$$

あとはラグランジアンを微分したりしたやつを求める作業。まずは $\dot \theta$ で偏微分。

$$\pdv{L}{\dot{\theta}} = m \qty( rl\dot{\phi} \cos( \phi-\theta ) + l^2 \dot{\theta} ) $$

これの時間微分。

$$\begin{align}
\dv{t} \qty( \pdv{L}{\dot{\theta}} )
&= m \qty[ rl \qty( \ddot{\phi} \cos( \phi-\theta ) - \dot{\phi}\qty( \dot{\phi}-\dot{\theta} ) \sin( \phi-\theta )) +l^2\ddot{\theta} ] \\
&= m \qty( rl \qty( \dot{\phi}\dot{\theta}-\dot{\phi}^2 ) \sin( \phi-\theta ) +l^2\ddot{\theta} )
\end{align}$$

$\phi = \omega t$ であったから、$\ddot \phi = \dot \omega = 0$ である。別に、$\theta$ で偏微分。

$$\pdv{L}{\theta}= m \qty( rl\dot{\phi}\dot{\theta} \sin( \phi-\theta ) - gl \sin \theta )$$

はい。ここからラグランジュの運動方程式 $\dv{t}\qty( \pdv{L}{\dot{\theta}} )-\pdv{L}{\theta}=0$ を作る。

$$\begin{align}
\dv{t}\qty( \pdv{L}{\dot{\theta}} )-\pdv{L}{\theta} &= m \qty( rl \qty( \dot{\phi}\dot{\theta}-\dot{\phi}^2 ) \sin( \phi-\theta ) +l^2\ddot{\theta} - rl\dot{\phi}\dot{\theta} \sin( \phi-\theta ) + gl \sin \theta ) \\
&= m \qty( -rl \dot{\phi}^2 \sin( \phi-\theta ) +l^2\ddot{\theta} + gl \sin \theta ) \\
&= m \qty( -rl \omega^2 \sin(\omega t-\theta ) +l^2\ddot{\theta} + gl \sin \theta ) = 0
\end{align}$$

よって、次の運動方程式を得る。

$$\ddot{\theta} = \frac{r \omega^2}{l} \sin(\omega t-\theta ) - \frac{g}{l} \sin \theta $$

普通の単振り子の強制振動の式に似た形。ぐわんぐわん揺らしたいので、微小振動どころの話ではなく、$\theta \approx 0$ とした近似は使えない。解析的には解けなそう。

解析力学 (物理入門コース 新装版)

解析力学 (物理入門コース 新装版)

Amazon
ならば数値計算だ。SageMath を使ってルンゲ・クッタ法による数値解を出してもらう。ただ、$\theta (t)$ の概形を見てもあまり面白くないので、条件を変えながら実際の振り子の動きを見てみようではないか。とりあえず $\omega$ と $\theta(0)$ のバリエーションで $\vec l$ の動き方の図を出した。各種定数や初期値を各図の左上に表示してある。w は $\omega$ のつもり。

気持ちよくぐわんぐわん揺れてくれたが、実際の乗客は気持ち悪くなること必至。

円の半径よりも棒が適度に長いほうが面白いかな。

今回使用した SageMath のコードを書いておく。

def f(r,l,w,g=9.8,t0=0,theta0=0,theta_d0=0):
    (t,theta,theta_d)=var('t theta theta_d')

    A=desolve_system_rk4([theta_d,r*w*w*sin(w*t-theta)/l-g*sin(theta)/l],
    [theta,theta_d],
    ics=[t0,theta0,theta_d0],
    ivar=t,
    end_points=20)

    coord_lim = (r+l)*1.3
    legend1='r = %.2f\nl = %.2f\nw = %.2f\ng = %.2f'%(r,l,w,g)
    legend2='t0 = %.2f\ntheta0 = %.2f\ntheta_d0 = %.2f'%(t0,theta0,theta_d0)

    def coord_r(t):
        return r*sin(w*t_),r*cos(w*t_)

    def coord_p(t,th):
        return vector(coord_r(t))+vector((l*sin(th),l*cos(th)))

    an=animate([
    text(legend1+'\n'+legend2,(0,1),axis_coords=True,vertical_alignment='top',horizontal_alignment='left')+
    text('t = %.1f\nphi= %.2f\ntheta = %.2f\n'%(t_,w*t_,theta_),(1,1),axis_coords=True,vertical_alignment='top',horizontal_alignment='right')+
    circle((0,0),r,linestyle=':')+
    arrow(coord_r(t_),coord_p(t_,theta_),color='red') for t_,theta_,_ in A],
    xmin=-coord_lim,xmax=coord_lim,ymin=coord_lim,ymax=-coord_lim,aspect_ratio=1,figsize=10,dpi=70)
    an.gif(delay=10,savefile=(legend1+legend2).replace(' = ','_').replace('\n',''),show_path=True)

根はアウトドア派の近況

日常 はてな

ガキのころは外で遊んだり、スポーツに明け暮れたりしていたものだが、やはりおっさんにもなると外出がおっくうになったりもして、休みの日は部屋でひねもすのたりしているのがよかったりもして。

そんな昨今、もやっとした心境の中、外に出ては美しい風景が見たい気分になっているのであって、さらにそれを写真に収めたいわけである。日本海はすぐそこにあるし、鳥海山だって軽いドライブ感覚で行けるというロケーションを活用しないわけにはいかないし、そのすばらしさを再認識すべきであるのだ。




こんな感じで、ちょっと休眠していたカメラを本格稼働させつつある日々。写真楽しい。

https://www.satofull.jp

＃私のふるさと

特別お題「私のふるさと」キャンペーン
Sponsored by さとふる

東京遠征

日常 野球

ヤクルトファンの後輩氏に連れられる形で神宮へ。セの試合は初。そして初の外野席。なかなかひどい試合ではあった。

こちらを 0 時に出て、3 時ころから久々の再会に話を盛り上がらせつつ後輩氏の車で夜通し東京に向かい、昼から都内を散策した後に 17 時開始という形だったので、観戦中にこっくりしていたけど。とりあえず丸 2 日間寝ずに活動できることはわかった。運転好きの後輩氏とはいえ長距離運転をお任せしたのはさることながら、炎天下において自分の趣味や嗜好に付き合ってもらったりしてありがたい限り。よい心の洗濯ができた。

無線ルータ買った

日常 PC

先週のことですけどね。BUFFALO の WSR-2533DHP-CB というやつ。今使っているのが 10 年近く前に買ったやつで、軽く驚き。自室の無線環境の向上を目指して。部屋でスマホをいじるときに、不安定感抜群だったので。11ac に進化したのだが、どうかな。

BUFFALO 無線LAN親機 11ac/n/a/g/b 1733+800Mbps Giga ブラック【 Nintendo Switch動作確認済】 WSR-2533DHP-CB

BUFFALO 無線LAN親機 11ac/n/a/g/b 1733+800Mbps Giga ブラック【 Nintendo Switch動作確認済】 WSR-2533DHP-CB