スカラー非線形問題を解くには…¶

まずは Newton 法を学ぶ¶

Newton法の理屈は授業で学んだ通り．そこで、今回対象とするスキームを1ステップだけ時間発展させる式を解くという問題を Newton 法で解いてみよう．

例えば√2 を求めるには、$x^2 = 2$ を解けばいいわけなので、これを解いてみよう．¶

In [2]:
# 対象の、右辺がゼロとなる問題を記述する

f(x) = x^2 - 2
Out[2]:
f (generic function with 1 method)
In [3]:
# 式の微分形が必要だ．ここでは手計算で求めたものを入れてしまおう．

diff_f(x) = 2x
Out[3]:
diff_f (generic function with 1 method)

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

In [5]:
# 初期値．なるべく良い近似値がいいのだが．
x = 1.0

# 反復改善する．まあ、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)/diff_f(x)  # 肝心の Newton 反復改善部分．
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 を比べてみよう．

In [7]:
√2, x
Out[7]:
(1.4142135623730951, 1.4142135623746899)

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

In [8]:
√2 - x
Out[8]:
-1.5947243525715749e-12

近似値の値そのものも、12桁ぐらいあっていることがわかる．

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

In [9]:
# 初期値．なるべく良い近似値がいいのだが．
x = 1.0

# 反復改善する．まあ、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)/diff_f(x)  # 肝心の Newton 反復改善部分．
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桁以上の精度を求めようとするのはそのままでは無理ゲーだ、ということを知っておくことだ．何桁の精度が欲しいのか、常に意識しておくと健全だ．

自力でやってみよう¶

$sin(x) = 3/4$ の近似解を求めてみよう¶

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

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

問題は、 $\left\{ \begin{eqnarray} x + 2y & = & 0, \\ x^2 + 2y^2 & = & 0 \end{eqnarray} \right.$ の近似解を求めることにしておこう．

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

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