Cambridge Hall of Fame featuring Sir Isaac Newton (centre) flanked by the likes of Sir Alfred Tennyson and Francis Bacon in the Anti-chapel at Trinity College, Cambridge (Mar., 2008).
Photo by K. Mitch Hodge on Unsplash

ゼロ点を求める¹

コンピュータで計算する！ というとたいそう難しく聞こえる． しかしまあ，おおざっぱには多くの問題は以下のふたつの問題を組み合わせたものだと考えられるだろうから，なんとかなりそうな気もするね？

制約を満たす解を求める問題

例. $\sin(x) = 0.85$ を満たす $x$ を求める問題を考えてみよう²．

  なにか制約を満たすものを「解」と定義して，それもしくはその近似を求める問題だ．
このタイプの問題には，求めた近似解をその制約式に代入するなどで妥当性を確認できる，というありがたい特徴がある． つまり，近似計算結果が(それなりに)正しいかどうか，チェックできるのだ．

なにかを最大・最小にする解を探す問題(最適化問題)

例. $x \in [-0.5,0.5]$ で $y = - x^2 - 5\exp(- 10^{5} \times x^2)$ (下記の図)の最小値を求めたいという問題． $y$ のグラフを下図のようにこまか～く描けばすぐわかるが，それには大きなコストがかかるし，運が悪ければそうしたとしてもうまく見つからない.

  最大・最小値を求めたいというタイプの問題は，それ自体が本質的に困難だ．探し方が悪いと見つからない． しかも，本質的には，得た近似解が本当に最小値に近いのかを調べる方法が無い.

幸いなことに？ この授業の対象は基本的に前者に近い問題が多い． 得た近似解がどれくらい妥当か，確認できる世界にいるときは安心しよう．

そして今回はその中でも単純な，入力も出力も単なる実数である関数を考えて，その関数値がゼロになる入力値(ゼロ点)を求める，という問題, すなわち，一般的に \begin{equation} f(x) = 0 \,\, \mbox{ for } {}^{\exists} x \in \mathbb{R} \end{equation} な状況を想定してこの $x$ を求める問題を考えよう．

Newton 法を学ぶ

今回は，幅広いシーンで使えて，理解も容易で，かつ使いやすい Newton 法というものを学ぼう．

これは

「今持っている近似値より少しはマシな近似値を求める」=「改善する」

試みを繰り返す方法(反復法と呼ぶよ)³の一種だ． だから，「どうやって近似解を改善するのか」を知れば全体がわかることになる.

  上で「試み」と書いたように，Newton法では反復によって近似解が改善される保証はない．あくまで「期待する」だけだ．

さて，その肝心の Newton 法の反復の 1ステップ分の改善プロセスは，下記の図のようになっている．

1. 今もっている近似値 $x_{old}$ を用いて，関数値 $f(x_{old})$ とその微分値 $f'(x_{old})$ を「数値として」計算する．
   
   
2. $(x_{old}, f(x_{old}))$ を通って傾きが $f'(x_{old})$ な直線の式 \begin{equation} y - f(x_{old}) = f'(x_{old}) ( x - x_{old} ) \end{equation} を考える．
   
   
3. この直線が $x$軸と交わるところを「改善された次の近似値 $x_{new}$」だと定める⁴．
   つまり，この直線の式に $y=0$ を代入して求まる $x$ がそれだということになるので， \begin{equation} x_{new} = x_{old} - \frac{f(x_{old})}{f'(x_{old})} \end{equation} として新たな近似値を手に入れる．

  上の説明からわかるように，Newton 法は関数の「微分値」が計算できるときしか使えない． 微分計算ができない関数を相手にする場合は他の方法を考えないといけないな．

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

さて，では実際に Julia のプログラムを動かしてやってみよう． ここからはこの Newton法を「実際に動かせるか」頑張ってみるのだ．

  余裕があるならば，下記の内容を Julia を動かせる環境で試しながら聞こう．

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

1f(x) = x^2 - 2   # この関数 f のゼロ点が求めたい値だ．

さらに、Newton 法にはこの関数の微分形も必要だ． 今回は幸い容易にわかるので、具体的に与えてしまおう．

1df(x) = 2x

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

 1# 初期値．なるべく良い近似値がいいのだが, まあ適当に与えておこう．
 2x = 1.0
 3
 4# グラフにしたいので、途中経過を配列に保存していく
 5sq = [ x ]  
 6
 7# 反復改善する．うまくいかなくても 100回で停まる
 8for i in 1:100
 9    println("$i : $x, $(f(x))" ) # 状況を出力
10    if abs(f(x)) < 1.0e-06     # f(x) が 0 に十分近ければ OK！
11        break
12    end
13    x = x - f(x)/df(x)  # Newton 法による近似値更新！
14    push!(sq, x)        # 新しい近似値を保存
15end

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)

  Julia では √ は \sqrt と打ってから TABキー を押すと入力できるぞ． √2 とする代わりに sqrt(2) としても良い．このあたりは好みの問題だ．

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

1√2 - x

-1.5947243525715749e-12

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

誤差 abs$(\sqrt{2} - x)$ がどのように変化していくか、グラフにしてみよう． まずは下記のように Plots package の利用宣言をして,

1using Plots

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

1error = abs.( sq .- √2 )  # 近似解から真値を引いて，誤差配列を作る．
2plot(error, marker = :circle)

上のプログラム中の，abs. や .- などの，オペレータに .(ドット) をつけた記述は，通常のスカラー計算を配列の各要素にそれぞれ作用させるという機能で，broadcast などと呼ばれる． 気になる人はマニュアルの broadcasting という箇所を読んでみよう．

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

1plot(error, marker = :circle, yaxis = :log10)   #  y軸が対数に．

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

今の方法で，さらに高精度な近似解が得られないか，欲張ってみよう． それには，下記のように収束の判定部分を厳しくしてみる，というのが一番手っ取り早いのだが，さてどうなるかな？

 1# 初期値．なるべく良い近似値がいいのだが, まあ適当に与えておこう．
 2x = 1.0
 3
 4# グラフにしたいので、途中経過を配列に保存していく
 5sq = [ x ]  
 6
 7# 反復改善する．うまくいかなくても 100回で停まる
 8for i in 1:100
 9    println("$i : $x, $(f(x))" ) # 状況を出力
10    if abs(f(x)) < 1.0e-50     # 50桁あうか？
11        break
12    end
13    x = x - f(x)/df(x)  # Newton 法による近似値更新！
14    push!(sq, x)        # 新しい近似値を保存
15end

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
…略…
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

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

  今回のプログラムは最悪でも 100回計算したら終わるように作ってあったので，実質的に無限ループに陥ってもプログラム自体は終了した．安全装置がうまく機能したわけだ． このように，プログラムには「安全装置をつけておく」癖をつけよう．

  さて，今回計算がうまく収束しなかったのは，特に指定しなければ Julia では浮動小数点に有効桁数が約 16桁のシステム(float64, 倍精度とも言うね)が使われるからだ． この制限は Julia に限らず多くの環境で共通だ．これをなんとかするのには「多倍長変数」という，高精度な変数に基づく計算をする必要がある． 詳しくは Julia の BigFloat のマニュアルなどを読むと良いが，まあ，簡単に体験してみよう．次のプログラムを動かしてみよう．

 1# 多倍長計算を使ってみる．big で簡単に変換できるが，できれば BigFloat 命令などを使ったほうが正確だ．
 2
 3# 初期値．なるべく良い近似値がいいのだが, まあ適当に与えておこう．
 4x = big(1.0) 
 5# 1.0 を BigFloat 型に変換して x に格納する．
 6# このように精度を指定しない場合，10進法でおおよそ77桁の精度になる．
 7
 8# グラフにしたいので、途中経過を配列に保存していく
 9sq = [ x ]  
10
11# 反復改善する．うまくいかなくても 100回で停まる
12for i in 1:100
13    println("$i : $x, $(f(x))" ) # 状況を出力
14    if abs(f(x)) < 1.0e-50     # 50桁あうか？
15        break
16    end
17    x = x - f(x)/df(x)  # Newton 法による近似値更新！
18    push!(sq, x)        # 新しい近似値を保存
19end

1 : 1.0, -1.0
2 : 1.5, 0.25
3 : 1.416666666666666666666666666666666666666666666666666666666666666666666666666661, 0.006944444444444444444444444444444444444444444444444444444444444444444444444417576
4 : 1.414215686274509803921568627450980392156862745098039215686274509803921568627451, 6.007304882737408688965782391387927720107650903498654363706266820453671663376355e-06
5 : 1.414213562374689910626295578890134910116559622115744044584905019200054371835385, 4.510950444942772099280764360710487611553552636972638583381341682774977026972629e-12
6 : 1.414213562373095048801689623502530243614981925776197428498289498623195824228933, 2.543584239585437205842792747002928940919398699981858849216171323264197378926519e-24
7 : 1.41421356237309504880168872420969807856967187537723400156101313311326525563035, 8.087275979834283525943226474966356195438464890457763565398783053640267780573221e-49
8 : 1.414213562373095048801688724209698078569671875376948073176679737990732478462102, -1.727233711018888925077270372560079914223200072887256277004740694033718360632485e-77

おお，良い近似解 $x$ が計算出来ている．なんと，関数 $f(x)$ の誤差が $10^{-77}$ まで小さくなるとは…

ちなみに，この近似値 $x$ を使って下記のように $x^2$ を計算してみるとその精度が実感できる(数えてみると，"9"が 76回続いている…)．

1x^2

1.999999999999999999999999999999999999999999999999999999999999999999999999999983

このように，多倍長計算を使えばかなり高精度な近似解を求められることがわかるだろう．

ただし，通常の倍精度計算の多くが CPU のハードウェアでできるために速いことに比べると，多倍長変数を用いた計算はソフトウェアでいろいろ工夫するためにどうしても計算時間はかかる．メモリーも余計に必要だ． この「計算コストの高さ」は覚悟しておこう．

だからどうしてもという時以外は16桁以上の精度を求めようとしない方が良いと認識しておこう⁵．

発展: 多次元問題の場合は？

気になる人もいるだろうから，$n > 1$ で $f$ が $\mathbb{R}^n \rightarrow \mathbb{R}^n$ な写像のときの Newton 法も本質的には同じ考えだよ，ということをさらっと書いておこう．

ただし図示による説明は困難なので Taylor 展開で説明しよう．

前提は一緒で，最初に近似値 $\boldsymbol{x}_{old} \in \mathbb{R}^n$ を持っているとする． そして，次の(改善)近似値を $\boldsymbol{x}_{new} \in \mathbb{R}^n$ と書こう．

このとき，関数 $f(\boldsymbol{x}_{new})$ を，$\boldsymbol{x}_{old}$ を中心として 1次項まで Taylor 展開すると， \begin{equation} f(\boldsymbol{x}_{new}) \cong f(\boldsymbol{x}_{old}) + J_f(\boldsymbol{x}_{old})\cdot( \boldsymbol{x}_{new} - \boldsymbol{x}_{old} ) \end{equation} となることは大学入ってすぐの 1年生の微積分学で学んだよね？
ちなみに $J_f(\boldsymbol{x}_{old})$ は $(i,j)$成分が $\partial f_i/\partial x_j$ であるヤコビ行列だね．

そして，$\boldsymbol{x}_{new}$ は「改善された」ゼロ点近似値であることを期待している． これは $f(\boldsymbol{x}_{new}) \cong \boldsymbol{0}$ だと期待していることであるので，上の式は \begin{equation} \boldsymbol{0} \cong f(\boldsymbol{x}_{old}) + J_f(\boldsymbol{x}_{old})\cdot( \boldsymbol{x}_{new} - \boldsymbol{x}_{old} ) \end{equation} と期待することになる． だからあとは上の式を単純に変形して，連立一次方程式 \begin{equation} J_f(\boldsymbol{x}_{old})\cdot ( \boldsymbol{x}_{new} - \boldsymbol{x}_{old} ) \cong - f(\boldsymbol{x}_{old}) \end{equation} を解いて $\boldsymbol{x}_{new} - \boldsymbol{x}_{old}$ を数値的に求めて，それに $\boldsymbol{x}_{old}$ を足せば $\boldsymbol{x}_{new}$ が求まるということになる．
数式で書くなら \begin{equation} \boldsymbol{x}_{new} \cong \boldsymbol{x}_{old} - \left( J_f(\boldsymbol{x}_{old}) \right)^{-1} f(\boldsymbol{x}_{old}) \end{equation} と変形されるわけだね．

ただし，この数式をみてプログラムを組むときに， 行列 $J_f(\boldsymbol{x}_{old})$ の 逆行列 $\left( J_f(\boldsymbol{x}_{old}) \right)^{-1}$ を計算してはいけない！

逆行列を直接求める計算量はとんでもなく大きいうえに誤差もたっぷりはいってくる⁶．良いことはないのだ． 思い出そう．「連立一次方程式を解く」のに逆行列は必要ないことを．

あくまで「連立一次方程式を解く」計算でプログラムを組もう．

実際にプログラムを組んでやってみる

実際に \begin{equation} \left\{ \begin{array}{rcl} x + 2y & = & -1, \cr x^2 + 2y^2 & = & 1 \end{array} \right. \end{equation} の近似解を Newton 法で求めてみよう． 初期値を決めないといけないが，まあ適当に $(1.0, -1.0)$ あたりで．

Julia で連立一次方程式を解くにはどうすれば良い？

プログラムを書く前に，Julia で連立一次方程式を解く文法だけ簡単に説明しておこう．

Julia では連立一次方程式 $A\boldsymbol{x} = \boldsymbol{b}$ を解くには

1x = A \ b

と書くだけで良い⁷,⁸． これでこの連立一次方程式が解かれて，$\boldsymbol{x}$ がその解になる．

さて，話をすすめよう．まずはゼロになってほしいベクトル値関数とその微分であるヤコビ行列の定義だ．以下のような感じだ．

 1# ゼロになってほしい関数
 2function f(v)    # 関数を定義するたぶん最も正式な方法
 3    x, y = v     # ベクトル v を成分 x,y にわけている
 4    return [     # 出力は 2x1 の縦ベクトルで．
 5      x+2y+1     # 式中にスペースを入れるな！ 行列の別要素とされてしまう
 6      x^2+2y^2-1   
 7    ]
 8end
 9
10# f の微分であるヤコビ行列
11function J(v)
12    x, y = v
13    return [    # 出力は 2x2 の行列で．
14        1 2
15        2x 4y
16    ]
17end

関数の定義をしたら，まずは次のようにして動作確認をしておこう．

1f( [2.0, 3] ) # ← 関数 f はベクトルを受けるつもりで定義したので，ベクトルを渡す．

2-element Vector{Float64}:
  9.0  # この場合は x+2y+1 = 2+2*3+1 = 9 のはずなので合っている．
 21.0  # この場合は x^2+2y^2-1 = 2^2+2*3^2-1 = 21 のはずなので合っている．

1J( [2.0, 3] ) # ← 行列値の関数 J もベクトルを受ける．

2×2 Matrix{Float64}:
 1.0   2.0  # これはいつも 1 と 2.
 4.0  12.0  # ここは 2x = 2*2 = 4 と 4y = 4*3 = 12 のはずで，合っている．

というわけで，自分で定義した関数 $f$, $J$ が合っていそうなことを簡単に確かめることが出来た．

次に，ベクトルのノルム計算を楽にやりたいので，線形計算ライブラリ LinearAlgebra を読み込む． これは標準インストール済みのライブラリなので，インストールは不要で呼び出すだけで良い．

1using LinearAlgebra

あとは先の例と同様に Newton 法を適用するだけだ．細かい部分は異なるが，全体は同じなことがわかるだろう．

 1# 初期値．なるべく良い近似値がいいのだが, まあ適当に与えておこう．
 2x = [1.0, -1.0]
 3
 4# グラフにしたいので、途中経過を(2次元)配列に保存していく
 5sq = x 
 6
 7# 反復改善する．うまくいかなくても 100回で停まる
 8for i in 1:100
 9    println("$i : $x, $(norm(f(x)))" ) # 状況を出力
10    if norm(f(x)) < 1.0e-06     # f(x)のノルムが 0に近ければ OK！
11        break
12    end
13    x = x - J(x) \ f(x)  # Newton 法による近似値更新！
14    sq = hcat(sq, x)      # 新しい近似値を横につなげる形で保存
15end

1 : [1.0, -1.0], 2.0
2 : [0.5, -0.75], 0.375
3 : [0.35, -0.675], 0.03375000000000017
4 : [0.33353658536585357, -0.6667682926829268], 0.0004065660321237452
5 : [0.33333336430743143, -0.6666666821537157], 6.19481976826819e-8

結果を見ると，近似値の更新が 4回行われて，5回めループの最初で「終了」したことがわかる．

試しに，2次元平面上で近似値がどう動いたかも見てみよう． 次のようにすると見ることが出来る．

1using Plots
2
3plot(sq[1,:], sq[2,:], 
4     st = :path, marker = :circle, markersize = [4,5,6,7,8])

等高線グラフで納得してみる

上の近似解とその変化を，等高線グラフを描くことで納得してみよう．次のような感じだ．

1Rx = 0:0.05:2    # $x$ の範囲
2Ry = -1:0.025:0  # $y$ の範囲
3
4F = [ norm(f((x,y))) for y in Ry, x in Rx ]
5# 行列の添字の動き方は xy-平面の軸と 90度違うので，先に y を動かしておく．
6
7contour(Rx, Ry, F, fill = true )

レポート

以下の課題について作業と考察を行い，
2024-numerical-analysis-report-03
という題名をつけて e-mail にて教官宛にレポートとして提出しよう. なお，e-mail ではなく，TeX などで作成した書面の形で作成したレポートを教員に手渡しで提出してもよい.

  注意
　　近年はセキュリティ上の懸念から，実行形式のプログラムなどをメールに添付するとそのメールそのものの受信を受信側サーバが拒絶したりする． そういうことを避けるため，レポートをメールで提出するときは添付ファイルにそういった懸念のあるファイルが無いようにしよう．

まあ要するに，レポートは pdf ファイルにして，それをメールに添付して送るのが良い ということだと思っておこう．

  レポートやメールには，自分の 所属・学年・学籍番号・氏名 を必ず書こう (経験上，数年に一度，差出人不明のレポートが届くので気をつけよう).

課題

1. 自分で Newton 法を用いて $\sin(x) = 3 / 4$ の近似解を求めてみよう
   上の例を参考にすれば、ほとんど即座にできるはずだ． 可能ならば，上の例を見ないでやってみよう．
   
   
2. 次に、応用問題として \begin{equation} \left\{\begin{array}{rcl} x + 2y + z & = & -1, \cr x^2 + 2y^2 + z^2 & = & 10, \cr \sin(x+y+z) & = & 0.7 \end{array}\right. \end{equation} の近似解を Newton 法で求めてみよう． 考え方は上の例と同じだ．
   $[x,y,z]$ の初期値はそうだなあ，$[1.0, 0.0, 0.0]$ あたりでいいだろう．