Newton 法(非線形方程式の求解)
スカラーの非線形問題
Newton法の理屈は授業で教えよう.理屈は大変に簡単だ.
そして今回は、その理屈を「実際に動かせるか」というところを頑張ってみよう.
例えば, $x^2 = 2$ を解いて解 $\sqrt{2}$ を求めてみよう.
典型的な手法は、まずは問題を 左辺 = 0 の形に直すところから始まる.
|
|
さらに、Newton 法にはこの関数の微分形も必要だ. 今回は幸い手で求まるので、与えてしまおう.
|
|
さて実際に Newton 法を実行してみよう.Newton 法は近似解がだんだん良くなっていく(と期待する)反復解法なので、初期値を設定してから、ループを回すという形になる.
|
|
1: 1.0, -1.0 2: 1.5, 0.25 3: 1.4166666666666667, 0.006944444444444642 4: 1.4142156862745099, 6.007304882871267e-6 5: 1.4142135623746899, 4.510614104447086e-12
どうやら Newton 法による改善が 4回行われ、5回目のループに入っての判定で終了となったようだ. この場合は、$x^2 - 2$ が大変小さくなっているようだ($10^{-12}$ぐらい).
念のために、$\sqrt 2$ と x を比べてみよう.
|
|
(1.4142135623730951, 1.4142135623746899)
大変近いな.実際の差も見ておこう.
|
|
-1.5947243525715749e-12
関数だけでなく、近似値の値そのものも 12桁ぐらいあっていることがわかる. この計算はうまくいったと言ってよいだろう.
誤差 abs$(\sqrt{2} - x)$ がどのように変化していくか、グラフにしてみよう. まずはいつもの準備をして…
|
|
上で作った sq
ベクトルを使って次のようにすればよいだろう.
|
|
誤差が急速に小さくなっていることが分かる.
さらに、対数グラフにしてみようか.
|
|
対数グラフで見ると、誤差が急速に小さくなっていることがよりはっきりと見て取れる.たいしたものだ.
欲張って、20桁の精度が出るかやってみよう
|
|
1: 1.0, -1.0 2: 1.5, 0.25 3: 1.4166666666666667, 0.006944444444444642 4: 1.4142156862745099, 6.007304882871267e-6 5: 1.4142135623746899, 4.510614104447086e-12 6: 1.4142135623730951, 4.440892098500626e-16 7: 1.414213562373095, -4.440892098500626e-16 …略… 99: 1.414213562373095, -4.440892098500626e-16 100: 1.4142135623730951, 4.440892098500626e-16
どうやら実質的に無限ループに陥ってしまったようだ.
これは、(特に指定しなければ)有効桁数が約 16桁のシステム(float64)が使われるからだ.この制限は Julia に限らず多くの環境で共通だ.これをなんとかするのには少し知識が必要なので、またいずれにしておこう.
ここで重要なのは、多くの計算機・プログラム環境ではそもそも約16桁以上の精度を求めようとするのはそのままでは無理ゲーだ、ということを知っておくことだ.何桁の精度が欲しいのか、常に意識しておくと健全だ.
非線形「連立」方程式の場合、近似解を Newton 法で求めるには?
例題として、
\[ \left\{ \begin{array}{ccc} x + 2y & = & -1, \cr x^2 + 2y^2 & = & 1 \end{array} \right. \]
の近似解を求めてみよう.
考え方はこれまでと同じだ.Taylor 展開を1次まで行った近似式を使って関数を近似する、という格好だ. やってみよう.
まず、上の式を $f(x) = 0$ の形に整理することを考えよう. それには、ベクトル値 $\boldsymbol{x}$ を入力にもらって出力にベクトル値を返す関数 $\boldsymbol{f}$ として
\[ \boldsymbol{f}(\boldsymbol{x}) = \left( \begin{array}{c} f_1(x,y) \cr f_2(x,y) \end{array} \right) := \left( \begin{array}{c} x + 2y + 1 \cr x^2 + 2y^2 - 1 \end{array} \right) \]
を定義して、与えられた問題を
\[ \boldsymbol{f}(\boldsymbol{x}) = \boldsymbol{0} \]
を満たす $\boldsymbol{x}$ を求めること、とすれば良い.
今の近似解の値を
\[ \boldsymbol{x}_{{old}} = \left( \begin{array}{c} x \cr y \end{array} \right) \]
とし, これから求める新しい近似解の値を
\[ \boldsymbol{x}_{{new}} = \boldsymbol{x}_{{old}} + \Delta \boldsymbol{x} = \left( \begin{array}{c} x + \Delta x \cr y + \Delta y \end{array} \right) \]
として、新しい近似解で
\[ \boldsymbol{f}(\boldsymbol{x}_{{new}}) \cong \boldsymbol{0} \]
が満たされていると期待しよう. このとき、Taylor 展開によりこの式の左辺は次のようになる.
\[ \boldsymbol{f}(\boldsymbol{x}_{{new}}) = \boldsymbol{f}(\boldsymbol{x}_{{old}} + \Delta \boldsymbol{x}) \cong \boldsymbol{f}(\boldsymbol{x}_{{old}}) + J(\boldsymbol{x}_{{old}}) \Delta \boldsymbol{x}, \]
ただし
\[ J = \left( \begin{array}{cc} \partial f_1/\partial x & \partial f_1/\partial y \cr \partial f_2/\partial x & \partial f_2/\partial y \end{array} \right) = \left( \begin{array}{cc} 1 & 2 \cr 2x & 4y \end{array} \right) \]
である.
だから、この二式をあわせると全体としては、既知の値 $x,y$ と未知の値 $\Delta x, \Delta y$ に対して
\[ J(\boldsymbol{x}_{{old}}) \Delta \boldsymbol{x} \cong - \boldsymbol{f}(\boldsymbol{x}_{{old}}) \]
この場合はすなわち,
\[ \left( \begin{array}{cc} 1 & 2 \cr 2x & 4y \end{array} \right) \left( \begin{array}{c} \Delta x \cr \Delta y \end{array} \right) \cong - \left( \begin{array}{c} x + 2y +1 \cr x^2 + 2y^2 -1 \end{array} \right) \]
という連立一次方程式の形になり、 これを解けば未知の値 $\Delta x, \Delta y$ が求まり、 「より良いと期待される」近似値 $\boldsymbol{x}_{{new}}$ が求まることになる.
じゃあプログラムして動かしてみよう.まずは、上の $f$ や $J$ が必要だ.
ちなみに、ベクトル $\boldsymbol{x}_{{old}}$ については vx とでも書いておくか.
|
|
これで準備ができたので、Newton 法をいきなり動かしてみよう.
Julia の文法で細かいところは、授業時に説明するから心配しなくていい.
|
|
1: [1.0, 0.0], [2.0, 0.0] 2: [1.0, -1.0], [0.0, 2.0] 3: [0.5, -0.75], [0.0, 0.375] 4: [0.35, -0.675], [0.0, 0.03375] 5: [0.333537, -0.666768], [0.0, 0.000406566] 6: [0.333333, -0.666667], [0.0, 6.19482e-8]
ふむ、6回程度で、なかなか良い近似値が得られた様子だ.
まず、sq の最後の列が「良い近似解」なので、こいつを抜き出しておく.
|
|
2-element Array{Float64,1}: 0.333333 -0.666667
これを使って、(近似)誤差の履歴を求めておこう.
|
|
こいつをもとに、$x, y$ それぞれの誤差がどう変化したかグラフにする.
|
|
どうやら、「誤差が減り始めると、あとは順調に減っていく」様子が見える.まあ、Taylor 展開に基づく近似解法だから、Taylor 展開が1次オーダーで十分に近似できるところまでくればあとは OK だ、ということだな.
近似解の挙動を図でみてみよう
まず、近似解の変遷そのものを全部列挙してみよう.
|
|
2×6 Array{Float64,2}: 1.0 1.0 0.5 0.35 0.333537 0.333333 0.0 -1.0 -0.75 -0.675 -0.666768 -0.666667
ふむ、ではこの近似解の動いた範囲で、関数のノルムがどうなっているか、図にするためにおおよそ計算してみよう.
|
|
25×29 Array{Float64,2}: 2.34401 2.09667 1.85914 1.63157 … 1.48325 1.54932 1.61308 1.67523 2.31653 2.06803 1.82921 1.60017 1.51906 1.58682 1.6522 1.7159 …略… 3.21224 2.97129 2.74296 2.52799 2.27371 2.37483 2.47741 2.58157 3.32602 3.08662 2.86 2.64689 2.34265 2.44369 2.54661 2.65149
両方あわせてプロットしてみよう
|
|
近似解が $f(x,y)$ の値を小さくする方向へ動いていることがわかるだろう.
Report No.01 自力でやってみよう
上の例を参考にすれば、ほとんど即座にできるはずだ.そこで、できれば、見ないでやってみることにチャレンジしてみよう.
- $\sin(x) = 3 / 4$ の近似解を求めてみよう.
- 連立方程式
\[ \left\{ \begin{array}{ccc} x + 2y & = & -1, \cr \sin(x + y) & = & 0.7 \end{array} \right. \]
の近似解を求めてみよう.