第 2  回Fortran  勉強会

辻野 智紀

  

本日のお題

前回の Advanced  の解答

高次代数方程式の数値計算 Newton  法

二分法

多元連立 1  次方程式の数値計算 ガウスの消去法

ピボットの選択

  

実習の前準備

以下のページから ,  データファイルをダウンロードしておいてください .

http://rain.hyarc.nagoya­u.ac.jp/~satoki/main/calc/Nagoya­Fortran/seminar­2.tar.gz

  

高次代数方程式の数値計算

  

高次代数方程式の数値計算

みなさん , 以下の方程式を計算してください .はたして , 解析解が求められるでしょうか？

  

高次代数方程式の数値計算

みなさん , 以下の方程式を計算してください .はたして , 解析解が求められるでしょうか？

X = 0 は自明ですよね .

しかし , 次の図のように , この方程式は x = 0 以外に , x = 5 付近にもう１つの解が存在することがわかり

ます .

  

この点は自明 この点を求めたい

  

高次代数方程式の数値計算

実は , 先ほどの方程式は ,ウィーンの変位則

の証明に必要な方程式です .

あの方程式を計算すると ,

が得られます .時間があれば先の方程式を導出してみましょう .

( プランク関数を波長で微分すればよいです )

  

ここでは , x = 0  でないものを求めたい .

そこで ,  数値計算を用いて ,  数値解を求めて見ましょう .

  

Newton  法

1. 方程式を f (x) = 0  という形に直す .

2. 任意の点 x = x (1)  における関数の値 f (x (1)) および ,  その点での接線の傾き f ' (x (1))  を計算する .

3. この f (x (1))  から接線を x  軸に伸ばし ,  その接線と x  軸との交点を x = x (2)  とする .

4. この作業の繰り返し .

5.  f (x(i)) << 1  となる点を数値解と近似する .

  

図でみる Newton  法

  

式でみる Newton  法

x (i) を求める漸化式は

解となるための条件は , 先のもの以外に

  

Newton  法の短所

Newton  法は ,  図で見たように ,  すぐに収束するアルゴリズムである .

一方 ,  関数 f (x)  の極値付近（ f '(x) = 0  ）で計算を行うと ,  不安定になる性質がある .

また , f (x)  とその微分 f '(x)  をあらかじめ計算しておかなければならない .

  

そこで , Newton  法より極値の扱いが上手な（？）アルゴリズムを考える .

  

二分法

1. f (x) = 0  という形の方程式に直す .

2. 解となる値 x =  α の前後で任意の値 x = a, b を決める .

3. この 2  点間の中点を x = c  とする .

4. f (a) * f (c) < 0  なら , b  に c  を代入 . f (a) * f (c) > 0 なら , a  に c  を代入 .

5. 以上の繰り返し .

6. 解の条件は , Newton  法と同じものを使用 .  もしくは , f (c) << 1  であればよいとする .

  

図でみる二分法

  

Newton  法と二分法

2  つの計算手法を見たとおり , Newton  法のほうが ,  計算の手順が複雑ですが ,  すみやかに収束するので ,  計算時間は短いです .

一方 ,  二分法は計算時間は多少かかりますが , 手順が簡単で , Newton  法が苦手とする ,  極値付近での計算も得意とするため ,  こちらのほうが便利だと思います .

  

実習 1:  超越方程式の数値計算

先ほどの ,  ウィーンの変位則で必要な代数方程式（超越方程式）を Newton  法 ,  二分法のいずれかで計算して ,  数値解を導出してみましょう .

時間があれば ,  両方作成し , 1  ステップで厳密解への収束がどれほどであるかを計算して ,  グラフに描画してみてもよいかもしれません .

詳しくは ,  配布資料参照 .

  

実習 1:  超越方程式の数値計算

ちなみに , 実際これらの計算手法を比べてみると , 収束の度合いは , 一目瞭然！

  

実習 1:  超越方程式の数値計算

  

実習 1:  超越方程式の数値計算

収束判定条件は , x  を解の近似としたとき , ABS (F(x)) < 10^(­6)  となったとき .

開始位置は ,  二分法が , a = 3.0, b = 6.0 . Newton 法が x0 = 6.0  である .

横軸は ,  計算の回数 .  縦軸は , ABS (F (x))  の変化である .

明らかに ,  二分法の方が ,  計算回数が多いことがわかる .

  

作成にあたっての知識 ( 文関数 )

プログラムのはじめの ,  変数の定義の直後に ,  そのプログラムで使用する専用の関数を定義することができる .

宣言 (example) real F, x F (x) = x ** 2

  

Appendix:  ランベルトの W  関数

先ほどウィーンの変移則の方程式は , 解析解を解くのが難しいといいましたが ,

形式解は存在します .

これが , 先の方程式の形式解です .

この関数には , 名前があり ,ランベルトの W 関数

といいます .

  

参考文献

会田勝『大気と放射過程』東京堂出版 東京大学理学部地球物理学科

『数値計算入門』

Wikipedia http://ja.wikipedia.org/wiki/

%E3%83%A9%E3%83%B3%E3%83%99%E3%83%AB%E3%83%88%E3%81%AEW%E9%96%A2%E6%95%B0

http://ja.wikipedia.org/wiki/NaN