今後の実習では、処理を順次記述したプログラムをテキストファイルとして作成し、それを読み込むことによって処理を行うという手順を繰り返す。テキストファイルの作成には一般にテキストエディタとよぶソフトウエアを利用するが、この授業ではoctaveに組み込みのエディタを用いる。
Newton法は解析的に求まらない方程式の解を数値的に求める際の基本的な解法であり、これをoctaveでプログラムする。
$f\in C^2({\bf R })$は$u\in{\bf R}$で$f(u)=0$, $f''(u)>0$を満たすとする。点列$\{x_n\}$を次で定義する。\[ x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)} \]このとき、ある$\varepsilon$が存在し、$u$の$\varepsilon$近傍$U_\varepsilon(u)$の任意の点を初期値$x_0$にとると、$\{x_n\}$は$u$に収束する。
これは、$x_n$で引いた$y=f(x)$の接線と$x$軸との交点を$x_{n+1}$にとるという意味であり、定理の仮定からこれが解に収束するのはほぼ明らかである。
収束の程度を評価したい。$|x_{n+1}-x_n|$の関係を調べる。素直に変形し て \[ f(x_n)=f(x_{n-1}-(x_{n-1}-x_n))=f(x_{n-1})-f'(x_{n-1})(x_{n-1}-x_n)+\frac{f''(x_{n-1})}{2}(x_{n-1}-x_n)^2+o(|x_{n-1}-x_n|^2) \]となる。 \( x_{n}-x_{n+1} = f(x_n)/f'(x_n)\)から 左辺について\(f(x_n)=f'(x_n)(x_{n}-x_{n+1})\)である。 右辺に$x_{n-1}-x_n=f(x_{n-1})/f'(x_{n-1})$を代入すると$f(x_{n-1})$は キャンセルし、 \[ f'(x_n)(x_{n}-x_{n+1}) = \frac{f''(x_{n-1})}{2}|x_{n-1}-x_n|^2+o(|x_{n-1}-x_n|^2) \] \[ |x_{n+1}-x_n| = \frac{f''(x_{n-1})}{2f'(x_n)}|x_{n}-x_{n-1}|^2+o(|x_{n}-x_{n-1}|^2) \] を得る。 これを二次の収束と呼び、極めて速い収束である。
Newton法は与えられた方程式の解を具体的に求める手段である。このように、問題に対して具体的な解を与える手続きをアルゴリズムと呼ぶ。プログラムを作成する際にはアルゴリズムを決めなければならない。アルゴリズムの効率によって、作成したプログラムの良し悪しが決まる。
Newton法に従って$x^2-2=0$の正値解を求めるプログラムを組むと次のようになる。これを適当なファイル名でテキストファイルとして作成する。ただし、octaveのプログラムは.mで終わる名前をつける習慣がある。
#から改行まではプログラムの実行には影響しない。このような部分を注釈と呼び、プログラムを読みやすくするために使う。
x = 10.0 # 初期値の設定
eps = 0.00001; # 許容誤差の設定
err = 1; # 誤差の設定
while( err > eps )
y = x - (x*x-2.0)/(2.0*x); # Newton 法の漸化式
err = abs(x - y); # 現在の誤差
x = y
endwhile
画面への出力は、変数xへの代入文にセミコロンを付けないことで代用している。
ここではnewton1.mというファイル名にした。octaveからsource( "filename" )で読み込み、同時に実行する。
% octave
GNU Octave, version 2.0.17 (i386-vine-linux-gnu).
Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 John W. \
Eaton.
This is free software with ABSOLUTELY NO WARRANTY.
For details, type `warranty'.
octave:1> source( "newton1.m" )
x = 5.1000
x = 2.7461
x = 1.7372
x = 1.4442
x = 1.4145
x = 1.4142
x = 1.4142
octave:2>
誤差はどれだけ小さくできるだろうか。それには計算機の内部で数値がどのように扱われているかを知る必要がある。以前のプリントに書いた通り、計算機イプシロンよりも小さい数は扱えない。
先のプログラムは正しく動くが、機能の構成が一見してはわからない。Newton法を使っていることを明確にするため、関数を新たに定義して使う。次のようなプログラムになる。
function y=f(x)
y=x*x-2.0;
endfunction
function y=df(x)
y=2.0*x;
endfunction
x=10.0; # 初期値の設定
eps = 0.00001; # 許容誤差の設定
err=1; # 誤差の設定
while( err > eps )
y = x - f(x)/df(x); # Newton 法の漸化式
err = abs(x-y); # 現在の誤差
x = y
endwhile
f(x)とその微分を独立した関数として定義した。関数の内部で使用する変数は局所変数と呼び、関数の外部で同じ変数名を使っても独立に使うことができる。関数を定義することでプログラムの構造がわかりやすくなる。プログラムの読みやすさは何よりも重要な点の一つである。