Powered by SmartDoc

Introduction

数値計算 西谷 滋人

数値計算，講義・演習の進め方

数値計算が対象とするのは，関数の解，積分，微分方程式，固有値問題などで，数学問題を解析的(analytical)ではなく，数値的(numerical)に解く手法の集大成．数値計算には，料理と同じで，正しい調理法(レシピ，recipe)がある．

物理(あるいは科学技術計算)のちょっとした問題を解くときの戦術は，

  1. 手持ちのツールで解いてみる，
  2. 解析的な解を捜す，
  3. 既存の数値計算のサブルーチンを写す，
  4. 使いやすいライブラリを捜す，
  5. 自分でサブルーチンを考える，

という順でおおむね進める．数値計算のサブルーチンを実際に自分で考えるということはほとんどない．しかし，例え3.や4.の既存のサブルーチンに頼るときにもブラックボックスの中で何がおこなわれているかを大まかに理解しておかないと大失敗をしでかすことになる．

といっても，数値計算の全てを理解しておく必要はない．料理において基本となる，切る，煮る，焼くなどを知っていれば，recipeを見ながら調理を進めていくことが可能である．数値計算においても同じで，基本となる誤差，精度，収束性，安定性などの本質を理解しておけば新しい数値計算手法もだいたいの振る舞いを予測することが出来る．

本講義ではこのような視点にたって，いくつかの典型的な数値計算手法の基礎的な考え方と実際のプログラムを紹介する．教科書の例は全てFortranで記述されている．数値計算の商用プログラムの多くがFortranで書かれているため，本格的に数値計算をするときには知っていると便利である．C言語に比べて予約語の数が少ないので，文法を憶えるのも楽であろう．しかし，情報学科の学生の理解を素早くするために，プログラム例はC言語に書き換えても紹介する．また，前述の戦術1.の手持ちのツールとして強力な"Maple"を講義では頻繁に使用する．Fortran, Mapleの理解に必要な文法はその都度大まかに紹介するが，使いこなすには自習が不可欠である．特にMapleは，演習室にインストールされているだけでなく，各自所有のパソコン(Win, Mac, Linux)にインストールできる．希望者は西谷(4階，教授室45)まで相談に来るように．

本講義で扱う題材は，

  1. 誤差
  2. 代数方程式(二分法，Newton-Raphson法，Bairstow-Hitchcock法)
  3. 行列，固有値(LU分解，Gauss消去法，Jacobi法)
  4. 関数近似(FFT,非線形最小二乗法)
  5. 微分方程式(Euler法，Runge-Kutta法他)

を予定している．

成績評価は2-3回のレポートと期末の試験でおこなう．同期に開講する演習では講義と関連する問題を実際に手を動かして解き，理解を深める．講義だけを受講する学生は十分に復習およびプログラミングの自習をおこなっておくように．

具体的な例

多倍長計算：グーグル人材募集広告

2004年の夏にちょっとした話題になったGoogleの求人試験には自然対数の底eを数百桁で求める必要がある．これは通常の計算機演算ではできない．これには多倍長演算という特殊な数値計算手法を用いる必要がある．

グーグル、謎の人材募集広告--シリコンバレーのビルボードに　

Stefanie Olsen (CNET News.com)

2004/07/12 08:40

先週、シリコンバレーの中心を走るハイウェー101沿いのビルボードに、複雑な数学の問題を載せた広告が現れた。(中略)

この広告には、「｛eの値中の、最初の連続する10桁の素数｝.com」("{first 10-digit prime found in consecutive digits e}.com." )と書かれている。この答えの「7427466391.com」にアクセスすると、そのウェブページにはさらに別の問題(下記参照)が用意されているが、ここにもGoogleが関与していることを示すものは全くない。

この問題を解くと、Googleの研究開発部門「Google Labs」へのページに辿りつく。このページには、「Googleの構築を通して我々が学んだことの1つに、自分が何かを探しているとき、向こうも自分を探している場合のほうが見つかりやすいということがある。我々が探しているのは、世界最高のエンジニアであり、あなたこそその人なのだ」と書かれている。

7427466391.com
Congratulations. You've made it to level 2. Go to 
www.Linux.org and enter Bobsyouruncle as the login 
and the answer to this equation as the password. 
f(1)=7182818284
f(2)=8182845904
f(3)=8747135266
f(4)=7427466391
f(5)=__________

行列：ラジオシティ法

グラフィックスで物体の表面の光の反射を計算する方法の一つに，ラジオシティ法がある．そこでは，シーンの中の全ての物体の表面を「パッチ」と呼ばれる小さな領域に分割して，線形連立方程式を作る．このパッチの数は非常に多いため，数万行の行列式を解く必要がある．このような単品種大量計算に対しては，C言語やFortranと高速ライブラリーを用いた，並列クラスター計算が必要となる．一方，Mapleは，いわゆる多品種少量計算を得意とする．

コンパイルと実行

コンパイル

授業ではMac上のターミナルを使って実際のプログラム手順などを紹介する．linux上でのprogrammingとほとんど変わらない．プログラムの作成・実行は，1)emacsによる編集，2)gccおよびg77でのコンパイル，3)./a.outなどによる実行，4)gnuplotなどによる視覚化という手順でおこなう．最も簡単なcおよびfortranプログラムのコンパイル，実行結果を示す．

C
[BobsNewPBG4:~/NumRecipe/chap2] bob% cat c21.c
#include <stdio.h>

int main(void){
  printf("Hello world!\n");
}
[BobsNewPBG4:~/NumRecipe/chap2] bob% gcc c21.c
[BobsNewPBG4:~/NumRecipe/chap2] bob% ./a.out
Hello world!
Fortran
[BobsNewPBG4:~/NumRecipe/chap2] bob% cat c21.f
      program main
      write(*,*) "Hello world!"
      end
[BobsNewPBG4:~/NumRecipe/chap2] bob% g77 c21.f
[BobsNewPBG4:~/NumRecipe/chap2] bob% ./a.out
 Hello world!

fortranはプログラム中で大文字小文字を区別しない．伝統的には全て大文字でコードは書かれている．本講義ではC言語的な小文字を多用した書き方をする．各行のあたまには通常6文字分の空白をあける．例外は注釈，継続行，およびラベルである(後に詳しく述べる)．write(*,*)は「書式指定なしでdefaultの出力にそれ以降の文字や数字を書きだせ」という命令である．

結果の出力

出力結果をファイルに保存するには

media: a.out > res1000.dat

等とすればよい．

参考文献

  1. 宮下精二著「数値計算」(朝倉書店，2002)

    物理の問題を実際に解くときに役立つ数値計算法が，ちょっとしたコツとともに紹介されている．特に量子力学，モンテカルロ法などの詳しい解説は，日本語では希少で役に立つ．数式，プログラムには誤植が多いので注意．

  2. 河村哲也著「数値計算の初歩！」(山海堂，2002年)

    レベルは高くないが，重要な数値計算の初歩を丁寧に解説．

  3. William H. Press他著「ニューメリカルレシピ・イン・シーC言語による数値計算レシピ」(技術評論社，1993)

    数値計算のバイブル．原著Numerical recipeでは非常に広範な計算対象に対して，C, Fortran, C++, Pascal, Basic版が用意されており，数値計算プログラムをコーディングする際の洗練されたスタイルも提示している．記述は初学者には難しいが，ある程度経験を積んだプログラマには，手法を選ぶうえで非常に役に立つ情報である．

  4. 奥村晴彦著「Ｃ言語による最新アルゴリズム事典」(技術評論社，1991)

    いわゆる数値計算に限らず，いろいろな計算機問題とその解法が載っている．数値計算についても必要最小限の記述とプログラムがまとめられており，非常に便利．