皆さん、こんにちは。
今回から、コンピュータシミュレーションに関する数学である「数値解析」について紹介していきます。
数値解析が対象とする分野は
・方程式の解を見つける →ニュートン法
・定積分を計算する →台形公式、シンプソンの公式 etc
・微分方程式を解く →ルンゲ・クッタ法、差分法による陽解法・陰解法 etc
・巨大な物理モデルを解く →有限要素法
などがありますが、最後の「有限要素法」の概要紹介を最終目的に、上から順番に紹介していく予定です。
初回の今回は、方程式の近似解を求める「ニュートン法」について紹介します。
0. 数値解析とは？
数値解析とは、数学的な問題を手計算ではなく「コンピュータで解く」ことです。
ここで、知らない人からすると衝撃的な事実を言います。それは、
「コンピュータは単純な四則演算しかできない」
ということです。
つまり、
・複雑な方程式を解く
といった手計算では当たり前にやってきた操作が、コンピュータにはできません。
これらの操作をコンピュータにさせようと思ったら、コンピュータが計算できるように「四則演算」の形に操作を書き換えることが不可欠になります。もちろん100%正確に置き換えることは大概の場合不可能なので、四則演算の形で近似する、というのが実態になります。
数学的な複雑な操作を、いかに精度よく、そして効率のいい四則演算に近似できるか？それを研究する学問が「数値解析」という分野なのです。
1. ニュートン法
数値解析とはどういうものなのか一番手っ取り早く理解できる例の１つが、今回紹介する「方程式f(x)=0を解く」という作業です。
（※今回求めていく解は「実数解」のみで、虚数解は対象外とします）
もしf(x)が単純な一次方程式なら単純な四則演算だけで解を求めることができますが、f(x)が2次式や3次式といった多項式、果ては指数関数、三角関数、対数関数、それらが全部混じった式ともなれば、四則演算だけで解くことはできません。
なので、一般的に方程式f(x)=0を解こうと思ったら、四則演算で何かしら近似する方法を見つけないといけません。その最も単純でメジャーな方法が「ニュートン法」です。
まずは、解の値にある程度あたりを付けて、それに近い実数x0を好きに選びます。
f(x)の符号を調べるなどして、実数解αを不等式で評価することは大学入試でも頻出するテクニックなので、それを使ってあたりを付けるわけです。
次のステップからは、y=f(x)のグラフと睨めっこしながら読み進めると理解が速いと思います。
x座標がx0となるy=f(x)上の点を取り、その点における接線を求めるのがステップ2です。
次のステップ3では、ステップ2で調べた接線とx軸との交点x座標を調べ、その値をx1とするわけです。
次のステップ4では、x1に対してステップ2,3と同じ作業をしてx2を導出します。
あとは同じ手順を何度も繰り返してx1, x2, x3,・・・・と数列xnを作っていきます。
上の図を見ると、xnに比べてxn+1の方がαに近くなってることが分かりますので、この操作を繰り返せば繰り返すほどxnの値はαに近づいていき、n→∞ではαに収束するはずです。
もちろん、本当に無限回計算することなどできないので、途中で計算を打ち切りそのタイミングで求まったxnの値をαの近似値と見なすわけです。
これがニュートン法の発想です。
実際上は漸化式④を使って次から次へとxnの項を求めていく作業がニュートン法になります。
ニュートン法に限らずですが、このように漸化式を使って繰り返し計算をする作業がコンピュータの得意分野です。
例として、√2の近似値を求めてみましょう。
そのために
(A)を漸化式④に当てはめると、
という漸化式が出来ます。
1<√2<2が２乗することで容易にわかるので、今回は初期値x0として
を採用します。
この下でx4まで計算した時点で、
小数点以下8桁まで正確に√2を再現できていることが分かります。
この例のニュートン法をPythonコードで表現すると次のようになります。
ちなみに、初期値としてx0=-1を採用すると、
負の方の実数解-√2が求まります。
2. ニュートン法の近似精度
次にニュートン法で求まる近似解xnと真の解αとの差、つまり近似精度がどの程度かを評価したいと思います。
そのために、漸化式④の右辺をαの周りでテイラー展開したいと思います。
関数を多項式の形で書こう！ ～テイラー展開・ローラン展開～ - ちょぴん先生の数学部屋 (hatenablog.com)
f(x)とf'(x)をそれぞれテイラー展開すると、
となるので、④の右辺の形を
とすれば、
g(x)はx-αの２次式で近似できることが分かります。
ここでO(△^m)という記号が登場していますが、これはランダウの記号と呼ばれるもので、「△のm次以上の項だけが存在する多項式」という意味です。
通常△に入る数（今回の場合はx-α)は非常に小さい数なので、次数が上がれば上がるほどどんどん小さくなっていきます。なので、今回のような誤差を議論したい場合には△の次数の低い項だけが重要になり、次数の高い項は無視してよいわけです。
そんな無視してよい高次の項の部分はランダウの記号でまとめてしまい、最終的に無視して切り捨てます。
というわけで⑨においては、分子は2次まで近似して3次以降は切り捨て、分母は1次まで近似して2次以降は切り捨てる、という近似を行なっています。
⑨を利用すれば、
xnとαとの差|xn-α|は、１回の計算ごとに2乗される形でどんどん小さくなっていくことが分かります。
但し、⑩の係数の分母にf'(α)があり、これが0に近い（つまりy=f(x)の傾きが緩やか）と係数全体の値が大きくなってしまい、誤差が小さくなりにくく計算回数を増やしていかないと精度が良くならないことが分かります。
よって、xnの近似精度は
・x=αが極値に近い→αへの収束が遅い ⇒計算回数を増やさないと精度が上がらない
・x=αが極値から遠い→αへの収束が速い ⇒計算回数が少なくても精度が出やすい
と言えます。
また、初期値x0をαからあまりにも遠い値で取ってしまうと、途中の計算で極値を通ってしまって分母のf'(x)の値が0にいって全体が発散してしまうリスクを孕んでしまいます。
なので、少ない計算で十分な精度を確保したいと思ったら、初期値x0はなるべくαに近い値にすべきなのです。
次回は、定積分を数値解析で計算する「台形公式」「シンプソンの公式」について紹介します。
