2003年度は11/5(水), 11/12(水), 11/19(水)の3回を大石が担当いたします。よろしくお願いいたします。
開講時限: (水曜5限 16:20〜17:50)
場所: 大宮キャンパス内システム工学部棟3階351教室 (コンピュータ実習室を使う場合もある)
質問: 井戸川先生にお願いいたします。
コンピュータを使った数学の定理証明、理工学のさまざまな分野で使われている数値計算結果の誤差の解析などを目的に精度保証付き数値計算という分野が近年爆発的に発展している。この分野はまさに数学と計算工学の境界分野であり、応用数学の最前線の一分野である。精度保証付き数値計算について、どの分野に進む人にも有益と思われる部分について解説する。数学的な厳密さこそが、最もユニバーサルでポータブルなプログラミング言語となることを理解するのが目標である。数値計算の道具として、MATLABやScilabについて紹介する。これらは一生つきあえるプロの道具である。
数値計算は浮動小数点数を用いて行うのが現代では一般的である。これは，浮動小数点数がハードウエア的にサポートされ，例えば，有理数演算に比べて圧倒的に高速であるからである。
浮動小数点数の表現は，IBM方式など様々なものがあったが，1985年IEEE754が制定され，標準化された。現在のほとんどのCPUはIEEE754に従う，浮動小数点数が実装されている。IntelのPentiumもそうである。
IEEE754とはどんなものか。IEEE754によって標準化されているのは2進の倍精度浮動小数点数である。まず，規格を読んでおく必要がある。浮動小数点数の規格IEEE754はつぎのサイトでその規格書を読むことができる。
IEEE Standard for Binary Floating-Point Arithmetic
この規格書を読んでみましょう。PentiumIII,IVやSunのCPUなどほとんどすべての計算機はこの規格を満たしています。
(A) 佐賀大皆本さんのIEEE754に関する資料(pdf)が詳しくまた丁寧に書かれていて参考になります。
(B) 大石の講義ノート。C言語を利用することを原則とする。環境としてはWindows+Cygwinを考えよう。Cでかかれた数値計算ツール（著者の開発したもの）Slabも用いることにする。Slabの最新版(Slab121_cyg.tzg 約4M)をここにおくので、利用してほしい。次のようなことがWindows上のMATLABでも起きることがわかります。
>> a=0/0
警告: ゼロ割です.
a = NaN
1. IEEE754の企画書の日本語訳 及びIEEE Standard 754 Floating Point Numbersもある。
2. 当事者であるW. Kahanの解説 もある。 制定の歴史を語った文書もある (An Interview with the Old Man of Floating-Point)。
３.また、IEEE754にしたがってどのように浮動小数点演算が実装されているかの例は、例えば、 IntelのPentiumIIIのマニュアルなどを見るとよい。
４.与えれたCPUがIEEE754規格を満たしているかどうかをチェックするプログラムが存在する (TestFloat)。TestFloat1はソフトウエアでシミュレートされた浮動小数点システム SoftFoat と対象としているハードウエアの実行結果を比較して、違いが あるとバグの候補として出力するようなソフトウエアである。
数値計算を行う環境として、各種の数値計算ツールが利用可能です。これを紹介するのも今日の目的です。
ここでは、Scilab(サイラボと読む)とMATLAB(マットラボと読む)について紹介します。また、筆者が作ったSlabもあります。今日はこれらのツールの使い方も勉強しましょう。
Scilabはインタプリタと呼ばれる種類のプログラム言語です。まず、その使い方を勉強してみましょう。
ここにScilabの使い方があります。これを読んでみましょう。ktermを立ち上げて
> scilab &
と打ち込むとScilabの入力画面が立ち上がります。
参考文献
ポアソン分布のグラフ\ref{fig:fig-poi.eps}を描くことを考えよう．
まず，ポアソン分布
P_k(t)=\frac{(\lambda t)^k}{k!}e^{-\lambda t},~(k=0,1,2,\cdots)
を計算するための関数を定義する方法を示そう．
エディタで次を内容とするファイルpoisson.sciを作る．
function [p]=poisson(l,n)
p=zeros(n+1,1);
p(1)=exp(-l);
for i=1:n,
p(i+1)=l*p(i)/i;
end;
このファイルがScilabのパスの通っているフォルダにあるものとしてScilabに呼び出す．
--> getf('poisson.sci')
その後，次のようにすると，図が描かれる．
-->p1=poisson(2,10);
-->p1=poisson(1,10);
-->p2=poisson(2,10);
-->p3=poisson(3,10);
-->n=0:1:10
n =
! 0. 1. 2. 3. 4. 5. 6. 7.
8. 9. 10. !
-->plot2d([n' n' n'],[p1 p2 p3])
レポート、論文、書籍などを書くのにtexが有用である。texの使い方については、つぎのメモがある。
さて、以上のようにScilabにグラフを書かせることは簡単であるが、これを文書の形に残すためのリテラシを簡単にメモしておこう。
Scilabでグラフを描かせた後、ScilabGraphの中のFIleメニューを開く。この中にExportがあるので、これを選択する。Export Typeを
Postscipt-Latexに選びOKとすると、保存するデイレクトリとファイル名を聞かれるので、例えばexp2.epsというファイル名
にして保存する。platex2etexのファイルでは、
\documentclass{jbook}
\usepackage{graphicx}
とスタイルファイルとしてgraphicx.styを利用する。図を描かせる位置で以下のように引用する。
\begin{figure}[htbp]
\begin{center}
\includegraphics[keepaspectratio=true,height
=60mm]{exp2.eps}
\end{center}
\caption{保留時間}
\label{fig:exp2.eps}
\end{figure}
こうして、図がtexファイル中に引用できる。図を描かせるには様々な方法があるので、これはあくまで一つのスタイルである。
Scilabはフリーのツールである。MATLABと呼ばれる商用のソフトウエアもある。これはグラフィックスもより高度で、ディジタル信号処理ツールボックスなどもある。また、数値計算ツールのデファクトスタンダードとなる文法を提供してきたことも見逃せない。
参考文献
Calcの使い方について述べて行く。最初は電卓のように感じるが、後にその奥の深さがわかると思う。
calcをインストールすると、コマンドラインから、
$ calc
と入力するとCalcが起動する。helpと入力するとhelpにどのようなトピックスがあるか表示される。そこで、例えば、help fullとすると、すべてのhelp情報を見ることができる。
こうして、立ち上がった状態で
> 5*(2+9)
> 55
のように、数式をC言語のように記述すると、電卓のように答えを返す。ここでさらに、
> .*10.5
> 577.5
とすると、最後に計算された値が'.'というシンボルによって呼び出され(今の場合は55)計算を続行することができる。Calcでは表示上は数は小数点で表され
るが内部ではすべて有理数で記憶されている。変数を導入することもできる。
> a=.
とすると。現在の値が変数aに記憶される。
ここで、階乗の計算を行ってみよう。
> a=100
とする。a!を計算するには次のようにする。
> fact(a)
93326215443944152681699238856266700490715968264381621468592963895
21759999322991560894146397615651828625369792082722375825118521091
6864000000000000000000000000
この計算は正確であり、浮動小数点演算ではできない種類のものである(桁数が多すぎる!)。
小数点以下何桁まで精度を保つかも指定できる。
> config("display", 50)
> epsilon(1e-50)
としておけば、次のようにしてcos(1)の値が小数点以下50桁まで正しく計算される。
> cos(1)
~.54030230586813971740093660744297660373231042061792
初等関数を計算する精度が制御できることがわかる。
参考文献
連立一次方程式
1.1x + 2.2y - 3.3z = 1
2.1x - 3.2y + 3.3z = 2
-3.1x + 4.2y + 4.3z = 3
をScilabとMATLABとCalcで解け。
>> b=a^0
b = 1
64台のクラスタを利用して10000次元までの蜜係数行列をもつ連立一次方程式を解いた結果を述べた。
PCクラスタ上での連立一次方程式の精度保証付き数値計算(大石進一、荻田武史)
その結果、Oishi-Rumpの提案した高速解法は5000次元の連立一次方程式(良条件)ぐらいまでは有効であることがわかる。また、10000次元でも高速化していないアルゴリズムなら解けることもわかる。
などについて論じた。
用語 MPI (Message Passing Interface) 、BLAS(Basic Linear Algebraic Subroutines)
Bauer-Fikeの定理を利用して固有値の精度保証を行う方法について解説する。次の解説（電子情報通信学会2001年1月号掲載)
大石進一:実用段階に到達した精度保証付き数値計算(dviファイル、約100k)
も参考にされたい。また、講義では、成分ごとのBauer-Fikeの定理を示し、それを利用する方法についても解説する。この方法では局所的に固有値が一意であることの十分条件も示される。
まず、次の小島政和氏(東京工業大学教授)の資料（数理計画法の最近の発展 --- 内点法を中心にして ---）(pdf,100k)を見てみよう(またはより詳しい解説)。小島氏は内点法の研究で大変優れた成果をあげておられるが、その主双対内点法はまさに精度保証の理論を利用しているといっても過言ではない。ここでは双対定理を利用した計画問題の精度保証法について議論する。
微分方程式の解の存在証明(常微分、偏微分)
ケプラー問題など証明の最後のつめに使う。ケプラー予想の証明の紹介を情報科学演習(情報学科3年生必修科目)で，河邊昌彦君，内藤一兵衛君，八幡淳君がしてくれた。彼らは実際に区間演算プログラムを作り，計算機援用証明の部分も追試(pdf.0.5M)を行っている。以下はその発表の時の河邊昌彦君の草稿である。河邊君の許可を得てここにpdfファイルをおく。pdfで約1Mのファイル(河邊昌彦氏pdfファイル)である。(2001/6/9)
カオス(Conley Index(Mischaikowのサーベイpdf)など構造安定な量の計算) 談話会終了後、石井(九大)さんからLorenzアトラクタの存在検証がW.Tuckerによってなされていることを教えていただいた。これはSmaleの第14問題(Lorenz微分方程式のアトラクタはWilliams, Guckenheimer, とYorkeが提案したgeometric Lorenz attractorであるか?)と関連する。そして、 Tucker (2002) はこれを肯定的に解いた(pdf,180k)。（中尾先生のご指摘、関連のページ)
最近の精度保証関係の国際シンポジウム(2002 SIAM Workshop on Validated Computing)のWorkshopプログラムを見ることができる。
まず、Scilabによる非線形方程式の解法について述べる。
Newton法の解法ルーチンの例(NEWTON2.sci)
を示す。これはScilabのプログラムである。これには自動微分が使われている。自動微分のプログラムは次のようになる。
さらに、要素としてなんでも型が入る行列クラスが必要となる。それは、つぎで定義される。
これらをダウンロードしてフォルダSciproに入れておこう。そして、非線形方程式としてつぎの例を考える。
非線形方程式をNewton法で解くにはつぎのようにする。
-->getf('E:\Documents and Settings\Administrator\デスクトップ\九大\SCIPRO\NEWTON2.SCI');
-->getf('E:\Documents and Settings\Administrator\デスクトップ\九大\SCIPRO\MAT.SCI');
-->getf('E:\Documents and Settings\Administrator\デスクトップ\九大\SCIPRO\AUTODIF.SCI');
-->getf('E:\Documents and Settings\Administrator\デスクトップ\九大\SCIPRO\FUNC.SCI');
-->x=newton(2,'func')
x =
x(1)
!mat dim !
x(2)
! 2. 1. !
x(3)
0.9995920
x(4)
0.0285618
このようにScilabにおけるクラス（オブジェクト）はt-listを用いて定義される。
もう少し複雑な例を解いてみよう。つぎの非線形関数のゼロ点を求めてみる。
function y=func(x)
y=x;
y(3)=x(3)*x(3)-x(4)*x(4)-3*x(3)+2;
y(4)=2*x(3)*x(4)-3*x(4)-1;
y(5)=3*sin(x(3)+x(4)+2*x(5))+x(5)*x(5)-3;
y(6)=x(6)*x(6)+3*sin(x(3))-3;
y(7)=cos(x(7))*x(7)+x(3)*x(3)-x(5)*x(5)+x(4)-2;
y(8)=x(8)+2*x(3)-x(7)+x(6)+x(3)*x(4)-5;
endfunction
この関数定義をつぎのファイルに入れた。これくらいになると、手で微係数を求めるのは少々めんどくさくなる。自動微分のありがたみがわかる。
さて、Scilabで解いてみよう。
-->getf('E:\Documents and Settings\Administrator\デスクトップ\九大\SCIPRO\FUNC3.SCI');
Warning :redefining function: func
-->x=newton(6,'func')
x =
x(1)
!mat dim !
x(2)
! 6. 1. !
x(3)
2.3002426
x(4)
0.6248105
x(5)
- 0.4799244
x(6)
0.8737101
x(7)
14.396064
x(8)
12.484653
いよいよKrawczyk作用素を定義してみよう。それはつぎのようになる。
著者が考案した幾つかの方法を議論する。
つぎの論文に示した方法である。
Shin’ichi Oishi: Numerical verification of existence and inclusion of solutions for nonlinear operator equations, Journal of Computational and Applied Mathematics, Vol.60, 1995, pp.171-185
非線形楕円型方程式等を含む非線形作用素方程式の解の存在・唯一性・存在範囲を精度保証付き数値計算にもとづく計算機援用証明によって証明する方法を示した。また、そのためのソフトウェアライブラリを作製し、検証例を示した。
この論文を書いた頃(1993)は著者は有理数演算を用いていたが、現在では倍精度浮動小数点数を用いている。この方法でDufiing方程式を解くには次のようにする。以下、Scilabのプログラムである。
Duffing方程式の決定方程式と精度保証プログラム(DUFFING.sci)
これにはフーリエクラス(FOURIER.sci)が用いられている。
常微分方程式のホモクリニックオービットやヘテロクリニックオービットの数値的存在検証のために開発した方法。
Shin'ichi
OISHI: Two Topics in Nonlinear System Analysis through Fixed Point Theorems,
IEICI Trans. Fundamentals Vol.E77-A,
No.7, (1994) pp.1144-115
に一般の境界値問題の場合に対する検証法を書いた。また、次の論文に常微分方程式のホモクリニックオービットやヘテロクリニックオービットの数値的存在検証法を示した。
Shin'ichi Oishi: "A Numerical Method of Proving the Existence and Inclusion of Connecting Orbits for Continuous Dynamical Systems", Journal of Universal Computer Science, vol.4 no.2, (1998), 193-201
参考文献 非線形解析入門 大石 進一 著 コロナ社 2,800円 1997/04発行 ISBN4-339-02600-X
これについてはつぎの文献を参照しながら議論する。
参考文献 精度保証付き数値計算 大石 進一 著 コロナ社 2,200円 1999/12発行 ISBN4-339-02605-0
最近は村重氏(東京大学)と特異積分方程式の解の数値的検証法について研究を進めている。
このページのURIはhttp://www.oishi.info.waseda.ac.jp/~oishi/siba/siba.html
最終更新 2003/11/5