数値計算(情報学科3年必修)

数値計算の講義資料並びにレポート提出状況(大石担当分)など講義情報を掲載するページです。講義の予習,復習,課題作成に有効に活用してください。また,このページの内容が正式シラバスとなります。

教室変更

教室に入りきれなかったので,次のように教室を変更します。講義の時間は変更しません。

水曜3限 54-201(123) -> 52-101(180)

金曜4限 52-104(150) -> 52-202(180)

教科書

次の2冊の教科書を使います。必ず用意してください。既に,生協には入荷しています(2001/4/9)。

[1] 大石進一:数値計算(裳華房) 1999

[2] 大石進一:Linux数値計算ツール(コロナ社) 2000

講義予定

  1. 4/11,13 (大石) 1 浮動小数点数(IEEE754)
  2. 4/18,20 (柏木) C言語(I,II)
  3. 4/25,27 (柏木) C言語(III,IV)
  4. 5/9,11 (柏木) C言語(V),C++言語と2.区間演算
  5. 5/16,18 (大石) 3.連立一次方程式とその数値計算ツール(1,2)
  6. 5/23,25 (柏木) 4,5 連立一次方程式の数値解法とその精度保証
  7. 5/30(柏木), 5 連立一次方程式 6/1 (大石) 6,7 多項式・補間
  8. 6/6,8 (大石) 7. 補間 8 積分 とその数値計算ツール(4,7)
  9. 6/13,15 (柏木) 9 ニュートン法
  10. 6/20(柏木)自動微分,22 (大石)数値積分
  11. 6/27,29 (大石) 数値計算ツール(,5,6,8)
  12. 7/4(柏木),7/6 (大石) 12 常微分方程式
  13. 7/11(大石),13 (柏木) 数値計算ツール(3,9,10)

左側の数字nは第n週の講義を表します。担当教員の右側に現れる数字mは教科書[1]の第m章を講義することを表します。また,数値計算ツール(i,j)という表示は教科書[2]の第i章とj章を講義することを表します。また,C言語について基礎から講義することにしました(2001/4/13)。

使用言語

C及びC++を用いる。FORTRANに言及することがある。Javaについても言及することがある。数値計算ツールMatlab,Scilab,Octaveなどのインタープリタ言語も用いる。

CはRedhat6.2JまたはVine2.1のgccを標準とする。C++はRedhat6.2JまたはVine2.1のg++を標準とする。これらは第4,5端末室での標準である。第4,5端末室はリニューアルされ,Redhat6.2JベースのLinux環境であるので,課題の作成はこの両端末室の利用を前提としている。第4端末室は現在の所,水曜3,4限,金曜,3,4限が実験のため使用不可となっている以外は利用可である。また,第5端末室は現在のところいつでも利用できるようである(ただし,講義での利用申請があると制限が増えることがあるので各自確認のこと)。Matlabは第7端末室で利用可(リモートも可)。

課題

各週に必ず課題が出題される。課題はすべて提出しないと単位は出ない。また,期末には試験を行う。その合計で成績の評価を行う。 尚,大石分の課題はこのページに掲載する。柏木先生分は柏木先生より別途指示する。

第1回課題

第一回課題(4/13出題) 第一回課題の解説(dvi-file)

第2回課題

第2回課題(優秀作を決定(6/2))

第3回課題

今回は力作揃いだった。例として,優秀なレポートを幾つか掲載する。 すべてdviファイルで100k程度。

福田 浩章 君

河邊 昌彦 君

益子 理絵 君

斎藤 卓也 君

宮之前 謙一 君

森野 岳宏 君

第4回課題

今回の課題は,出題の意図がつかみにくかったようだが,実際には意図は講義の際に説明している。

問1は最適化BLASを作る場面を想定し,その一つの関数であるdgemmを作ることにある。したがって,目標は標準BLASのdgemmより速いものをつくることである。既存の最適化BLASにかなわなくても,それと比較しうる速度をもたせることができれば目標達成である。

問2の意図はBLAS関数(dgemmなど)を高速化したら,連立一次方程式の解法関数が高速化される仕組みを理解することにある。 これは分割解法を利用する点にあるのだが,それがどれだけ理解できたかを問うている。

以上の意図が完全に理解された見事なレポートの例を示す。河邊昌彦氏のレポート(dvi)である。また,斎藤卓也氏のレポート(dvi)もいいセンスである。

永田貴彦氏のVisual C++にインストールしながらのレポート(dvi)は推理小説を読むようで,引き込まされた。インストールはどうもゲームをクリアすることに似ているところがある。LAPACK ver3は最近完成したLAPACKの最新バージョンでver2よりもいろいろ改良されているプロのお薦めバージョンである。

任意(提出しなくてよい)の問題

「Rlabをglibc2.1系にインストールせよ。もしできたら,一番最初にできた人は試験を免除にする。」という問題を出したところ,早速,匂坂君から,4,5端末室のRedhat6.2では何の問題もなくできるとのメールをもらった。確認したところやっぱりできる。 Vine 2.1ではどうであろうか(6/23)。

第5回課題

今回の課題は「数値計算用インタプリタを作れ 」である。このインタプリタを記述する言語は原則としてCとする。GUIなどにJavaを利用して,ポータビリティを持たせることは問題ない。一部C++なども利用してよい。

以下,Linux,Windows+cygwin,Mac OSX(Developing Toolはインストール済みとする)などUnix系OSを利用していると想定する。これらの環境では,通常bison(yacc),flex(lex)などのコンパイラ作成ツールがもともとインストール済みである。また,メモリのガーべージコレクションのためのgcもGPLライセンスで利用できる(のでダウンロードして利用できる)。ここでは,これらのツールを使うことを前提に数値計算用インタープリタを作っていく方法を説明する。尚,今回は,最低40点から最高300点までレポート点がつく。例えば,200点がつけば,(大石の出題した)他のレポートを1つ出していなくても,それは100点になり,提出しなくてもよい。更に,点が余れば,試験点にも補充される。

C言語を用いて数値計算用のインタプリタを作るとき,flexは字句解析用のCプログラムを自動生成するツールで,bisonは構文解析用のCプログラムを自動生成するツールである。それぞれ,lexおよびyaccの上位コンパチである(使い方はbisonとyaccで若干異なる)。flexとbisonについてはweb上に解説があるのでそれらをまず読んでみて欲しい。

(解説1)flexの説明

(解説2)bisonの説明

これらには簡単な電卓用のプログラムが附いている。これを実行させて解説すれば40点となる。今回のレポートの最低ラインはここである。すなわち,

(1) bisonの機能を用いて,かけ算わり算を足し算引き算より優先させて計算する倍精度浮動小数点数を用いた電卓を作れ(40点こここまでだったら数時間もかからずにできるであろう)。

(1+) 実部,虚部がともにdoubleの複素数が扱える,電卓をつくれ(以下との兼ね合いがあるが+30点程度,楽な演習)。

(1++) 行列(double)が扱えるようにする。ただし,加減乗算にはBLASの関数を用いること。(+30点程度,これがいままでの講義との関連の最初)。

(1+++) 行列(複素数)が扱えるようにする。ただし,加減乗算にはBLASの関数を用いること。(+30点程度)。

(1++++) 変数を扱えるようにする(ここまででメモリ電卓となる+10点)。

(1+++++) sin,cosなどやべき乗もできるようにする(ここまででメモリ付き関数電卓となる+10点)。

(2) 制御文(if文やwhile文やfor文)を扱えるようにする。(+40点)

(3) ユーザが関数を定義できるようにする。(+40点)

(2),(3)の問題を解くためには,解説1,2では不十分となる。これらができるようになるための説明のあるwebをリンクしよう。

ずばり,面白いページは次である。

前橋和弥氏のページ ここの「電卓を作ってみよう」がとても面白い。

これを丁寧に解説して,コピーすると上の評価基準ではほぼ満点かそれを越える。もちろんコピーしただけでは満点は保証できない。ということで

(M+) 前橋氏の電卓にsin,cos,べき乗などの関数を組み込んでみよう。また,ユーザ定義関数の返す値を指定できるようにしよう。(これぐらいできると100点)。

というような見当だろう。これは,一日ぐらい考えればできる。

(M++) 前橋氏の電卓に行列ができるようにしてみよう。もちろん,加減乗算はBLASの関数を使う。また,ついでだから,

Ax=b

> x=A\b

でできるようにしてみよう。ここまでできたら200点。複素数も扱えれば+20で,複素行列も扱えれば+30とする。複素行列も扱えるようになったら,固有値問題も解けるようにしてみたい。

> [P D] = eig(A)

で行列Aの固有値(を対角要素に並べた行列)Dと対応する固有ベクトルを並べた行列Pが計算できるようにしたい(これができたら+30点)。

(M+++) 前橋氏の電卓でgcでメモリの確保をしたらどうなるだろう(これができたら+30点)。

(4) ファイルも読み込めるようにするというのは前橋氏の電卓はすでにできる。スクリプトと関数定義を読み込めるようにしてみよう。(+30点)。

(5) オブジェクトが定義でき,演算子多重定義ができるようにするのも面白い。Scilabのtyped listは演算子多重定義できるオブジェクトである(多分参考になるはず)。MATLABのクラス定義の仕方も変わっていて面白い(中がみれないので仕様しか参考にならない)。+50点

(6) 絵が書けるようになりたい。gnuplotを呼び出す仕組みにするのは簡単だろう。ただ,これだけでは寂しい。Javaで絵を描くのはいいかもしれない。+50点。

第5回解答解説

今回のレポートは大変であったが数値計算の講義の総集編となる課題であることは理解してくれている感想が多かった。250点以上のレポートを優秀レポートとして以下に3つのレポートを紹介する。

まず,河邊昌彦 君の作品である。これは上記方針から完全に自由な数値計算インタプリタの提案になっていて,素晴らしい。

次は,永田貴彦君の作品である。これは,上記方針をたどっていったときに,かなりがんばるとこうなるという例である。

最後は,永山陽平君の作品である。BLASは使っていないが,行列周りを頑張っている。

ここに掲載して,その努力(と才能)を称える。尚,今回は,私もこの課題に取り組んだ。現在のところ,ユーザがオブジェクトを定義できるようにする機能以外は実装されており,総ステップ数が1万を優に超えながらも快適に動作している。研究室紹介の際など適当な機会に紹介したい。尚,前橋和弥氏は氏のwebの参照やプログラムの利用を快諾して下さいました。ここに深く感謝の意を表します。ありがとうございました。また,bisonやflexの参考webを作成されている方にも謝意を捧げる。

尚,今回は最後のレポートであったので,感想をつけてくれた人が多かった。それをまとめてみた。文字化けが多いのは,sift-jisのままtthを使ったためで,勘弁してほしい。

感想

また,課題の(1+),(1+++++)に対する解答例を示そう。全般に永田貴彦君のレポートが参考になる。また,中島求君のレポートも面白い部分がある。

定期試験問題の解説

試験問題(定期試験) 教科書,ノート,電卓持込可

次の問題の空欄を埋めよ。尚、解答は空欄を埋めた部分しか見ないので、注意のこと。

問題1


cos0.0001 - cos0.00011
の値をXとする。Xの値を有効数字3桁まで正しく求めると、


X=
となる。

問題2

非線形連立方程式


x2 + 4y2 + 4xy - 2x - y
=
0
2x2 + 3xy - 2y2 -1
=
0
に対して、初期値を(x, y) = (1, 1)としてNewton法による反復を 一回行って得られる値を有効数字3桁まで書くと


x=
となる。

問題3

1階のスカラー変数常微分方程式


dx(t)
dt
= 1-x2(t) + cost,  x(0)=1
(1)
の解x(t)の2回導関数のt=0での値x(0)を有効数字3桁まで正しく求めると


x(0)=

となる。

解答

問題1

これは,桁落ちを避けるように,数式の変形をすることによって解く。加法定理により,sinに直して値を計算する。sinの値はsinのテイラー展開を利用して計算する。

答え 1.04 x 10^(-9) または1.05 x 10^(-9)

問題2

これは,ヤコビ行列を計算するところが山。途中,

[1;1]-inv([4 11;7 -1])*[6;2]

などの計算を経て(これはMATLABの記号,実際はこれを手と電卓で計算する),

答え 0.654 または 0.655

問題3

アダムスーバッシュフォースの公式をそのまま適用すればよい。

答え -2 または -2.00

 


©大石進一

このページのURIはhttp://www.oishi.info.waseda.ac.jp/~oishi/lec2001/num.htm

最終更新:2001/9/13 | トップページ