ディレクトリ以下は次のようになっています。 equation/ Pfaffian 方程式生成プログラム hgd/ ホロノミック勾配法プログラム それぞれについて説明します。 ○Pfaffian 方程式生成プログラム ディレクトリ equation にはプログラム mydso9.rr があります。 入力: tmp-dso3.ab 出力: save9/pf.ab pf.ab は asir の bsave でセーブされたオブジェクトであり、 そのフォーマットは次のようになっています。 Pfs = bload("pf.ab")$ Pfs[0]: 係数行列 9 個のリスト Pfs[1]: Pfaffian のベース (B = [1,dx_3_1,dx_3_2,dx_3_3]) Pfs[2]=HashA, Pfs[3]=HashB: ハッシュ表 Pfs[0][I] は変数 GV[I] (=x_*_*) で縦ベクトル B.u (u は未知関数)を 偏微分したときの係数行列です。 各行列の成分は、すべて変数 b_0 〜 b_136 で書かれています。 変数 b_0 〜 b_136 は、それぞれ a_0 〜 a_225 の多項式、 a_0 〜 a_225 は、x_1_1 〜 x_3_3 の多項式です。 ハッシュ表の構造は次のようになっています。 HashA = [ [ a_0, x_* の既約多項式], [a_1, x_* の既約多項式], ...]: 連想リスト HashB = [ [ b_0, a_* の既約多項式], [b_1, a_* の既約多項式], ...]: 連想リスト 以下、覚え書き 計算の手順は次の通りです。 G = tmp-dos3.ab は多項式が多すぎるので、次のように間引きます。 (1) リーディングモノミアルが一致する多項式は、一つだけ残して他は捨てる。 (2) リーディングモノミアルが dx_1_1 * dx_3_1 より大きな多項式は、 パッフィアンの左辺には現れないので、全て捨てる。 これで、545 個の多項式を 30 個まで減らせます。 残った多項式の係数を降順に全て並べて行列を作ります。 行列なので連立方程式と思っても構いません。 (dx_1_1 * dx_3_1 より小さな全てのモノミアルを変数と思う) 既に掃出しが終わっているので方程式が解けることは容易に確認でき、 実際にこの連立方程式を解いてやればパッフィアンの係数が全て得られます。 ここは線形代数なので適当に変数の置き換えをしても安全です。 ○ホロノミック勾配法 ディレクトリ hgd には以下のファイルがあります。 hgd.rr ホロノミック勾配法プログラム oh_runge_kutta.rr 部分プログラム(４次 Runge-Kutta) oh_hash.rr 部分プログラム(Pfaffian処理) hgd_test.rr テストプログラム hgd_test16.rr テストプログラム(小数点表示16桁) check.rr 清さんの結果との比較検証プログラム etrcalc.r 清さんによる数値積分プログラム Log.txt hgd_test.rr の実行ログ Log16.txt hgd_test16.rr の実行ログ hgd.rr がホロノミック勾配法の本体です。S^d 論文の Algorithm1, すなわち、 grad 方向に探索するアルゴリズムを実装しています。データは hgd.rr の関数 hgd_test 内に埋め込んであります。積分値、微分値の生成には etrcalc.r を 利用します。 hgd_test.rr は comets データについて、差分幅 1.0*10^(-4) で極小値を探索したものです。 asir -f hgd_test.rr と打ち込んで実行します。出力ログが、Log.txt です。 もともとの清さんの結果では、極小点が見掛けの特異点(x31=0)上にあり、 このプログラムでは計算不可能なので、perturb した点から出発してあります。 以下、check.rr の説明 Xs=[0.333,-0.033,0.448, 0.367,-0.279,-0.752, 0.01,-0.081,-0.247] Ys=0.8367516489 (etrcalc.r で求めたもの) を初期点として、ホロノミック勾配法を実行すると、 Xp=[0.32825,-0.0311452,0.44599, 0.36469,-0.278664,-0.752192, 0.00146575,-0.0791232,-0.247411] Yp=0.8367354795509646 で極小点となります。Xp での関数値は etrcalc.r では Yp=0.8367411574 となり、誤差は、0.56*10^(-5) です。