23.2.14. 科学技術計算の簡単な例

Excel を用いると、科学技術の分野で出てくる数値計算を行うことができます。Excel では高度な計算や精密な誤差評価まではできませんが、ちょこっと実験的に計算をするだけなら簡単です。

今回は例として √2 = 1.41421356… の近似計算の実験を、Excel でやってみたいと思います。√2 は、放物線 y = x2 – 2 と x 軸の交点のうち x>0 の領域にあるものです。ですから √2 の計算は「関数 f(x) = x2 – 2 に対し、f(x) = 0 を満たす点 x を求める」という問題に言い換えられます。このような方程式の解を数値的に求める手法として、二分法と Newton 法というものが良く知られています。

二分法 #

二分法は、ある関数の値が 0 になる場所を、区間幅を半分に狭めながら探っていく手法です。

二分法の原理 #

関数 f(x) = x2 – 2 は区間 [0, 2] で短調増加です。さらに f(0) = -2 と f(2) = 2を合わせると、区間 [0, 2] 上に f(x) = 0 となる点が唯一存在することが分かります。これを使うと「区間を左半分と右半分とに分けて、どっちに f(x) = 0 が含まれるかを調べる」という操作を十分な回数繰り返すことで、f(x) = 0 となる x の範囲を狭めていくことができます。いまの場合

区間 [0, 2] の中点では f(1) = -1 < 0 なので、f(x) = 0 の解は [1, 2] 上のどこかにある
区間 [1, 2] の中点では f(3/2) = 1/4 > 0 なので、f(x) = 0 の解は [1, 3/2] 上のどこかにある
……

といった操作を繰り返すことになります。このステップを 1 回行うごとに解の存在し得る範囲を半分に狭められますから、n ステップ進めることで区間の幅を 2-n にできます。

実際の計算 #

それでは実際に、二分法で √2 の計算をしてみましょう。今回は計算の技法を理解することが目的ですから、予めセル A1 に「真の √2」、セル A2 に 1.41421356 を書き込んでおきましょう。その上でセル A3, B3, C3, D3 の順にそれぞれ「ステップ数」「区間の左端」「区間の右端」「区間の中点」と書き込みます。以下、この下に計算の情報を書き込むことになります。

まずは初期値の設定です。最初のステップ数は 0 なので、セル A4 に 0 を入力します。また最初は区間を [0, 2] に取るので、区間の左端 B4 には 0, 区間の右端 C4 には 2 を設定します。そして最後に、区間の中点 D4 に式

=(B4+C4)/2

を書きます。ここまでの状況は、次の図のようになっているはずです。

続いて、1 ステップ進めた後の計算を書きましょう。セル A5 には 1 を書き込みます。そして前のステップにおける区間の中点で x2 – 2 の値を見て

符号が負だったら、次のステップにおける区間の左端は、前のステップにおける区間の中点へと更新する
符号が正だったら、次のステップにおける区間の左端は、前のステップにおける区間の左端から動かさない

というようにします。これは Excel の IF 関数を用いると

=IF(D4*D4-2<0,D4,B4)

と表せます。これをセル B5 に書き込みましょう。

区間の右端も、同様に考えます。右端の場合は、前のステップにおける区間の中点での x2 – 2 の値を見て

符号が正だったら、次のステップにおける区間の右端は、前のステップにおける区間の中点へと更新する
符号が負だったら、次のステップにおける区間の右端は、前のステップにおける区間の右端から動かさない

という操作をします。ですから Excel の式で書けば

=IF(D4*D4-2>0,D4,C4)

となります。これをセル C5 に書き込みましょう。

区間の中点 D5 の式は D4 の行番号をずらすだけなので、D4 の式をコピー & ペーストすれば埋まります。ここまでできれば、次のステップ移行もやることは全く同じです。したがって行全体を下方向にコピー & ペーストすれば、自動的に行番号が 1 つずつズレた式が入り、計算が行われます。

ちゃんと √2 の値に近づいていることを、確認してください。

計算の後は、誤差の評価をしましょう。たとえば 0 ステップ目での誤差は、D4 – $B$4 です。そこでセル E3 に「誤差」と書き、その下のセル E4 に

=ABS(D4-$B$1)

と書きます。ここで、実際に知りたいのは「誤差の大きさ」だけですから、絶対値を意味する ABS 関数を施しておきました。

あとはオートフィルなりコピー & ペーストなりで、セル E4 の中身を下の方まで埋めましょう。そうすると、各ステップにおける誤差が数値で現れます。

ここで現れた E-05 や E-06 は指数表記の略記法で、それぞれ 10-5, 10-6 を表します。0.00… のように 1 の位から書き始めると値が長くなりすぎるので、このように書かれます。これで確かに、ステップを進めるごとに誤差が 0 へと近づいていることが分かりますね。

Newton 法 #

Newton 法は、接線を使って f(x) = 0 の解の場所を調べる方法です。

Newton 法の原理 #

既に述べたように、関数 f(x) = x2 – 2 は区間 [0, 2] 上で 1 回だけ値が 0 になります。そして f”(x) = 2 なので、y = f(x) のグラフは下に凸です。そこで x = 2 を出発点として「接線を引き、それの x 切片を求める」という操作を繰り返すことで、f(x) = 0 となる x の値に近づくことができます。

今の場合、放物線 y = x2 – 2 の x = t における接線は、微分を用いると y = 2tx – t2 – 2 と求まります。したがって

x = 2 における接線 y = 4x – 6 の x 切片は 3/2 である
x = 3/2 における接線 y = 3x – 17/4 の x 切片は 17/12 である
……

という操作を繰り返していくことになります。

実際の計算 #

この計算もやってみましょう。後で二分法と Newton 法の性能を比較したいので、二分法の計算をしたのと同じシート上で Newton 法の計算を行いましょう。

セル G3, H3 にそれぞれ「ステップ数」「値」という見出しを作ります。その後初期値として、セル G4 に 0 を、セル H4 に 2 を書き込みます。

続いて、ステップを進めるにはどうすれば良いか考えましょう。放物線 y = x2 – 2 の x = t における接線は y = 2tx – t2 – 2 で与えられます。そしてこの直線の x 切片は、y = 0 とおくことで x = (t2 + 2) / (2t) と求まります。ですから t に対して (t2 + 2) / (2t) を対応させることが、1 ステップ進める計算に相当します。これを行うため、セル G5 に 1 を入力した後、セル H5 に