S-TCI プログラム - なぜか数学者にはワイン好きが多い

日本VR学会で使った，S-TCIを使ったグラフのプログラムを公開します．Mathematica 6.0で使いました．
対象は，これ以上簡単な例は無いというくらい簡単な，単なる調和振動子です．外力も無し，ダンパ定数なんかも無し，バネも質点もイッコです．質点の質量も1とします．つまり，ハミルトニアンで言うと

$H=\frac{1}{2}(p^2+k q^2)$

です．

STCI[tau_] :=
Module[{r = 0.6, v = 0, k = 1, list = {}},
Do[rr = r; vv = v;
r = (1 + 1/4 k tau^2)^-1 ((1 - 1/4 k tau^2) rr + vv*tau + 1/2 k tau^2);
v = (1 + 1/4 k tau^2)^-1 ((1 - 1/4 k tau^2) vv - k*(rr - 1) tau);
list = Append[list, r], {i, 0, 10/tau}]; list]

この関数を使って，例えば次のようにグラフをプロットできます．

ListPlot[STCI[0.1], Joined -> True, PlotRange -> {0.4, 1.4}]

これだけじゃ面白くないので，普通の1次の陽的オイラー法も関数化しました．

EULER[tau_] :=
Module[{r = 0.6, v = 0, k = 1, list = {}},
Do[rr = r; vv = v; r = rr + vv*tau; v = vv - k*(rr - 1)*tau;
list = Append[list, r], {i, 0, 10/tau}]; list]

そして，Mathematica 6.0の新機能を使うと，次のようにアニメーション化できます．

Manipulate[{ListPlot[STCI[t], Joined -> True,
PlotRange -> {0.4, 1.4}],
ListPlot[EULER[t], Joined -> True, PlotRange -> {0.4, 1.4}]}, {{t, 0.001, "? t"}, 0.001, 1, 0.01, ImageSize -> Large}]

これを見ると，オイラー法だと刻み幅?tが0.001くらいじゃないとバネの振幅の1周期も持たないのが，S-TCIだと刻み幅が1.0くらいでも大丈夫なことが分かります．これくらいは普通のシンプレクテイック数値解法でもイケるとは思うのですが，ノートPCでのプレゼンではグラフを3つ重ねるのはきついと思って二つにしました．

今年の応用数理学会でも，東大に行ったんでしたっけ，松尾さんが「シンプレクティック性の保存とエネルギー保存は同時に出来ないという宿命が．．．」と発言なさっておりましたが，確かに現在のS-TCIでもシンプレクティック性は保存できません．
で，その時，私は「シンプレクティック性の保存は，そもそも厳密にエネルギーを保存しなくても大丈夫なところに特徴があるので，両立しなくても良いのでは」と思ったのですが，発言する時間がありませんでした．
しかし，私は近似的に両立が可能だと思っています．アイディアは秘密ですが．．．