東京工業大学 ロボット技術研究会

└(՞ةڼ◔)」＜やっほ〜したろうだよォ〜！

最近色んな人からドヤ顔白衣と呼ばれるようになりました。
アノ人の影響力はすごいですね…。

i-engineer.jp

ついに今日、卒業論文発表が終わりました！
これで学部生としてやるべきことは尽くしました。今はただ天命を待つだけです。
というわけで、新年度に向けてしたろうが生まれ変わりました！→@shitaro2016
心機一転、ゼロからの再スタート！これからの活躍に期待してください！

概要

従来では、FortranでOpenGLを使うにはGLUTという選択肢しかなかった。
しかし、モダンなFortranではC言語とより密な連携がとれるようになったことに着目したしたろうは、
自分でGLFWとGLEWのラッパーを書くことで、モダンなOpenGLをモダンなFortranで使用することを実現したのだ！

続きを読む

これはrogy Advent Calendar 2015の17日目の記事です.

こんにちは、この度初めてrogyブログの記事を書かせて頂く事になりました、rogy15のくろすけ(@AstroBoy2030)と申します。
幸福の物理に配属されてますが宗教色は(今のところ)ないです。

はじめに

自分はrogy希少種である理学部の人間なので今回の記事もどちらかと言うと理学系のお話です。
しかし、理学系と言えば工学部やその他の一般の方からしたら象牙の塔の住人の巣窟と思われがちです(個人的な偏見)。
なので工学部の方にも興味を少しでも持っていただけるような話題としてFPU問題を取り上げることにしました。

Fermi-Pasta-Ulamの問題

本題のFermi-Pasta-Ulamの問題ですが、ざっくり言うとFPU問題とは「何個かある格子点を弦のように繋いで両端を固定させたときに隣り合う格子点同士の引っ張り合い方を複雑化(非線形化)させて、ある1つのモードにエネルギーを与えて動かすとモード毎のエネルギーがどのように分配されていくか」といったある種のエネルギー分配問題です。
ここでいう「モード」は弦の基準振動、要するに山が何個かあるsin波の形のシンプルな振動の事を指します。
この問題は現代的コンピュータが登場した1950年代に研究されていたもので、当時のロスアラモス研究所が保有していた高性能計算機MANIACというコンピュータで数値解析されていたものです。
60年前は研究所で保有する大型コンピュータで行われていた実験が今や個人で持てるノートPCで出来るようになったわけです。
さてFPU問題の数式モデルを作ってみます。
まずは普通の弦を下の図のように$N+1$個の格子点の集合体に捉えたとして、波動方程式を作るときと同じ要領で格子点$P_{i}$の変位$x_{i}$($i=1,2,...,N-1$)に関する運動方程式を立ててみます。
なお$P_{0}$と$P_{N}$は端点なので固定されていて運動方程式を考える必要はありません。


このとき簡単のため、格子点の質量と復元力の係数を両方$1$としています。
\[\frac{{\rm d}^{2} x_{i}}{{\rm d} t^{2}} = (x_{i+1}-x_{i})-(x_{i}-x_{i-1})\]
次に運動方程式に非線形の項を挿入します。
ここでは非線形項として隣合う格子点同士の相互作用として$1$次の項に続く$2$次の項を採用し、その係数を$\alpha$とします。
\[\frac{{\rm d}^{2} x_{i}}{{\rm d} t^{2}} = (x_{i+1}-x_{i})+(x_{i}-x_{i-1})+\alpha\left[(x_{i+1}-x_{i})^{2}-(x_{i}-x_{i-1})^{2}\right]\]
これが今回取り扱うFPUの運動方程式となります。
因みに元の論文で取り扱われている非線形項はこの$2$次の項に限らず、$3$次の項であったり折れ線グラフで表される複雑な$1$次関数だったりします。
要するに非線形項であれば良いのです。

Runge-Kutta法

まずは微分方程式を数値計算で解く方法として(多分)メジャーなものであるRunge-Kutta法という手法を使います。
Runge-Kutta法は一階の微分方程式を解く方法なので先程の運動方程式を格子点の速度$v_{i}$を使って下の式のように分解します。
\[ \frac{{\rm d} v_{i}}{{\rm d} t} = \left( x_{i+1} - x_{i} \right) - \left( x_{i} - x_{i-1} \right) + \alpha \left[ \left( x_{i+1}-x_{i}\right)^{2} - \left( x_{i}-x_{i-1} \right)^{2} \right] \\
   \frac{{\rm d} x_{i}}{{\rm d} t} = v_{i}\]
この2つの運動方程式を$i=1,2,..,N-1$に関して解けば格子点の位置が求まります。
Runge-Kutta法の手法は何種類かありますが効率よく計算できる4次の手法が使われることが多いので、ここでも4次の方法で計算します。
以下、細かい原理は省いた全体的な計算の流れです(沿え字$i$は省略)。

1. 時刻$t_{0}$における値$v_{0}$、$x_{0}$を用いて$k_{(1)}=\displaystyle\frac{{\rm d} v}{{\rm d} t}{(v_{0},x_{0},t_{0})}$、$j_{(1)}=\displaystyle\frac{{\rm d} x}{{\rm d} t}{(v_{0},x_{0},t_{0})}$を計算する。
2. $k_{(2)}=\displaystyle\frac{{\rm d} v}{{\rm d} t}{(v_{0}+k_{(1)}\cdot\frac{dt}{2},x_{0}+j_{(1)}\cdot\frac{dt}{2},t_{0}+\frac{dt}{2})}$、$j_{(2)}=\displaystyle\frac{{\rm d} x}{{\rm d} t}{(v_{0}+k_{(1)}\cdot\frac{dt}{2},x_{0}+j_{(1)}\cdot\frac{dt}{2},t_{0}+\frac{dt}{2})}$を計算
3. $k_{(3)}=\displaystyle\frac{{\rm d} v}{{\rm d} t}{(v_{0}+k_{(2)}\cdot\frac{dt}{2},x_{0}+j_{(2)}\cdot\frac{dt}{2},t_{0}+\frac{dt}{2})}$、$j_{(3)}=\displaystyle\frac{{\rm d} x}{{\rm d} t}{(v_{0}+k_{(2)}\cdot\frac{dt}{2},x_{0}+j_{(2)}\cdot\frac{dt}{2},t_{0}+\frac{dt}{2})}$を計算
4. $k_{(4)}=\displaystyle\frac{{\rm d} v}{{\rm d} t}{(v_{0}+k_{(3)}\cdot dt,x_{0}+j_{(3)}\cdot dt,t_{0}+dt)}$、$j_{(4)}=\displaystyle\frac{{\rm d} x}{{\rm d} t}{(v_{0}+k_{(3)}\cdot dt,x_{0}+j_{(3)}\cdot dt,t_{0}+dt)}$を計算
5. 時間dt後の$x$、$v$の値$x_{dt}$、$v_{dt}$を$v_{dt}=v_{0}+\displaystyle\frac{dt}{6}\cdot(k_{1}+2k_{2}+2k_{3}+k_{4})$、$x_{dt}=x_{0}+\displaystyle\frac{dt}{6}\cdot(j_{1}+2j_{2}+2j_{3}+j_{4})$と計算する
6. $(v_{dt},x_{dt})\to(v_{0},x_{0})$として手順1.に戻る

以上の操作で微分方程式の解が数値的に解かれます。
イメージとしては、次の時刻に進むとき現在時刻の周りで再計算して誤差を小さくしている感じです。

Runge-Kutta法の欠点

ある程度精度良く微分方程式の解を計算出来るRunge-Kutta法ですが、FPUのシミュレーションをする上でのある欠点を抱えています。
幾ら細かく計算したとしても、連続的な関数を離散的に計算している以上ある程度の誤差は生じます。
ですが、この誤差がシミュレーションを実行していく中で蓄積されていき、ついには計算結果が発散してしまう事が起きてしまうのです。
これはエネルギー保存則がある物理のシミュレーションでは非常にまずい問題であり、実際にFPUシミュレータを作っていて簡単にエネルギーが発散してしまいシミュレーションにならない事も多々ありました(というか長めにシミュレーションしてると大体発散しました)。

Symplectic積分

Runge-Kutta法だとエネルギーが発散し易いので別の手法を使います。
ここで先程立てたFPUのモデルに関してHamiltonianという新しい関数を導入します。
このHamiltonian(以下$H$)は運動量$p_{i}$と座標$q_{i}$の関数($i=0,1,...,N$)として表された運動エネルギー$T(p,q)$とポテンシャルエネルギー$V(p,q)$で
\[H=T+V\]
と定義されます。

※注意※

このHamitonianの定義は一般的なものではありません。
ここではあたかも$H$が全エネルギー$E=T+V$そのものに見えてますが、必ずしもHamiltonianが全エネルギーEを表したものとは限らないのです。
ですが、ここではHamiltonianが全エネルギーに一致する場合で話を進めます。

さて、FPUのモデルをHamiltonianで書き直してみます。
ある格子$P_{i}$に対応する運動量と座標をそれぞれ$p_{i}$、$q_{i}$とします。
ここではHamiltonianが全エネルギーの式と一致するとしているので、格子毎の運動エネルギーと弾性エネルギーを足し合わされば$H$を計算できます。
ここでは簡単のため格子点の質量$m$とバネ係数を両方$1$としてあります。
\[H=\sum_{i=0}^{N}\frac{p_{i}^{2}}{2}+\sum_{i=1}^{N}\frac{1}{2}(q_{i}-q_{i-1})^{2}+\sum_{i=1}^{N}\frac{\alpha}{3}(q_{i}-q_{i-1})^{3}\]
さて、Hamlitonianを導入したところで次に下の式で表される正準方程式なる式を使います。
\[\frac{\partial H}{\partial p_{i}}=\frac{{\rm d} q_{i}}{{\rm d} t} \\
  \frac{\partial H}{\partial q_{i}}=-\frac{{\rm d} p_{i}}{{\rm d} t}\]
この2式はHamiltonianを用いて表された運動方程式で、この方程式の解が格子点の運動を決定します。
そこで正準方程式を使って格子点の運動を数値的に解析する事を考えます。
単純に正準方程式の両辺を時間$t$で積分すると$p$、$q$の数値解を得ることが出来ますが、あまり精度が良くないのでより効率的な手法を用います。
先程の式をベクトルの形で書き直します($D$：Hamiltonianに依存する演算子)。
\[ D \left ( \begin{array}{ccc} p_{i} \\ q_{i} \end{array} \right ) = \displaystyle\frac{\rm d}{{\rm d} t}\left ( \begin{array}{ccc} p_{i} \\ q_{i} \end{array} \right ) \]
さて、上の式の通りベクトル$\left (\begin{array}{ccc} p_{i}\\q_{i}\end{array}\right)$に$D$を作用させると$\left (\begin{array}{ccc} p_{i}\\q_{i}\end{array}\right)$を時間微分したものが得られるので、$D$を$n$回だけ作用させると$n$回時間微分したものが得られるわけです。
ここで$\left (\begin{array}{ccc} p_{(\Delta t)}\\q_{(\Delta t)}\end{array}\right)$を$t=0$周りでTaylor展開すると
\[
\left (\begin{array}{ccc} p_{(\Delta t)}\\q_{(\Delta t)}\end{array}\right)=
\left (\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+
\displaystyle\frac{1}{1!}\frac{\rm d}{{\rm d} t}\Delta t\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+\cdots+
\displaystyle\frac{1}{m!}\frac{{\rm d}^{m}}{{\rm d}^{m} t}(\Delta t)^{m}\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+\cdots
\]
ここで$D^{m} \left ( \begin{array}{ccc} p \\ q \end{array} \right ) = \displaystyle\frac{{\rm d}^{m}}{{\rm d}^{m} t}\left ( \begin{array}{ccc} p \\ q \end{array} \right )$であるので$D^{0}={\rm i.d.}$と約束すると
\[
\left (\begin{array}{ccc} p_{(\Delta t)}\\q_{(\Delta t)}\end{array}\right)=
\frac{1}{0!}D^{0}(\Delta t)^{0}\left (\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+
\displaystyle\frac{1}{1!}\left(\Delta t D\right)^{1}\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+\cdots+
\displaystyle\frac{1}{m!}\left(\Delta t D\right)^{m}\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)+\cdots
\]
と表すことができるので、結果的に
\[
\left (\begin{array}{ccc} p_{(\Delta t)}\\q_{(\Delta t)}\end{array}\right)={\rm Exp}\left(\Delta t D\right)\left(\begin{array}{ccc} p_{(0)}\\q_{(0)}\end{array}\right)\\
\left({\rm Exp}\left(\Delta t D\right)=\frac{1}{0!}(\Delta t D)^{0}+\frac{1}{1!}(\Delta t D)^{1}+\cdots+\frac{1}{m!}(\Delta t D)^{m}\cdots\right)
\]
と書けます。
上記の${\rm Exp}(\Delta t D)$の部分を$p$による要素と$q$による要素に分解してそれぞれの厳密な積分を行う事で$\Delta t$後の$p_{(\Delta t)}$と$q_{(\Delta t)}$が得られます。
${\rm Exp}(\Delta t D)$を$p$による要素と$q$による要素に分解するとき、これらの要素による作用子が可換であると仮定すると先程の正準方程式を単純に両辺時間積分した形となりますが、実は可換であるとは限らず一般には次式で作用子${\rm Exp}(\Delta t D)$が近似的に展開されます。
\[{\rm Exp}(\Delta t D)\simeq\prod_{i=1}^{N} {\rm Exp}(u_{i} \Delta t D_{p,i}) {\rm Exp}(k_{i} \Delta t D_{q,i})\]
$u_{i}$と$k_{i}$はそれぞれ総和が$1$となる$\Delta t D_{p,i}$と$\Delta t D_{q,i}$の係数です。
このとき係数は一意に定まらないので計算しやすいものを選ぶ余地があります。
ここで$N\to\infty$の極限をとると誤差が無い厳密解となります。
今回作成したシミュレータでは$N=2$として係数$\left(u_{1},u_{2},k_{1},k_{2}\right)=\left(\displaystyle\frac{1}{2},\displaystyle\frac{1}{2},1,0\right)$で計算しています(これは2次のSymplectic積分と呼ばれる)。

(外見だけ)完成した(気持ち悪い)シミュレータ

以下は山が1つのsin波を初期状態としてシミュレーションした結果です

この通りです。
見た目キモすぎてポです。

これがFPUシミュレータで計算された系のエネルギーを散布図化したものです。
青が運動エネルギー、橙がポテンシャルエネルギー、灰色が全エネルギーを表しています。
みた通り気持ち悪いほど全エネルギーの値がぶれてません。


因みに同じ条件を使ってRungeKutta法で計算した系のエネルギーも散布図化してみました。
気持ち悪いぶれ方をしてます。

関連資料とか 

1. E.Fermi, J.Pasta, S.Ulam and M. Tsingou (1955), "Studies of nonlinear probems I", Los Alamos preprint LA-1940 → 本文(ミラノ大学)
2. Joseph Ford (1992), "THE FERMI-PASTA-ULAM PROBLEM:TURNS DISCOVERY", School of Physics, Georgia Institute of Technology, Altanta, Ca 30332, USA → 本文(CBPF)
3. Haruo Yoshida (1990), "Construction of higher order symplectic integrators", Physics Letters A, Volume 150, number 5,6,7 → 本文(CERN)  
4. Symplectic数値積分法
   

1番はFPU問題の最初の論文で1955年に出されたものです。
そして2番はFPU問題のモデルをHamiltonianで書き換え、更なる解析を行ったものです(最期まで読んでませんが大体そんな感じかと)。
今回用いたHamiltonianもこの論文から引用したものを使わせて頂きました。
3番は高次のSymplectic積分に関する論文です。
前述の通りSymplectic積分では非可換な操作を行うため、作用子の係数は簡単に決められるものではありませんが、作用子の係数にフラクタル的構造を見出して漸化式を立て係数を求める事が出来るそうです。
今回は2次のSymplectic積分を適用しましたが、この論文では8次のSymplectic積分の係数まで求められています。  
4番はSymplectic積分を実装する際に参考にしたサイトです(このサイトをChromeで開くと文字化けしました。Firefoxで開くと何ともなかったので文字コードとかの問題かもしれません)

まとめ

Symplectic積分は†光†
あとモードのエネルギーがなんちゃらかんちゃらとかほざいておきながらモード別のエネルギーをプロットしてないのでFFT実装してスペクトル展開的なのをやりたいです(粉蜜柑)。 

自己紹介

"計算機"を作る

携帯ゲーム機を作る

～もしもハードウェアが完成し、ソフトウェアの規格が決まったら～

　・・・まあ、まずはモノを完成させなきゃね。

└(՞ةڼ◔)」＜やっほォ～！したろうだよォ～！

ついに春到来！桜で化粧した東工大構内を散策するのは気持ち良い！
入学式や研究室所属など、様々な「はじまり」や「出会い」があるのもちょうどこの季節ですね。

さて今週、私のもとに一封の封筒が届きました。 
封されていたのは7ヶ月前に参加したサマーチャレンジの報告書でした。
これでもって第8回サマーチャレンジが終了しました。

サマーチャレンジに参加できたのはロボット技術研究会のおかげです。
ついにサマーチャレンジが終わったのでこの機会にこの場を借りて振り返り、感謝をしたいと思います。

サマーチャレンジって？

Fig1.サマーチャレンジで訪れた高エネルギー加速器研究機構。

サマーチャレンジとは夏休み期間中に開催される物理学のイベントです。
主な対象は大学３年生であり、全国の物理学愛好者が茨城県つくば市にある高エネルギー加速器研究機構(KEK) に集結します。
サマーチャレンジはコースが２つあり、「素粒子・原子核コース」と「物質・生命コース」のどちらか一方を選択して応募します。私は素粒子・原子核コースを応募し当選しました。
私が参加した第8回サマーチャレンジは高エネルギー加速器研究機構の主催、
高エネルギー物理学研究者会議、原子核談話会、PF-UA、日本中間子科学会の共催、
日本物理学会、日本加速器学会、日本放射光学会、日本中性子科学会、日本数学会の後援を受けて
2014年8月19日から27日の9日間行われました。 

サマーチャレンジで行われる主な内容は３つあります。
１つ目は「講義」です。原子核や素粒子、宇宙論などを研究されている様々な大学の先生方が自身の研究されているテーマについて講演します。
特に2008年に「小林・益川理論」でノーベル物理学賞を受賞された小林誠先生の講演は興味深かったです。
２つ目は「演習」です。事前にいくつかの実験テーマが与えられ、選択した１つの実験を9日間で取り組みます。
実験では加速器や磁気スペクトログラフなど、普段扱えないような高価な機材や興味深いテーマを扱えるのでとても楽しかったです。
３つ目は「見学」です。 高エネルギー加速器研究機構構内にある加速器やJ-PARCを一日かけて見学するツアーがあります。
現場で活躍している科学者や最先端技術の詰まった装置を間近で触れることができるのでとても興奮しました。

Fig2.研究施設見学ツアーで見学した高エネルギー加速器研究機構の研究施設の一つ、Belle II。

ロ技研で得た“キッカケ”

サマーチャレンジ素粒子・原子核コースの定員は50名であるのに対し、サマーチャレンジの応募は160名ほどありました。
高い競争率でしたが私が当選されたのはロボット技術研究会のおかげだと考えています。

サマーチャレンジの応募にあたって「大学生活で取り組んだ物理に関する機械工作、電子工作、プログラミングなどを高校生にも分かるように800字で説明せよ」という小論文的なものが課せられます。
この大学生活で取り組んだ物理に関する機械工作、電子工作、プログラミングはまさにロボット技術研究会で取り組んできたことではないか！
ロボット技術研究会で幸福の物理として取り組んだ物理に関する取り組みをそのまま800字にまとめ、提出しました。
そのおかげでサマーチャレンジ素粒子・原子核コースに当選し、最先端物理学交流への「はじまり」となったのです。


また、サマーチャレンジに参加することができたおかげで物理学を愛する者たちとの「出会い」 も得ることができました。
今でもサマーチャレンジで交流した人たちと交流があり、先月はサマーチャレンジの演習をともに取り組んだ仲間と一緒に
サマーチャレンジでお世話になった大阪大学の素粒子実験系研究室の人たちのもとへ訪問しました。


こうして、物理学の最先端に触れたり、物理学の壮大な謎に日夜挑んでいる研究者たちと交流できるのもロボット技術研究会で好き勝手やらせてもらえたからです。

最後に

物理学は自然科学の一つです。神秘的であり混沌とした活動の中に驚くほど麗美な秩序が隠されている自然を注意深く扱います。
だから紙の上だけでなく実際に現象を観察したり起こしてみせたりすることも重要になります。
しかし、自分で準備できるのはせいぜい紙とペン。実験装置を作成する環境は大学の研究室に所属しない限りは難しく、諦めていました。

そんな自分にとって、モノをつくる環境と知識を提供し、あとはご自由にどうぞと言ってくれたロボット技術研究会はまさに救いの手でした。

「ロボット」って名前がついているのにロボットとは全く関係のない理学一色に染まった物理学に関する活動を快く受け入れてくれて応援してくれたロボット技術研究会には感謝しています。

これからも宇宙すら支配する物理学を理論と実験の両面から探究したいと思います。

補足

　サマーチャレンジに関することは私のブログの「サマーチャレンジ2014」カテゴリーでまとめられています。
特に開催期間中のことは高エネルギー加速器研究機構の宿舎でリアルタイムに書き連ねたものなので新鮮なデータだと思います。
興味があればぜひ閲覧してみてください。
幸福の物理ブログ：カテゴリー「サマーチャレンジ2014」

　サマーチャレンジの演習をともに取り組んだ仲間たちと大阪大学の素粒子実験系研究室に訪問したときの話は私のブログに書きました。
大阪大学観光名所などを取り上げています。興味があれば閲覧してみてください。
幸福の物理ブログ：大阪大学素粒子実験研究室を訪ねて～したろう二泊三日京都大阪の旅～
 


執筆者：したろう
ブログ：幸福の物理
Twitter：shitaro2012 

こんにちは。幸福の物理の13わつかです。
現在したろうさんと一緒に作成中の「イオントラップ」の進捗状況を報告します。

まずイオントラップとはなんぞや？
イオントラップとは、イオン（電荷を帯びた微粒子） を空気中のある一点に固定させる装置です。つまり物質を宙に浮かすことができる摩訶不思議な装置なのでございます。
原理は簡単で、上下左右に置かれている4つの電極間に交流を流し、電極間に発生する交流電場によって常にイオンの進行方向を変えてやります。するとイオンは運動方向を常に変えながらも、4つの電極のどれにもたどり着かない状況、つまりイオンが空気中で右往左往する状況が完成します。これがイオントラップのおおまかな原理です。
例として物を空中で離すと下に落ちていきますよね。でもその後重力を反転させて反転させて…を繰り返せば結果的に物は空中にとどまるということができます。もちろん、地球の重力を反転させることは無理なので代わりに電場を使ってやろうというわけです。面白いですね。
実際にはイオンに下方向に重力がかかっているので、横方向だけでなく上下方向にも電位差が生じるように電極の形をいろいろ考えてやらないといけないのです。そこで電極の形を決定するために電位のシミュレーションが必要なのですが、シミュレーションにおいて進捗があったので報告します。Mathematicaでの結果が以下の通りです。
 
画像の上の結果で、左右の高くなってるのが電位1vの電極部分、上下の低くなってるのが 接地された電極部分です。正電荷を帯びたイオンはこの鞍型の電位上で、高いところから低いところに向かって移動します。そして高周波でこの鞍の高低が逆転していると思ってください。イオンが中心部分にとどまり続ける理由がわかると思います。
今回の結果は参考にしている論文中の電極形状（双曲線）を再現したまでですが、矩形や円形電極などを一通りシミュレーションしてみて、一番トラップに適すると思われる形状を用いてイオントラップ制作を進めていきます。 

以上13わつかでしたー！