AWAはPascal-XSCというマイナーな言語で書かれていますから、メジャーなMATLABで書かれていた方がありがたい、という人もいるでしょう。解を初期値に関して高次に展開することによる非常に高性能な精度保証付きアルゴリズムとして知られています。今回は、この両者を実際に動かしてみたレポートです。MATLABはLinux上でR2018aを使いました。PayPalで支払うと、メールで6時間有効なダウンロードリンクが送られてくるので、INTLAB.zipをダウンロードしておきましょう。中身のファイルそのものは特にプロテクトがかかっているわけではなく自由に使えます。Intlabのインストールは簡単で、zipを解凍して好きなディレクトリに置くだけです。初回はいろいろ調査するため時間がかかります。toolboxはoctaveに対応していません。octaveでのテストを行わずに開発していたそうで、またいくつかのoctaveの非互換の部分が問題ですぐに対応するのは難しいそうです。次に実際に動かしてみます。(x(0),y(0))=(1,1)でt\in[0,20]の範囲で計算する問題にしました。toolboxの方の使い方を説明します。toolboxのページを見て下さい。最小限のコードはこんな感じだと思います。ファイルを作成し、MATLAB上でhogeで実行できます。mu*(1-x(1)^2)];endこのように、ODEの右辺を定義する関数(ここではvdp)と、そのヤコビ行列を返す関数(ここではvdp_J)を定義し、関数awaを呼びます。awaの計算したデータを分かりやすく表示してくれるのがawa_dispです。ヤコビ行列を定義するのが面倒ですし、Intlabの持っている自動微分機能でなんとかなりそうな気がするのですが、高速化のためにやむを得なかったそうです。typeadjustという関数が使われていますが、関数の戻り値に定数があるときはこうやって型を補正する必要があります。関数awaの第5引数にawaset('order',10,'h0',0,'EvalMeth',4,'AbsTol',1e-16,'RelTol',1e-16)のような関数を使って、様々なパラメータを変更することができます(これらの値はデフォルト値)。変更したいパラメータの名前を奇数番目の引数に、値を偶数番目の引数に書きます。mu*(1-x(1)^2)];endを実行してみました。すると、次のような解が得られました。8.7E-012でした。どちらも精度保証付きソルバーなので、両者の共通部分に真の解があることになります。計算時間はオリジナルが0.81secだったのに対してこちらは20.97secと、かなり遅くなってしまっています。Bünger先生によるとかなり高速化の工夫をしたとのことですが、これがMATLABの限界でしょうか。toolboxこちらの最小限のコードはこんな感じだと思います。toolboxを見て下さい。endendこちらはawaの方にあったヤコビ行列を書く必要はありません。xに続く第3の変数が増えており、そこに自然数が与えられたらその番号に対応する右辺の式だけを計算して返すことが要求されています。これも高速化のためです。Xrと3つに増えており、またverifyode_dispには初期値も渡す必要があります。'h0',0,'loc_err_tol',1e-10,'shrinkwrap',1,'precondition',0,'blunting',0)のようなオプションを第4引数に指定することによって動作をカスタマイズできます。AWAと違って非常にオプションが多く、またオプションの指定によって大きく性能が違うようなので、まずどんなオプションがいいか調べてみました。Polでt=15まで計算、orderは18に固定して、区間幅の変化を調べてみました。それを行った結果が次のグラフです。Toolboxの計算結果で、下3つは他のもの(参考値)です。カッコ内の数値は計算時間(秒)です。グラフの説明は、それぞれ次のようなオプションと対応しています。横軸がt、縦軸は計算結果の区間幅を表しています。(QRのみ)が一番いいことが分かります。これがデフォルト値なのか理解に苦しみますが。設定して、以下を実行してみました。endendすると、次のような解が得られました。2.71e-12少し精度が落ちているものの、同様に精度保証された解が得られました。計算時間は25秒弱とやや遅めです。0)は試行錯誤の結果ですが、もっと良い設定があるかもしれません。方程式や初期値が違えばまた最適値も変わる可能性が高いです。modelは前評判からするとやや性能が出ておらず、少しがっかりしています。kvライブラリと比較してみました。\end{align*}です。kv::ode_param<itv::base_type>().set_verbose(1));}です。toolboxは、関数のヤコビ行列を書くのが面倒で試しませんでした。計算結果を図示します。このようになりました。精度保証付きソルバー、当たり前ではありますがグラフがぴったり重なっています。toolboxはt=10までは計算できましたが、t=11では計算が止まりませんでした。AWAはt=12過ぎで解が発散しました。kvはt=17過ぎまで計算できました。解の区間幅を図示すると、このようになりました。kvが長い時間区間幅が発散せずに持ちこたえていることが分かります。toolboxが5時間20分17秒、AWAが17.9秒、kvが22.9秒でした。kvは、他の2つには不可能な、内部の演算型を全てdouble-double(疑似4倍精度)に変えて精度保証付き計算を行うこともでき、この場合は軽々とt=30以上まで計算することができます。kvの圧勝なわけですが、他のソルバーかかってこいや。他の例題、特にstiffな例題や、他のソルバーを試してみたいと思います。応用数理セミナーInformation精度保証付き数値計算のことなど、いろいろつぶやく予定。