kashiの日記

2018/11/29(木)kv-0.4.47

kvライブラリを、0.4.47にアップデートしました。

実は、日記には書きませんでしたが、昨日、kv-0.4.46にアップデートしています。

kv-0.4.46では、精度保証付きODE solverのstep sizeの制御部分が自分の意図通りに書かれていなかったいくつかのバグを修正しました。以前のversionでも精度保証が破れていることはないはずですが、step sizeの戦略が変わるので計算結果は0.4.45以前とは少し変わります。特にstiffな問題に関しては、大きく変わる可能性もあります。

kv-0.4.47もODE solverに関する変更です。ユーザから見れば、ode-maffine0.hppが削除されたこと以外はほとんど変わっていません。しかし、内部には大きな変更が入りました。kv-0.4.35で行った変分方程式の解き方の改善は、ode-maffineに対するものでしたが、それをode-autodifに直接組み込むようにしました。これにより、(ほとんど使われていないでしょうけど)ode-qrがstiff ODEに強くなっていることが予想されます。

変更が大きいので気分的にスナップショットを残しておきたくて、短期間でのversion upになりました。いろいろチェックはしたつもりですが、何か不具合がありましたらお知らせ下さい。

コメント（0件）

2018/11/24(土)LohnerのAWAで遊ぼう

精度保証

AWAという常微分方程式の初期値問題を精度保証付きで計算する、一部ではとても有名なソフトウェアがあります。Rudolf J. Lohner氏がアルゴリズムを考案し、ソフトウェアとして実装したものです。極めて先駆的な研究であり、その方法は通称Lohner法と呼ばれ、その後の研究の基礎になったものです。また、実装が大変優秀で完全に実用レベルに達しています。ソフトウェアそのものは2005年で更新が止まっていますが、現在でも一線級の実力を持ちます。なお、AWAは、「anfangswertaufgabe」の3文字で、ドイツ語で「初期値問題」の意味だそうです。

古いソフトウェアなので手元で動かした経験のある人は少ないと思いますが、今回は、この知る人ぞ知るソフトウェアを手元で動かして、実際に初期値問題の精度保証を行うまでの方法を書きたいと思います。

AWAのホームページは AWAにありますが、ここには古い機種でのバイナリしかありません。awa.src.tar.Zはありますが、ここには拡張子.pのファイルしか入っていません。AWAは、精度保証付き数値計算のために開発されたPascal-XSCという言語で書かれています。よってAWAを使うには、まずPascal-XSCコンパイラを準備しなければなりません。幸いなことにPascal-XSCのパッケージはその中に丸ごとAWAを含んでいて、Pascal-XSCをmakeすればAWAが手に入ります。

Pascal-XSCは、XSC LanguagesのページにあるPascal-XSCの欄の、Free Download (Compiler and Toolbox)にあります。そこに、ソースコード ( p-xsc-3.6.2.tar.gz (released 2005-12-19) ) があるのでそれをダウンロードしましょう。

以下、makeの仕方を示します。Ubuntu 18.04 (64bit)で行いました。32bitのbinaryを作成するので、

sudo apt install build-essential
sudo apt install gcc-multilib

が必要だと思います(後者が重要)。archiveを展開するとMakefileがあり、まずMakefileを書き換えます。どういうわけだか単にmakeするだけでいわゆるmake installに相当する作業を行い更にpathまで切ろうとするというお行儀の悪さなので、少し手を入れます。まず、Makefileの最初の方の

PREFIX=/usr/local/p-xsc

の部分を書き換えて、書き込んでもよい場所を絶対パスで指定します。自分は念の為ホームディレクトリの中の適当な場所を指定しました。また、Makefile中の

        @echo "****************************************************"
        @echo " Install System PATH in /etc/profile.local as root *"
        @echo " or in ~/.bashrc as user                           *"
        @echo "****************************************************"
        cd profile; ./uidtest.sh ${PREFIX}

の部分で/etc/profile.localというシステムファイルに書き込もうとするので、ここは#でコメントアウトしておきます。

makeそのものは自分の環境では問題なく出来ました。PREFIXで指定した場所の bin/awa が目的のAWAのコンパイル済み実行ファイルです。Pascal-XSC言語を使うつもりがないならば、

生成されたawaの実行ファイル
prts/rts/o_msg1.h (Pascal-XSCの"runtime message"。current dirに置く)
awa/examples/* (AWAのサンプルファイル)

だけどこかに保存して、残りは消してしまってもいいでしょう。

さて、実際に初期値問題を精度保証する方法を示します。例題は、van der Pol方程式を1階に直した

\begin{align*} \frac{dx}{dt} &= y \\ \frac{dy}{dt} &= \mu (1-x^2)y - x \end{align*}

を、\mu=1とし、初期値 (x(0),y(0))=(1,1)でt\in[0,20]の範囲で計算することにしました。まず、test.inを

2
test.dat
test.out

の内容で作成します。「2」は方程式の次元、test.datは方程式を記述するファイル、test.outは出力するファイル名です。次に、test.datを、

F1 = U2
F2 = 1 * (1-U1*U1) * U2 - U1
;;

1 0 20 0 24

4 0

1 1

1e-16 1e-16

の内容で作成し、

awa < test.in

とすると、test.outが生成されます。この場合は、

ode read from file test.dat 


                 out = 1
             T_start =  0.000000000000000E+000
              T_end  =  2.000000000000000E+001
initial stepsize  h  =  0.000000000000000E+000
        order     p  = 24

enclosure method :  intersection of interval-vector and QR-decomposition

output :  function values


t =  0.000000000000000E+000
[ 1.00000000000000E+000, 1.00000000000000E+000]  0.00E+000  0.0E+000
[ 1.00000000000000E+000, 1.00000000000000E+000]  0.00E+000  0.0E+000
          h =  1.911222934722900E-001

t =  1.911222934722900E-001
[ 1.16970857077089E+000, 1.16970857077090E+000]  6.67E-016  5.7E-016
[ 7.61501247518094E-001, 7.61501247518095E-001]  3.34E-016  4.4E-016
          h =  1.831474304199219E-001

(中略)

t =  1.993603903055191E+001
[ 2.00246009621560E+000, 2.00246009621563E+000]  2.18E-014  1.1E-014
[ 1.69827498867334E-001, 1.69827498867578E-001]  2.43E-013  1.5E-012
          h =  6.396096944808960E-002

t =  2.000000000000000E+001
[ 2.00848791779841E+000, 2.00848791779843E+000]  8.00E-015  4.0E-015
[ 2.32898543064796E-002, 2.32898543066811E-002]  2.02E-013  8.7E-012
Number of integration steps: 134

Time:       1.64 sec

のようなファイルが生成されました(計算速度はLet's Note RZ4で)。時刻毎のx,yの値が区間で表示されていますが、真の軌道がこの区間内にあることが数学的に厳密に保証されています。区間の後の2つの数値はそれぞれ区間幅、相対区間幅を表しているようです。x,yそれぞれの上限、下限をプロットしてみましたが、このくらいの簡単な問題だと差はグラフ上では見えません。