情報地球科学演習 (2024年前学期) 第10回

今日のテーマ: 非線形方程式の数値解法

今回から演習の後半戦。後半戦では、物理現象の簡単な数値シミュレーションの入口にあたる実習を行う。その手始めとして、手計算で解くのが大変な 非線形な方程式を、コンピュータに代わりに解いてもらう ことをやってみる。

今日の演習では、2つの方法を試してみる。そのいずれも「一発で」解を求めるものではなく、目星をうまくつけながら少しずつ (「逐次的に」という) 求める解に近づいていくための方法である。

解法その1: 二分法

簡単な例題として、$ f(x)=x^2-2=0 $ という方程式の2つの解のうち、正のほう ($x=\sqrt{2}$) を二分法により求めるプログラムを作成してみよう。必要な手順 (の大部分) を以下に示す (別の画面で開きたい人はここを)。ただし 肝心なところは (わざと) 抜いている ので、各自で考えて埋めてほしい。

なお、このプログラムのミソ (というか動作の流れ) は以下の通り。

必要な「器」を宣言する。この例では、整数を入れる変数 i (この用途は後で説明する) の他に、実数を入れる変数として、
- 解くべき変数の値を入れる x
- $ f(x)=x^2-2 $ の値を入れる y
- x が存在しうる範囲の最小値 a と最大値 b
- 繰り返し計算を終了するかどうかの判定に使う小さな数 verysmall
- x=a や x=b のときの y の値を入れておくのに使う ya と yb
を用意している。
二分法による繰り返し計算を始める前に、解 x を探す範囲の最小値 a と最大値 b を指定する。講義のノートを参考にして、a と b の値を各自でうまく設定 しておくこと。とはいうものの、「絶対にダメ」な選び方をしない限り、二分法は正しい答を返してくれる ので、さほど深刻に考える必要はないが、a と b はいずれも正 にとっておいてください。
繰り返し計算に本当に突入する前に、念のため a と b の値が適切に設定されているかを調べる。 x=a のときの y の値 ya と、 x=b のときの y の値 yb とを考え、ya と yb の積が正の場合は a と b の設定が不適切 だと判定して処理を中止するようにしているのだが、余力のある人はその根拠を考えてみよ。
繰り返し計算を終了するかどうかの判定に使う verysmall を 1.0e-2 ($1.0\times10^{-2}=0.01$のこと) と指定している。同様に 1.0e-6 とすると $1.0\times10^{-6}$ の意味になる。
二分法による繰り返し計算を do ループにより行わせる。 do ループにはカウンタ変数 i を用いて、繰り返し回数の最大値を設定している (ここでは99回)。この理由は、(プログラムの設計ミスなどの仕業で計算の終了条件が満足されないことにより) 繰り返し計算が無限に続くことを防止したい からである。
二分法のカギの1つである、いま考えている範囲のちょうど中間の x の値での、関数 $ y=f(x) $ の値 を調べる。 この y の符号によって、次の繰り返し計算で解を探す範囲の最大値または最小値を変更 することが、二分法のもう1つのカギである。講義のノートを参考にして、以下の サンプル中のの部分に a または b を正しく書き入れ よ。なおここの設定で間違うと、期待している答に全く近づかない ので注意せよ。
二分法の繰り返し計算を終了する条件として、解を探す範囲の幅 b-a が十分小さくなったか かどうかをチェックする。 verysmall は「探す範囲がこれより小さくなったら止めろ」と判定するために用意したものである。
繰り返し計算の各回で、計算されている値を表示するようにしている。なおここでは、出力の書式仕様を write(*,*) とすることで、「小数点以下の表示に何桁使うか?」などを gfortran におまかせ にしている。

(

)

(

)

(

)

(

)

(

)

(

)

(

)

(

)

これを基にして、$ x^2-2=0 $ を満たす x の値を二分法により求めるプログラムを作成せよ [みほん]。このプログラムを (例えば) enshu11a.f90 という名前で保存し、例の如く gfortran enshu11a.f90 とした後に、

$ ./a.out

と実行し、期待通りの答 (もちろん「ひとよひとよに....」のアレです) が得られているか確認せよ。またこのプログラムの出力を注意深く観察し、二分法の反復のたびに 実にいい感じで a と b が近づいていく ことを確認せよ。なおこれ以降の演習では、これまでに使ってきた matsuyama-kion.txt などを使うことは (基本的には) ない。

プログラムが正しく動いていることが確認できたら、本日の演習の出席確認のため、このプログラムを末尾のフォーム経由で提出せよ。さらに (せっかくなので) 次の課題を行ってみよ。

上のプログラムで求まった答は、我々がよく知っている値と比べて、小数点以下の数字があまりよく合っていないはずだ。それを克服するには verysmall の値を変えて計算してやるとよい のだが、その理由は何か。また実際に verysmall の値をしかるべく変えて計算 を行い、より「真の値に近い」解を求めてみよ。ただし (計算機中の「有効数字」の問題のため)、verysmall はだいたい $10^{-6}$ 程度より小さくとっても意味がない。さらにその場合、verysmall のとり方によって、繰り返し計算の回数がどう変化したか も確認せよ。
繰り返し計算の最初に設定する、解 x を探す範囲の最小値 a と最大値 b を変化 させて計算を行ってみよ。その場合でも、期待通りの答が出てくることを確認せよ。 (ただし a と b はいずれも正 にとっておいてください。) さらにその場合、a や b のとり方によって、繰り返し計算の回数がどう変化したか も確認せよ。
(ここまでで終わると面白くないので)、期待通りの解が出てこない a と b の組み合わせを少なくとも1つ 見つけ、その場合になぜうまくいかないかの理由を考えよ。

余力のある人向けのおまけ: 二分法をもうちょっと高機能にしてみよう

1限の授業でも述べた通り、二分法が正しく動くためには、区間 $ a\le{x}\le{b}$ で関数 $f(x)$ が単調 (単調増加あるいは単調減少) であることが条件の1つであった (もう1つの条件は、$f(x=a)$ と $f(x=b)$ が異符号であること)。一方、既に気付いている人もいるであろうが、上で作成した二分法のプログラムでは、区間 $ a\le{x}\le{b}$ で関数 $f(x)$ が単調増加 (要するに「$y=f(x)$ のグラフが右上がり」) であることを前提にしている。では関数 $f(x)$ が単調減少な区間に存在する解を二分法で求めるには、上のプログラムのどこをどのように直せばよいかを考えよ。これに基づき、$f(x)=x^2-2=0$ の解として $x=-\sqrt{2}$ が求まるような二分法のプログラムを作成してみよ。

さらに余力のある人は、$f(x)$ が単調増加だろうが単調減少だろうが、とにかく単調であれば使える ような二分法の仕掛けを考えてみよ。それには上のプログラム例における、ya や yb に相当するものを積極的に利用するとよい [みほん]。具体的には以下の例のように、ある繰り返しの時点で ya と y の積が正であったならば、次の繰り返し計算で解を探す範囲をどう狭めることができるかを考えてみよ。

(

)

(

)

(

)

(

)

(

)

(

)

解法その2: ニュートン法

次に、同じ $ f(x)=x^2-2=0 $ という方程式の解をニュートン法により求めるプログラムを作成してみよう。必要な手順 (の大部分) を以下に示す (別の画面で開きたい人はここを)。ただし今回も 肝心なところは (わざと) 抜いている ので、各自で考えて埋めてほしい。

なおこのプログラムで、特に注意すべきミソは以下の通り。それ以外の箇所はだいたい、二分法のプログラムと同じである。

必要な「器」を宣言する。この例では、二分法でも使っている変数 i、x、y、verysmall の他に、実数を入れる変数として、
- 関数 $ f(x) $ の微分 $ f\prime(x) $ を入れる dydx
- 繰り返し計算の途中で得られる新しい推定値を入れておくのに使う x_new
を用意している。
ニュートン法による繰り返し計算を始める前に、解 x の最初の推定値を x に入れておく。講義のノートを参考にして、最初に入れる x の値を各自でうまく設定 しておくこと。とはいうものの、あまり近すぎると面白くないし、あまり遠すぎるのも困るので、真の解にそこそこ近い値 をうまく入れてみてください。
ニュートン法による繰り返し計算を do ループにより行わせる。 いま考えている x の値での、関数 $ y=f(x) $ の値とそこでの微分 $ f\prime(x) $ の値を計算 しておき、さらに これらを用いて、x の値を新しく推定し直す ことが、ニュートン法のカギである。講義のノートを参考にして、必要な式を正しく書き入れよ。なおその際、$2x$ は 2*x ではなく 2.0*x と書くほうが安全である (その理由は第4回で紹介済み)。
ニュートン法の繰り返し計算を終了する条件として、繰り返し計算で x を変更した量 |x_new-x| が十分小さくなったか どうかをチェックする。 verysmall は「x を変更した量がこれより小さくなったら止めろ」と判定するために用意したものである。また abs は絶対値をとるための組み込み関数である。詳細はここを参照のこと。
繰り返し計算の do ループの最後で、x に x_new の値を入れ直し、次の繰り返し計算に備える。なおこれを忘れると、いくらやっても期待している答に全く近づかない ので注意せよ。

(

)

(

)

(

)

これを基にして、$ x^2-2=0 $ を満たす x の値をニュートン法により求めるプログラムを作成せよ [みほん]。このプログラムを (例えば) enshu11b.f90 という名前で保存し、例の如く gfortran enshu11b.f90 とした後に、

$ ./a.out

と実行し、期待通りの答が得られているか確認せよ。

プログラムが正しく動いていることが確認できたら、次の課題を行ってみよ。

ここでもやはり verysmall の設定がよくないおかげで、上のプログラムで求まった答は小数点以下の数字があまりよく合っていないはずだ。 verysmall の値をしかるべく変えて計算 を行い、より「真の値に近い」解を求めてみよ。さらにその場合、verysmall のとり方によって、繰り返し計算の回数がどう変化したか、また 繰り返し計算の回数は二分法のときと比べてどうなったか も確認せよ。
関数電卓でもそう簡単に教えてくれないであろう、2の3乗根、4乗根、5乗根.... を求めてみよ。
x に関する次の方程式の解を求めよ。ここで $\log$ は自然対数である (底が10の常用対数ではない)。 \[ \log(x)-1=0 \] この方程式の解は指数関数で出てくる e (「ネイピア数」という) であるが、二分法またはニュートン法により e のおおよその値を求めてみよ。ちなみにネイピア数の値をプログラム中で参照するには exp(1.0) などとすればOK。さらに (せっかくなので)、期待通りの解が出てこない設定 (二分法ならば a と b の組み合わせ、ニュートン法ならば最初に入れる x の値) を少なくとも1つ見つけ、その場合になぜうまくいかないかの理由を考えよ。

上の例からも分かるように、ニュートン法には 正解から余りにもかけ離れた初期値からスタートしてしまうと、正解にたどり着かないことがある という難点がある。ということで、$y=f(x)$ のグラフをイメージしておくなどして、どのあたりに求めたい答えがありそうか をあらかじめ知っておくことが重要である。

二分法とニュートン法の「速さ」の比較

上でみたように、二分法とニュートン法では、同じ方程式を解かせる場合でも、「十分もっともらしい答」に至るまでの計算回数が異なっている。

ニュートン法が最終的な答に近づく直前のあたりでは、繰り返し計算のたびに 真の答からのずれが前回の値のほぼ2乗で小さくなっている ことが見てとれるだろう。このような収束のしかたを 「2次収束」 という。おおざっぱに言うと、「2次収束」する場合には、繰り返し計算のたびに正しい桁数が倍に増える ことになる。

これに対して二分法では、たとえ最終的な答に近づく直前のあたりをとったとしても、真の答からのずれが前回の値の数分の一程度にしか小さくなっていない ことが分かる。このような収束のしかたを 「1次収束」 という。

当然のことながら、「1次収束」の二分法と比べて、「2次収束」のニュートン法のほうが (もし収束するならば) 速く答が求まる。より高次の収束性をもつ方法を使えばさらに速く収束するはずではあるが、一般的には「2次収束」で十分実用的な計算スピードであるといえる。

まだまだ余力のある人向けのおまけ (その1): はさみうち法

今日の最初で紹介した二分法で $f(x)=0$ という方程式を解く際に、$x=a$ と $x=b$ の間に解が存在すると (あらかじめ) 分かった上で、次に $x=\displaystyle\frac{a+b}{2}$ での関数 $f(x)$ の符号を調べる、という手順が含まれていた。この手順には「解が存在しうる区間の幅を半分に縮める」ための情報を取得するという目的があったのだが、別の見方をすれば、$x=\displaystyle\frac{a+b}{2}$ という その区間のど真ん中の値が求めたい解になっているかどうか を調べている、と言えなくもない。この見方をもう一歩進めて二分法を考え直してみると、ただ単純に「その区間のど真ん中」を調べるよりも、その区間の中で より正解に近そうなところをどうにかして推定し、そこでの値を調べる、という作戦でもいいはずだし、場合によってはこっちのほうが効率的になりそうだ。そのような方法の1つとして、「はさみうち法」 (false position 法とか regula falsi 法とも) と呼ばれるものを紹介する。

二分法が適用可能な場合と同様に、区間 $a\le{x}\le{b}$ で関数 $f(x)$ が単調であり、かつ $f(x=a)$ と $f(x=b)$ の符号が異なっていると仮定する。この区間には $f(x)=0$ を満たす $x$ がただ1つ存在するので、その値を求めることを考えよう。その際、二分法では $x=\displaystyle\frac{a+b}{2}$ における関数 $f(x)$ の値 (というか符号だけ) を考えていたが、はさみうち法では \[ x=\frac{a\times{f(x=b)}-b\times{f(x=a)}}{f(x=b)-f(x=a)} \] という点における関数 $f(x)$ の値を求め、それが十分0に近いかどうかを調べていく。なおこの点は、$(a,f(a))$ と $(b,f(b))$ の2点を通る線分 \[ y=\frac{f(b)-f(a)}{b-a}(x-a)+f(a) \] とx軸 ($y=0$) との交点のx座標である。これにより、$f(a)$ と $f(b)$ の符号だけでなく、$f(a)$ と $f(b)$ のどっちがどれだけ0に近いか、という情報まで含めて、次に値を調べる点を探す ことができるようになる。

まだまだ余力のある人は、二分法のプログラムを改良することによって、はさみうち法のプログラムを作ってみよ [みほん]。同じ問題を二分法で解いた場合とはさみうち法で解いた場合では、どちらが速く答に到達したか。

多くの場合は二分法よりもはさみうち法のほうが速く答に到達しそうなものなのだが、関数 $f(x)$ の性質や、最初に指定する a と b の値の組み合わせによっては、二分法のほうが速いこともある。ということでここでもやっぱり「答のありそうな範囲をあらかじめ押さえておく」ことが大事、という結論になるのだった。

まだまだ余力のある人向けのおまけ (その2): 割線法

今日2つめに紹介したニュートン法は (うまくいく場合には) とても強力な方法ではあるのだが、x_new を計算するのに関数 $f(x)$ の微分 (=接線の傾き) を使う必要があるので、微分が簡単に計算できない (or 微分の計算なんてしたくない) 場合には不向きである。しかし、微分をどうにかして近似することができるなら、ニュートン法と同じような手順で x_new を計算することができるはずだ。

ここで紹介する「割線法」(セカント法とも) は、そのような作戦に基づいて考えられた方法であり、具体的には以下のような手順で計算していく。

解 x の最初の推定値として、$x_1$ と $x_0$ の 2つを用意し、くり返し計算を開始する。ただし二分法やはさみうち法とは異なり、解が $x_1$ と $x_0$ の間になくても (原理的には) OK。
$i$ 回目のくり返し計算 ($i=1,2,\cdots$) において、「近似的な」解を $x_i$ $x_{i-1}$ の2つとる。
次の式により、$i+1$ 回目の「近似的な」解 $x_{i+1}$ を計算する。 \[ x_{i+1}=x_{i}-f(x_{i})\dfrac{x_{i}-x_{i-1}}{f(x_{i})-f(x_{i-1})} \]
十分正確な解が求まったらOK、求まっていなければもう一度くり返す。

ニュートン法と割線法での手順の違いに注意せよ。ニュートン法では曲線 $y=f(x)$ 上の点 $(x_i,f(x_i))$ における接線を考えて x_new を求めていたが、その接線を 曲線 $y=f(x)$ 上の2点 $(x_i,f(x_i))$ と $(x_{i-1},f(x_{i-1}))$ を通る直線で代用 したのが割線法である。

まだまだ余力のある人は、ニュートン法のプログラムを改良することによって、割線法のプログラムを作ってみよ [みほん]。同じ問題をニュートン法で解いた場合と割線法で解いた場合では、どちらが速く答に到達したか。

出席確認用プログラム提出フォーム

Moodle のコースページ以下の課題ファイル提出により、各自で作成した (enshu11a.f90 に相当する) 二分法を用いて方程式 $f(x)=x^2-2=0$ の解のうち正のほう ($x=\sqrt{2}$) を求める Fortran 90 プログラムの最終版を提出せよ。