Newton 法(非線形方程式の求解)

スカラーの非線形問題

Newton法の理屈は授業で教えよう.理屈は大変に簡単だ.

そして今回は、その理屈を「実際に動かせるか」というところを頑張ってみよう.

例えば, $x^2 = 2$ を解いて解 $\sqrt{2}$ を求めてみよう.

典型的な手法は、まずは問題を 左辺 = 0 の形に直すところから始まる.

1
f(x) = x^2 - 2   # この関数 f に対し、 f(x) = 0 を満たす値こそ、求めたい値だ.

さらに、Newton 法にはこの関数の微分形も必要だ. 今回は幸い手で求まるので、与えてしまおう.

1
df(x) = 2x

さて実際に Newton 法を実行してみよう.Newton 法は近似解がだんだん良くなっていく(と期待する)反復解法なので、初期値を設定してから、ループを回すという形になる.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# 初期値.なるべく良い近似値がいいのだが, まあ適当に与えておこう.
x = 1.0

# あとでグラフにしてみたいので、途中経過を全て保存しておくベクトルを用意しよう.
sq = [ x ]  

# 反復改善する.まあ、100回ぐらいで停まるようにしておこう.
for i in 1:100
    println(i, ": ", x, ", ", f(x)) # 状況を見たいので、逐一 出力しよう.
    if abs(f(x)) < 1.0e-06          # 6桁ぐらいあってればいいや、という判断ならばこんな感じだ.
        break
    end
    x = x - f(x)/df(x)  # 肝心の Newton 反復改善部分がこれだ.
    push!(sq, x)        # 新しい近似値をベクトルの要素として増やしておこう
end

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
2, x

(1.4142135623730951, 1.4142135623746899)

大変近いな.実際の差も見ておこう.

1
2 - x

-1.5947243525715749e-12

関数だけでなく、近似値の値そのものも 12桁ぐらいあっていることがわかる. この計算はうまくいったと言ってよいだろう.

誤差 abs$(\sqrt{2} - x)$ がどのように変化していくか、グラフにしてみよう. まずはいつもの準備をして…

1
2
using Plots
gr()

上で作った sq ベクトルを使って次のようにすればよいだろう.

1
plot(abs(sq - 2), marker = :auto)

svg

誤差が急速に小さくなっていることが分かる.

さらに、対数グラフにしてみようか.

1
plot(log10(abs(sq - 2)), marker = :auto)

svg

対数グラフで見ると、誤差が急速に小さくなっていることがよりはっきりと見て取れる.たいしたものだ.

欲張って、20桁の精度が出るかやってみよう

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# 初期値.なるべく良い近似値がいいのだが, まあ適当に与えておこう.
x = 1.0

# あとでグラフにしてみたいので、途中経過を全て保存しておくベクトルを用意しよう.
sq = [ x ]  

# 反復改善する.まあ、100回ぐらいで停まるようにしておこう.
for i in 1:100
    println(i, ": ", x, ", ", f(x)) # 状況を見たいので、逐一 出力しよう.
    if abs(f(x)) < 1.0e-20          # 20桁あってて欲しい!
        break
    end
    x = x - f(x)/df(x)  # 肝心の Newton 反復改善部分がこれだ.
    push!(sq, x)        # 新しい近似値をベクトルの要素として増やしておこう
end

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 とでも書いておくか.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
function f(vx)
    x,y = vx  # ベクトル vx の第一要素と第二要素をそれぞれ x,y として受ける
    return [  # [] の中に要素をスペース区切りで並べて書けばベクトルや行列になる
        x+2y+1
        x^2+2*y^2-1
    ]
end

function J(vx)
    x,y = vx
    return [
        1 2
        2x 4y
    ]
end

これで準備ができたので、Newton 法をいきなり動かしてみよう.

Julia の文法で細かいところは、授業時に説明するから心配しなくていい.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# 初期値.なるべく良い近似値がいいのだが, まあ適当に与えておこう.
vx = [ 1.0, 0.0 ]

# あとでグラフにしてみたいので、途中経過を全て保存しておく入れ物を用意しよう.
sq = copy(vx)

# 反復改善する.まあ、100回ぐらいで停まるようにしておこう.
for i in 1:100
    println(i, ": ", vx,", ", f(vx)) # 状況を見たいので、逐一 出力しよう.
    if norm(f(vx)) < 1.0e-06   # f(vx) はベクトルなので、大きさは norm で.
        break
    end
    vx = vx - J(vx)\f(vx)  # 肝心の Newton 反復改善部分がこれだ.
    sq = hcat(sq, vx)        # 新しい近似値を「横に加える」形で記録していこう.
end

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 の最後の列が「良い近似解」なので、こいつを抜き出しておく.

1
sol = sq[:,end]

2-element Array{Float64,1}: 0.333333 -0.666667

これを使って、(近似)誤差の履歴を求めておこう.

1
2
3
4
error = similar(sq)
for j in 1:6
    error[:,j] = sq[:,j] - sol
end

こいつをもとに、$x, y$ それぞれの誤差がどう変化したかグラフにする.

1
plot(error', marker = :circle)

svg

どうやら、「誤差が減り始めると、あとは順調に減っていく」様子が見える.まあ、Taylor 展開に基づく近似解法だから、Taylor 展開が1次オーダーで十分に近似できるところまでくればあとは OK だ、ということだな.

近似解の挙動を図でみてみよう

まず、近似解の変遷そのものを全部列挙してみよう.

1
sq

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

ふむ、ではこの近似解の動いた範囲で、関数のノルムがどうなっているか、図にするためにおおよそ計算してみよう.

1
2
3
4
X = 0:0.05:1.2
Y = -1.2:0.05:0.2

fz = [ norm(f([x,y])) for x in X, y in Y ]

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

両方あわせてプロットしてみよう

1
2
plot(X, Y, fz', aspect_ratio = 1)  # fz の等高線図を、x,y のアスペクト比1 で描く
plot!(sq[1,:], sq[2,:], marker = :circle)  # 近似解の挙動はこちら

svg

近似解が $f(x,y)$ の値を小さくする方向へ動いていることがわかるだろう.

Report No.01 自力でやってみよう

上の例を参考にすれば、ほとんど即座にできるはずだ.そこで、できれば、見ないでやってみることにチャレンジしてみよう.

  1. $\sin(x) = 3 / 4$ の近似解を求めてみよう.
  2. 連立方程式

\[ \left\{ \begin{array}{ccc} x + 2y & = & -1, \cr \sin(x + y) & = & 0.7 \end{array} \right. \]

の近似解を求めてみよう.