ニュートン法の視覚化

この記事は, 旧ブログから移植された記事です. よって, その内容として, 旧ブログに依存した文脈が含まれている可能性があります. 予めご了承下さい.

久しぶりにまた1なにか d3.js で視覚化してみたくなったのだが, このエントリがポストされる次の日はアイザック・ニュートンの誕生日らしいので, 今回はニュートン法 (Newton Raphson 法) を視覚化してみることにした. 早速であるが以下がその成果物である2. f(x)=0f(x)=0 となる関数 f(x)f(x) とその導関数 f(x)f'(x) 及びニュートン法の初期値を受け付け, 実行をクリックすると関数とニュートン法の計算過程における接線がプロットされる. デフォルトでは, 2\sqrt{2} の値を計算するように設定してある. 入力された関数を元に数値微分をしても良かったのだが, なんとなく導関数を入力したかったので, そのようなことはしなかった.





これで終わってしまうと何とも寂しいので, 一応簡単にニュートン法について書く. ニュートン法はとても有名な方程式の近似解を求める方法の 1 つである. 連続的関数3 f(x)f(x)f(x)=0f(x)=0 となるような xx を求めるときに, 予め決めた, あるいは事前の計算で求まった切片 xnx_{n} における関数 f(x)f(x) への接線 f(xn)f'(x_{n}) の切片 xn+1x_{n+1} を用いて再帰的に f(x)=0f(x)=0 に近づけていくことで求める. 同手法は非常に単純ながらも効率的な近似解法であり, 同手法から発展されたいくつかの手法が考えられている. いまニュートン法の漸化式を導出することを考える. 計算で必要となるのは f(xn)f'(x_{n}) に対する切片 xn+1x_{n+1} であるから, いま三点 (xn,0),(xn,f(xn)),(xn+1,0)(x_{n},0), (x_{n},f(x_{n})), (x_{n+1}, 0) の成す直角三角形について考えると, f(xn)=(xn)xnxn+1f'(x_{n})=\frac{(x_{n})}{x_{n}-x_{n+1}} より xn+1=xnf(xn)f(xn)x_{n+1}=x_{n}-\frac{f(x_{n})}{f'(x_{n})} 例えば 2\sqrt{2} を例にとると, x2=2x22=0x^2=2\Leftrightarrow x^2-2=0 なので f(x)=x22,f(x)=2xf(x)=x^2-2,f'(x)=2x とおいて xnxn222xnx_{n}-\frac{x^2_{n}-2}{2x_{n}} を計算すればよい. なお f(xn)f'(x_{n})00 であると, ゼロ除算になってしまうため計算することができない. この事実は, 傾きのない場合にどちらに進んでいけば良いのかわからないという直感的な考えにおいても筋が通る.

また, f(x)f(x) の解が複数あるとき, 初期値によっては望まない解が導かれることがある. いま 2\sqrt{2} の正の解を得たいとき, 初期値を 1-1 で実行してしまうと, 得られる解は 2-\sqrt{2} となる. 直感的には決められた初期値の傾き f(x0)f'(x_0) によって近づいていく方向が定まってしまうからといえる. 従って同手法を適用する際は, できる限り求めたい解に近い初期値を設定するのが望ましい.

なおニュートン法は 1 変数関数のみならず多変数関数に対しても同様にして解を求めることができる. あまり厳密には書かないが, まずは簡単のために 2 変数の関数 f1(x,y),f2(x,y)f_{1}(x,y),f_{2}(x,y) を用いてそれを示すこととする. f1(x,y),f2(x,y)f_{1}(x,y),f_{2}(x,y) は曲面の定義そのものであり, この 2 つの曲面の交わる曲線を辿っていくことで解が求まる. 1 変数のニュートン法の場合と同様に漸化式を求めていけば, 多変数関数に対するニュートン法の漸化式も関数の値をその傾きで割る部分が出てくるが, いま変数は複数であるので, 各方向への微小変化に対する変化量を求める必要がある. これを求めるには全微分をすればよいので, 結局

(xn+1,yn+1)=(xn,yn)Tf(xn,yn)1(f1(xn,yn),f2(xn,yn))T (x_{n+1},y_{n+1})=(x_{n},y_{n})^T-{\partial f(x_{n},y_{n})}^{-1}(f_{1}(x_{n}, y_{n}),f_{2}(x_{n},y_{n}))^T

ここで

f(x,y):=(f1(x,y)xf1(x,y)yf2(x,y)xf2(x,y)y) \begin{aligned}\partial f(x,y):= \left(\begin{array}{cc} \frac{\partial f_{1}(x,y)}{\partial x} & \frac{\partial f_{1}(x,y)}{\partial y} \\ \frac{\partial f_{2}(x,y)}{\partial x} & \frac{\partial f_{2}(x,y)}{\partial y} \end{array}\right) \end{aligned}

なお f(x,y)\partial f(x,y) はヤコビ行列といわれる. 実際にコンピュータで計算する際には, f(xn,yn)1(f1(xn,yn),f2(xn,yn))T{\partial f(x_{n},y_{n})}^{-1}(f_{1}(x_{n}, y_{n}),f_{2}(x_{n},y_{n}))^T を求めるのは計算量と誤差の観点から見て困難なので, f(xn,yn)a=(f1(xn,yn),f2(xn,yn))T\partial f(x_{n},y_{n}){\boldsymbol a}=(f_{1}(x_{n}, y_{n}),f_{2}(x_{n},y_{n}))^T を LU 分解などで解き (xn,yn)Ta(x_n,y_n)^T-{\boldsymbol a} と解くことになる.

因みに, 上記で描画される接線は, 単純に ラグランジュの補完公式より yy1=y2y1x2x1(xx1)y-y_{1}=\frac{y_{2}-y_{1}}{x_{2}-x_{1}}(x-x_{1}) を用いて描いている. 具体的には接線の関数を g(x)g(x^{\star}) としたとき, 接点と y2=0y_{2}=0 であるときの 2 点 (x,f(x)),(xf(x)f(x),0)(x,f(x)),(x-\frac{f(x)}{f'(x)},0) を使って

g(x)=f(x)f(x)f(x)(xx+f(x)f(x))=f(x)(xx+f(x)f(x))=f(x)xf(x)x+f(x) \begin{aligned} g(x^{\star})&=&\frac{f(x)}{\frac{f(x)}{f'(x)}}(x^{\star}-x+\frac{f(x)}{f'(x)}) \\ &=&f'(x)(x^{\star}-x+\frac{f(x)}{f'(x)}) \\ &=&f'(x)x^{\star}-f'(x)x+f(x) \end{aligned}

と導ける.


  1. 以前のエントリ, ベジェ曲線では d3.js を用いて二次ベジェ曲線が描かれていく過程を書いた.↩︎

  2. 実装. ここで懺悔すると, 実はグラフの描画の実装についてはそこそこ手抜きをしている. 例えば解が第 1, 2 象限であるものの関数 f(x)f(x) の値の多くが第 3, 4 象限にあると接線が見えない. 勿論計算結果そのものは影響しない.↩︎

  3. 連続性に関する論法 \to ϵδ\epsilon-\delta 論法.↩︎