まず手始めに、スカラーの非線形問題が解けるか頑張ってみよう

Newton 法を学ぶ

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

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

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

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

1

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

f (generic function with 1 method)

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

1

df(x) = 2x

df (generic function with 1 method)

さて実際に 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          # f(x) が 0 と比べて 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

using Plots

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

1
2

error = abs.( sq .- √2 )
plot(error, marker = :auto)

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

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

1

plot(log10.(error), marker = :auto)

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

欲張って、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          # f(x) が 0 と比べて 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 8: 1.4142135623730951, 4.440892098500626e-16 9: 1.414213562373095, -4.440892098500626e-16 10: 1.4142135623730951, 4.440892098500626e-16 11: 1.414213562373095, -4.440892098500626e-16 12: 1.4142135623730951, 4.440892098500626e-16 13: 1.414213562373095, -4.440892098500626e-16 14: 1.4142135623730951, 4.440892098500626e-16 15: 1.414213562373095, -4.440892098500626e-16 16: 1.4142135623730951, 4.440892098500626e-16 17: 1.414213562373095, -4.440892098500626e-16 18: 1.4142135623730951, 4.440892098500626e-16 19: 1.414213562373095, -4.440892098500626e-16 20: 1.4142135623730951, 4.440892098500626e-16 21: 1.414213562373095, -4.440892098500626e-16 22: 1.4142135623730951, 4.440892098500626e-16 23: 1.414213562373095, -4.440892098500626e-16 24: 1.4142135623730951, 4.440892098500626e-16 25: 1.414213562373095, -4.440892098500626e-16 26: 1.4142135623730951, 4.440892098500626e-16 27: 1.414213562373095, -4.440892098500626e-16 28: 1.4142135623730951, 4.440892098500626e-16 29: 1.414213562373095, -4.440892098500626e-16 30: 1.4142135623730951, 4.440892098500626e-16 31: 1.414213562373095, -4.440892098500626e-16 32: 1.4142135623730951, 4.440892098500626e-16 33: 1.414213562373095, -4.440892098500626e-16 34: 1.4142135623730951, 4.440892098500626e-16 35: 1.414213562373095, -4.440892098500626e-16 36: 1.4142135623730951, 4.440892098500626e-16 37: 1.414213562373095, -4.440892098500626e-16 38: 1.4142135623730951, 4.440892098500626e-16 39: 1.414213562373095, -4.440892098500626e-16 40: 1.4142135623730951, 4.440892098500626e-16 41: 1.414213562373095, -4.440892098500626e-16 42: 1.4142135623730951, 4.440892098500626e-16 43: 1.414213562373095, -4.440892098500626e-16 44: 1.4142135623730951, 4.440892098500626e-16 45: 1.414213562373095, -4.440892098500626e-16 46: 1.4142135623730951, 4.440892098500626e-16 47: 1.414213562373095, -4.440892098500626e-16 48: 1.4142135623730951, 4.440892098500626e-16 49: 1.414213562373095, -4.440892098500626e-16 50: 1.4142135623730951, 4.440892098500626e-16 51: 1.414213562373095, -4.440892098500626e-16 52: 1.4142135623730951, 4.440892098500626e-16 53: 1.414213562373095, -4.440892098500626e-16 54: 1.4142135623730951, 4.440892098500626e-16 55: 1.414213562373095, -4.440892098500626e-16 56: 1.4142135623730951, 4.440892098500626e-16 57: 1.414213562373095, -4.440892098500626e-16 58: 1.4142135623730951, 4.440892098500626e-16 59: 1.414213562373095, -4.440892098500626e-16 60: 1.4142135623730951, 4.440892098500626e-16 61: 1.414213562373095, -4.440892098500626e-16 62: 1.4142135623730951, 4.440892098500626e-16 63: 1.414213562373095, -4.440892098500626e-16 64: 1.4142135623730951, 4.440892098500626e-16 65: 1.414213562373095, -4.440892098500626e-16 66: 1.4142135623730951, 4.440892098500626e-16 67: 1.414213562373095, -4.440892098500626e-16 68: 1.4142135623730951, 4.440892098500626e-16 69: 1.414213562373095, -4.440892098500626e-16 70: 1.4142135623730951, 4.440892098500626e-16 71: 1.414213562373095, -4.440892098500626e-16 72: 1.4142135623730951, 4.440892098500626e-16 73: 1.414213562373095, -4.440892098500626e-16 74: 1.4142135623730951, 4.440892098500626e-16 75: 1.414213562373095, -4.440892098500626e-16 76: 1.4142135623730951, 4.440892098500626e-16 77: 1.414213562373095, -4.440892098500626e-16 78: 1.4142135623730951, 4.440892098500626e-16 79: 1.414213562373095, -4.440892098500626e-16 80: 1.4142135623730951, 4.440892098500626e-16 81: 1.414213562373095, -4.440892098500626e-16 82: 1.4142135623730951, 4.440892098500626e-16 83: 1.414213562373095, -4.440892098500626e-16 84: 1.4142135623730951, 4.440892098500626e-16 85: 1.414213562373095, -4.440892098500626e-16 86: 1.4142135623730951, 4.440892098500626e-16 87: 1.414213562373095, -4.440892098500626e-16 88: 1.4142135623730951, 4.440892098500626e-16 89: 1.414213562373095, -4.440892098500626e-16 90: 1.4142135623730951, 4.440892098500626e-16 91: 1.414213562373095, -4.440892098500626e-16 92: 1.4142135623730951, 4.440892098500626e-16 93: 1.414213562373095, -4.440892098500626e-16 94: 1.4142135623730951, 4.440892098500626e-16 95: 1.414213562373095, -4.440892098500626e-16 96: 1.4142135623730951, 4.440892098500626e-16 97: 1.414213562373095, -4.440892098500626e-16 98: 1.4142135623730951, 4.440892098500626e-16 99: 1.414213562373095, -4.440892098500626e-16 100: 1.4142135623730951, 4.440892098500626e-16

どうやら実質的に無限ループに陥ってしまったようだ．

これは、(特に指定しなければ)有効桁数が約 16桁のシステム(float64)が使われるからだ．この制限は Julia に限らず多くの環境で共通だ．これをなんとかするのには少し知識が必要なので、またいずれにしておこう．

ここで重要なのは、多くの計算機・プログラム環境ではそもそも約16桁以上の精度を求めようとするのはそのままでは無理ゲーだ、ということを知っておくことだ．何桁の精度が欲しいのか、常に意識しておくと健全だ．

Report No.01: Newton 法を自力でやってみよう

自分で Newton 法を用いて $\sin(x) = 3 / 4$ の近似解を求めてみよう

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

次に、応用問題として、非線形連立方程式の近似解も求めてみよう．

問題は

\[ \left\{ \begin{eqnarray*} x + 2y & = & -1, \cr x^2 + 2y^2 & = & 1 \end{eqnarray*} \right. \]

の近似解を求めるものとしよう．

一見どうしたら良いかわからないかもしれないが、実は、Newton 法がほぼそのまま使える． 数式を手で書いて、「Newton 法の成り立ち」をよく考えてみればわかるだろう．

よくわからない人は、わかっていそうな人や、教員に遠慮なく尋ねてみよう．