2015-05-31

数式処理システムSageMathの使い方（数体の計算編）

SageMathというフリーかつオープンソースの数式処理システムをご存知でしょうか。

http://www.sagemath.org/ からダウンロードできます（無料！）

とても便利なソフトですので、使い方の簡単な例を紹介しようと思います。

ここで次の計算問題を考えましょう：

体 f:id:hinabit:20150531205416p:plain

において、

f:id:hinabit:20150531210706p:plain

を簡単にせよ。また、x³-2x²+3x-1=0の根としてx=0.7849...+i 1.3071...を選び、

f:id:hinabit:20150531210710p:plain

として埋め込み（体の準同型）

f:id:hinabit:20150531210745p:plain

を定めたとき、 f:id:hinabit:20150531211040p:plain の値はおよそいくらになるか？

SageMathのノートブックという機能を使って上の問題を解いてみることにします。ノートブックとは、グラフィカルかつ対話的に計算を進めていくことができるSageMathの便利機能です。

f:id:hinabit:20150531223708p:plain

さて、SageMathでは有理数体上の既約多項式から数体を構成することができます。

f:id:hinabit:20150531211641p:plain

この状態では、xは f:id:hinabit:20150531205416p:plain の元とみなされます。例えば、x³を評価してみると、次のようになります。

f:id:hinabit:20150531211637p:plain

同じように f:id:hinabit:20150531212210p:plain を評価してみると、

f:id:hinabit:20150531212445p:plain

となります。 f:id:hinabit:20150531212604p:plain であることがわかりました。

問題の後半については、答えを出すだけならいま求めた式にx=0.7849...+i 1.3071...を代入すれば良いのですが、より汎用的な方法として、 NumberField関数のembeddingオプションを使って計算することにします。

f:id:hinabit:20150531213658p:plain

これで f:id:hinabit:20150531211040p:plain ≒-0.316158249 + 1.11729264iがわかりました。

浮動小数点の精度を上げることもできます。

f:id:hinabit:20150531213258p:plain

関数の使い方がわからなくなったら、NumberField?のように打ち込んで評価すると、ドキュメントを表示することができます。

f:id:hinabit:20150531212337p:plain

NumberFieldオブジェクトは様々なメソッドを持っています。一例として、ガロア群の計算をしてみましょう。

f:id:hinabit:20150531232948p:plain

今回はSageMathの、気軽な高機能電卓としての側面のみを紹介しましたが、実際には、SageMathを使えばもっといろいろなことができます。特に、Pythonを使ってプログラミングをすることができるので、Pythonを知っている人なら、新しいプログラミング言語を覚えることなくSageMathでプログラムを書くことができます。

2015-02-15

箱玉系について

授業で「箱玉系」というものを勉強して面白かったので、まとめておきます。
あと、初めてJavaScriptでプログラムを書きました。

「箱玉系」というのはある種のセル・オートマトンです。
無限の長さで一列に並んだ箱に、数字の書いてある玉を入れて、次の規則で時間発展をさせることを考えます。
ただし、玉に書いてある数字は2から9までの整数として、整数kの書いてある玉のことをk-玉と呼ぶことにします。

k=9とおく。
一番左にあるk-玉を箱から取り出して、それより右にある空の箱のうちで一番左にある箱に取り出したk-玉を入れる。
まだ移動していないk-玉のうち、一番左にあるものについて2.と同様の操作を行う。
すべてのk-玉を一度ずつ移動したら、kを1減らして、2.と3.を行う。
k=2になったら操作を終了し、時間を1つ進める。

具体例を下のプログラムで確認しましょう。

以下に挙げる数字の列は、箱玉の初期配置を表します。
ただし、2以上の整数kはk-玉が入っている箱を表し、1は空の箱を表すものとします。

53
443
4431111111111153
653211114431111111111153
6532111111111111443111153

これらの数字の列を「箱玉の初期配置（正の整数の列を入れてください）：」という入力欄にコピペして「input」を押すと、「箱玉の状態：」というウィンドウに数字の付いた箱玉が入ります。この状態で「時間発展」と書いてあるボタンを押してみましょう。

上の初期配置の例の他にも、自分で初期配置を考えていろいろと試してみると面白いと思います。

ここでは先ほどの3つ目の例「4431111111111153」について結果を見てみましょう。

443___________53_____________________________________________________
___443__________53___________________________________________________
______443_________53_________________________________________________
_________443________53_______________________________________________
____________443_______53_____________________________________________
_______________443______53___________________________________________
__________________443_____53_________________________________________
_____________________443____53_______________________________________
________________________443___53_____________________________________
___________________________443__53___________________________________
______________________________443_53_________________________________
_________________________________44_533______________________________
___________________________________44__533___________________________
_____________________________________44___533________________________
_______________________________________44____533_____________________
_________________________________________44_____533__________________
___________________________________________44______533_______________
_____________________________________________44_______533____________
_______________________________________________44________533_________
_________________________________________________44_________533______

他にもいくつかの例を確認すると、次のことがわかると思います：

「玉の番号が非増加な連なり」は時間発展をしても崩れないかたまりになっていて、時間を1つ進めると、連なりの長さだけ右に進む。
「玉の番号が非増加な連なり」同士で衝突が起きると、衝突しているときはかたまりが崩れるが、少し時間が経つと衝突前の連なりの大きさ（上の例だと、3つと2つ）が再び現れる。
衝突前後で連なりの大きさは変わらないが、箱の数字は入れ替わっている。（上の例だと、衝突前は443と53だったが、衝突後は44と533になっている）

実は他にも様々な興味深い現象が起きているのですが、とりあえず今回はこれくらいを挙げるにとどめておきます。

ここでは、最後の現象に注目して、次の問題を考えることにします。

衝突後の数字の状態は、時間発展で確かめることなく、衝突前の数字の情報だけから直接計算できるでしょうか。

例えば上の例だと、「443と53」という情報から、「44と533」を導出したいということです。

上にある時間発展のプログラムを使って他にもいくつか計算すると、
「433と53」→「44と533」
「4432と63」→「42と6433」
「5332と44」→「33と5442」
「6322と442」→「632と4422」
などになっています。
（上のプログラムを使って他にも様々な例でどうなっているかを確認することが出来ます。）

これらの数字の変化を見てみると、一見規則性がないように見えて、「衝突前の数字」→「衝突後の数字」という写像は具体的にはよくわからなさそうに思えます。

ですが、実はこの数字の変化を数学的に捉えられるということがわかっています。
具体的には「クリスタル」と「組み合わせR行列」という組合せ論的な道具を使うことによって、衝突前の数字から衝突後の数字を計算することができます。（[2]のTheorem 4.6）
この「クリスタル」というものは、アフィンLie 代数の量子群の有限次元表現を起源に持つ組み合わせ的な対象です。「組み合わせR行列」は通常のR行列の組合せバージョンと言えます。（もともとのR行列はヤン・バクスター方程式という方程式の解のことであり、これは統計力学や結び目理論などに現れることが知られています。）

単純な規則で時間発展を定義した箱玉系ですが、実はその背景には豊かな数学的な構造があって、非常に面白いなあと私は思いました。
少しでも興味を持った人がいれば、以下の文献をぜひ読んでみればと思います。
この記事は、私がとある講義のゼミで[1]の文献（のほんの一部）を読んだときの内容をまとめたものです。

[1]井上玲, 箱玉系の可積分性 - クリスタルとトロピカル幾何, Summer School 数理物理 2011.
[2]Kaori Fukuda, Masato Okado, Yasuhiko Yamada, Energy Functions in Box Ball Systems, International Journal of Modern Physics A 15 no. 9, 1379-1392 (2000), （arxiv : [math/9908116] Energy Functions in Box Ball Systems）.
[3]国場敦夫, ベーテ仮説と組合せ論, 開かれた数学 5, 朝倉書店(2011), （Amazon.co.jp：ベーテ仮説と組合せ論 (開かれた数学): 国場敦夫: 本）.
[4]Atsushi Nakayashiki, Yasuhiko Yamada, Kostka Polynomials and Energy Functions in Solvable Lattice Models, Selecta Mathematica 3, no. 4 , 547-599, (1997), （arxiv：[q-alg/9512027] Kostka Polynomials and Energy Functions in Solvable Lattice Models）.

2015-02-13

オイラー法の誤差について

期末試験でオイラー法に関する問題が出たんだけど、試験時間内に解けなかったので代わりにここで解くことにします。

オイラー法っていうのは、常微分方程式の近似解を数値的に求める方法で、ルンゲ・クッタ法のもうちょっと簡単バージョンみたいなやつです。具体的には、次のようなものです。

${\lbrack0,1\rbrack\times \mathbb{R}}$ 上の ${C^1}$ 級関数 ${f=f(x,y)}$ に対して、
微分方程式の初期値問題
\begin{align*}
\begin{cases}
y' = f(x,y) &(0\leq x\leq 1),\\
y(0) = y_0
\end{cases}
\end{align*}
が ${\lbrack0,1\rbrack}$ 上 ${C^2}$ 級の解 ${y=y(x)}$ を持つという状況を考えます。
いま、解 ${y(x)}$ の具体的な形はわかってないけど、分点 ${x_0,x_1,\dots , x_n}$ における近似値 ${Y_0,Y_1,\dots ,Y_n}$ を数値的に求めたいなあと思うことにします。
オイラー法と呼ばれる方法では、近似列 ${\{Y_i\}_{i=0}^{n}}$ を次で定めます。

まず、 ${\lbrack 0,1\rbrack}$ を ${n}$ 等分して、 ${x_i = i/n}$ とします。
また、 ${h=1/n}$ として、
\begin{align*}
\begin{cases}
Y_{i+1}= Y_i + hf(x_i , Y_i), \\
Y_0 = y_0
\end{cases}
\end{align*}
で近似列 ${\{Y_i\}_{i=0}^{n}}$ を定めます。

気持ちとしては、 ${y(0)}$ の値と ${y(x)}$ の傾きだけを頼りに、 ${x_{i}}$ での ${y(x)}$ の値を見積もっていこうという算段です。具体的な関数でグラフを書いたらわかったような気になると思います。
あ、Wikipediaに画像があったので貼っておきます。
Euler method - Wikipedia, the free encyclopedia
f:id:hinabit:20150213011059p:plain

それで、このオイラー法に関して期末試験で次のような問題が出題されました。

定数 ${K>0}$ が存在して、任意の ${x\in\lbrack0,1\rbrack}$ および ${y,z\in \mathbb{R}}$ に対し、
\begin{align*}
\lvert f(x,y)-f(x,z)\rvert \leq K\lvert y-z\rvert
\end{align*}
が成り立つとします。（以下、この条件をリプシッツ条件と呼ぶことにします。）
このとき、ある定数 ${B>0}$ （ ${h}$ に無関係）が存在し、
任意の ${h\geq 0}$ に対し、
\begin{align*}
\max_{i=1,\dots,n}\lvert y(x_i)-Y_i \rvert \leq Bh
\end{align*}
が成り立つことを示しなさい。

実際の試験ではもうちょっと誘導があったのですが、「そんなのがなくても解けるぜ！」という硬派な方のため、ここでは省略しました。
以下で証明を与えます。

証明

まず、 ${y(x)}$ に関するテイラーの定理より、各 ${x_i}$ について、
\begin{align*}
y(x_{i+1})=y(x_{i}) + hf(x_i,y(x_i)) + \frac{h^2}{2} y''(\xi)
\end{align*}
を満たす ${\xi \in ( x_{i} , x_{i+1} )}$ が存在することがわかります。

この ${\xi}$ を用いると、
\begin{align*}
y(x_{i+1})-Y_{i+1}=y(x_i) - Y_i + h \{f(x_i,y(x_i)) - f(x_i,Y_i) \} + \frac{h^2}{2}y''(\xi_i)
\end{align*}
となります。つまり、 ${y(x_{i+1})=y(x_{i}) + hf(x_i,y(x_i)) + \frac{h^2}{2} y''(\xi)}$ から ${Y_{i+1}= Y_i + hf(x_i , Y_i)}$ を引いて少し整理しただけです。

グラフで書くとこんな感じです。
f:id:hinabit:20150213014104j:plain

ここで、三角不等式を用いれば、
\begin{align*}
\lvert y(x_{i+1})-Y_{i+1}\rvert\leq\lvert y(x_i) - Y\rvert + h \lvert f(x_i,y(x_i)) - f(x_i,Y_i) \rvert + \frac{h^2}{2}\lvert y''(\xi_i)\rvert
\end{align*}
が成り立ちます。

いま、 ${\lvert y''(x) \rvert \leq M}$ ${(0\leq x \leq 1)}$ となる定数 ${M > 0}$ を取ります。
すると、リプシッツ条件（問題文の2行目に書いてあるやつ）と、ついさっき取った ${M}$ を用いて、
\begin{align*}
\lvert y(x_{i+1})-Y_{i+1}\rvert\leq(1+hK)\lvert y(x_i) - Y\rvert + \frac{h^2 M}{2}
\end{align*}
と評価できます。

これにより、 ${x_{i+1}}$ での誤差評価を ${x_i}$ での誤差評価に落としこむことが出来ました。
これを繰り返し用いることで、
\begin{align*}
\lvert y(x_{i+1})-Y_{i+1}\rvert&\leq(1+hK)\lvert y(x_i) - Y_i\rvert + \frac{h^2 M}{2} \\
&\leq(1+hK)\{(1+hL)\lvert y(x_{i-1}) - Y_{i-1}\rvert + \frac{h^2 M}{2}\}+\frac{h^2 M}{2} \\
&=(1+hK)^{2}\lvert y(x_{i-1}) - Y_{i-1}\rvert +(1+(1+hK)) \frac{h^2 M}{2}\} \\
&\leq(1+hK)^{3}\lvert y(x_{i-2}) - Y_{i-2}\rvert +(1+(1+hK)+(1+hK)^{2}) \frac{h^2 M}{2}\} \\
&\vdots \\
&\leq (1+hK)^{i+1}\lvert y(x_{0}) - Y_{0}\rvert +(1+(1+hK)+\dots (1+hK)^{i}) \frac{h^2 M}{2}\} \\
&=(1+(1+hK)+\dots (1+hK)^{i}) \frac{h^2 M}{2}\} \\
&=\frac{1}{hK}((1+hK)^{i+1} -1) \frac{h^2 M}{2}\} \\
&=\frac{hM}{2K}((1+hK)^{i+1} -1)
\end{align*}
と誤差評価ができます。長いですが、途中に ${\lvert y(x_{0}) - Y_{0}\rvert =0}$ であることや、等比数列の和の公式などを使って計算しています。

ところで、
\begin{align*}
(1+hK) \lt e^{hK}
\end{align*}
ですので、上の式はさらに次のように評価できます。
\begin{align*}
\frac{hM}{2K}((1+hK)^{i+1} -1) &\leq \frac{hM}{2K}(e^{hK(i+1)}-1) \\
&= \frac{hM}{2K}(e^{\frac{i+1}{n}K}-1) \\
&\leq \frac{hM}{2K}(e^{K}-1).
\end{align*}

まとめると、
\begin{align*}
\max_{i=1,\dots,n}\lvert y(x_i)-Y_i \rvert \leq Bh \leq \frac{hM}{2K}(e^{K}-1)
\end{align*}
がわかり、 ${\frac{M}{2K}(e^{K}-1)}$ は ${h}$ に無関係な定数なので、示したかったことが言えました。□

2014-04-30

ブログを始めました

ブログを作ってみました。