コンテンツにスキップ

非線形方程式の反復解法:二分法・不動点反復・Newton 法・割線法の収束を証明する

対象読者: 学部 1〜2 年の微積分(平均値の定理・Taylor の定理)と、距離空間の縮小写像の原理を学んだ人。本テーマ第 1〜2 章の誤差の考え方は既知としますが、反復法の知識は仮定しません。

前提となる記事:平均値の定理とテイラーの定理:微分から関数の姿を復元する、距離空間:完備化・Baire のカテゴリー定理・縮小写像の原理・コンパクト性

生 Markdown

電卓で 2\sqrt2 を求めるとき、中では「だいたいの答えを出し、それを少し直して、また直す」という繰り返しが行われています。式で解けない方程式でも、同じ発想で答えに近づけます。この記事では、そのやり方がいくつかあること、どれが速いか、そして本当に答えに近づいていると言える理由を、順に確かめます。読み終えると、状況に応じてどの方法を選べばよいかと、その根拠を説明できるようになります。

  • 方程式 f(x)=0f(x) = 0 の解を求める反復法は、誤差 ek=xk−x∗e_k = x_k - x^* がどんな速さで 00 に近づくかで比べます。速さは「次数」(誤差が次の回に何乗になるか)と「漸近誤差定数」(そのときの比例係数)の 2 つの数で表します。
  • 二分法は連続性だけで動き、kk 回で誤差が (b−a)/2k+1(b-a)/2^{k+1} 以下になります。遅いが必ず収束します。
  • 不動点反復 xk+1=g(xk)x_{k+1} = g(x_k) は、∣g′(x∗)∣<1|g'(x^*)| \lt 1 なら 1 次収束し、gg の p−1p-1 階までの導関数が x∗x^* で消えれば pp 次収束します。証明は平均値の定理と Taylor の定理です。
  • Newton 法は単根の近くで 2 次収束し、ek+1/ek2→f′′(x∗)/(2f′(x∗))e_{k+1}/e_k^2 \to f''(x^*)/(2f'(x^*)) です。重根では 1 次収束に落ち、比は (m−1)/m(m-1)/m です。凸関数では初期値を解の右に取れば単調に収束します。
  • 割線法は導関数を使わずに、次数が黄金比 (1+5)/2≈1.618(1+\sqrt5)/2 \approx 1.618 以上になります。
  • 多変数の Newton 法も 2 次収束し、Kantorovich 型の定理は初期値での計算だけから解の存在と一意性を保証します。証明は縮小写像の原理に帰着します。

x2=2x^2 = 2 の正の解、x=cos⁡xx = \cos x の解、ある固有値問題の特性方程式の根。解を式で書けない方程式は、近似解の列 x0,x1,x2,…x_0, x_1, x_2, \ldots を作って解に近づけるしかありません。ここで問うべきことは 3 つあります。

  1. 列は解に収束するか。どんな初期値から出発すればよいか。
  2. 収束するならどれだけ速いか。正しい桁数が 1 回の反復で 1 桁ずつ増えるのか、2 倍ずつ増えるのか。
  3. 計算をいつ止めてよいか。止めたときの誤差をどう見積もるか。

二分法は 1 に強く、Newton 法は 2 に強い、という対照がこの章の軸です。二分法は区間の端で符号が違えば必ず収束しますが、1 回あたり 1 ビットしか進みません。Newton 法は解の近くでは 1 回ごとに正しい桁数がおよそ 2 倍になりますが、初期値が悪いと発散したり周期軌道に入ったりします。実用の算法(Brent 法など。二分法のような確実な方法と、速い方法を切り替えて使う算法)は両者を組み合わせています。

3 の問いには、不動点反復の事後誤差評価と、Kantorovich 型の定理が答えます。前者は「直前の 2 つの反復値の差」から誤差の上界を出し、後者は初期値での値と導関数だけから、解が存在する球の半径を出します。後者は計算機援用証明(計算機の出力を数学の証明として使う方法)で使う区間 Newton 法(誤差の範囲を区間で追いながら行う Newton 法)の考え方の土台です(浮動小数と区間演算)。

以下、実数列 (xk)(x_k) が x∗x^* に収束するとき ek=xk−x∗e_k = x_k - x^* と書きます。

次の定義のねらいは、「誤差が 0 に近づく速さ」を 1 つの数で言えるようにすることです。次の回の誤差が今の誤差の何乗に比例するかを見て、その指数を次数と呼びます。

定義 2.1(収束の次数(Q 次数)と漸近誤差定数)↗

xk→x∗x_k \to x^* で、すべての kk で ek≠0e_k \ne 0 とする。p≥1p \ge 1 と 0<C<∞0 \lt C \lt \infty について

lim⁡k→∞∣ek+1∣∣ek∣p=C\lim_{k\to\infty} \frac{|e_{k+1}|}{|e_k|^p} = C

が成り立ち、さらに p=1p = 1 のときは C<1C \lt 1 であるとき、列は pp 次収束する(Q 次数が pp である)といい、CC を漸近誤差定数という。p=1,2,3p = 1, 2, 3 をそれぞれ 1 次(線形)・2 次・3 次収束という。また ∣ek+1∣/∣ek∣→0|e_{k+1}|/|e_k| \to 0 のとき超 1 次収束するという。

p=2p = 2 で CC がほどほどの大きさなら、∣ek+1∣≈C∣ek∣2|e_{k+1}| \approx C|e_k|^2 から、正しい桁数(−log⁡10∣ek∣-\log_{10}|e_k|)が 1 回ごとにほぼ 2 倍になります。1 次収束では、桁数が 1 回ごとに −log⁡10C-\log_{10} C ずつ増えるだけです。

Q 次数の定義は極限の存在を要求するので、誤差が規則正しく減らない列には使えません。そこで上界だけで速さを測る定義も用意します。2 種類あるのは、毎回きれいな比が取れる場合(Q 次数)と、上から押さえることしかできない場合(R 次数)の両方に使えるようにするためです。

定義 2.2(R 次数)↗

xk→x∗x_k \to x^* とする。p>1p \gt 1 について、定数 A>0A \gt 0 と 0<q<10 \lt q \lt 1 があってすべての kk で

∣ek∣≤A qpk|e_k| \le A\, q^{p^k}

が成り立つとき、列は R 次数 pp 以上で収束するという。p=1p = 1 では ∣ek∣≤Aqk|e_k| \le A q^k を条件とする。

qpkq^{p^k} は pkp^k が指数に乗るので、p=2p = 2 なら q,q2,q4,q8,…q, q^2, q^4, q^8, \ldots と桁数が倍々に増える上界です。以下の定理では、多くの場合に両方の意味での次数を示します。

定理 3.1(二分法の収束)✓

a<ba \lt b、f:[a,b]→Rf : [a,b] \to \mathbb{R} は連続で f(a)f(b)<0f(a) f(b) \lt 0 とする。a0=aa_0 = a、b0=bb_0 = b とし、k=0,1,2,…k = 0, 1, 2, \ldots について ck=(ak+bk)/2c_k = (a_k + b_k)/2 とおいて、f(ak)f(ck)<0f(a_k) f(c_k) \lt 0 なら (ak+1,bk+1)=(ak,ck)(a_{k+1}, b_{k+1}) = (a_k, c_k)、f(ck)f(bk)<0f(c_k) f(b_k) \lt 0 なら (ak+1,bk+1)=(ck,bk)(a_{k+1}, b_{k+1}) = (c_k, b_k) と定める(f(ck)=0f(c_k) = 0 なら ckc_k を解として止める)。止まらないとき、

  1. 各 [ak,bk][a_k, b_k] は ff の零点を含み、bk−ak=(b−a)/2kb_k - a_k = (b-a)/2^k である。
  2. (ak)(a_k) と (bk)(b_k) は共通の極限 x∗x^* に収束し、f(x∗)=0f(x^*) = 0 である。
  3. すべての kk で ∣ck−x∗∣≤(b−a)/2k+1|c_k - x^*| \le (b-a)/2^{k+1} が成り立つ。
証明(定理 3.1)

1. 帰納法で f(ak)f(bk)<0f(a_k) f(b_k) \lt 0 と bk−ak=(b−a)/2kb_k - a_k = (b-a)/2^k を示します。k=0k = 0 は仮定です。kk で成り立つとき、f(ck)≠0f(c_k) \ne 0 なので f(ck)f(c_k) は f(ak)f(a_k) と f(bk)f(b_k) のどちらか一方と異符号であり、選んだ半分の区間の端点で ff は異符号です。幅はちょうど半分になります。端点で異符号なら 00 は f(ak)f(a_k) と f(bk)f(b_k) の間にあるので、中間値の定理(定理 5.7[極限と連続性])から [ak,bk][a_k, b_k] に零点があります。

2. 作り方から ak≤ak+1a_k \le a_{k+1}、bk+1≤bkb_{k+1} \le b_k、ak<bk≤ba_k \lt b_k \le b です。(ak)(a_k) は単調増加で上に有界なので、単調有界定理(定理 2.1[級数と収束判定])から収束します。極限を x∗x^* とします。bk=ak+(b−a)/2kb_k = a_k + (b-a)/2^k なので bk→x∗b_k \to x^* でもあります。aka_k と bkb_k は [a,b][a,b] の点で、ff は x∗x^* で連続なので、数列による特徴づけ(定理 4.5[極限と連続性]。ak=x∗a_k = x^* となる項では f(ak)=f(x∗)f(a_k) = f(x^*) は自明)から f(ak)→f(x∗)f(a_k) \to f(x^*)、f(bk)→f(x∗)f(b_k) \to f(x^*) です。積の極限をとると、f(ak)f(bk)<0f(a_k) f(b_k) \lt 0 から f(x∗)2≤0f(x^*)^2 \le 0、すなわち f(x∗)=0f(x^*) = 0 です。

3. aka_k は増加して x∗x^* に、bkb_k は減少して x∗x^* に収束するので、ak≤x∗≤bka_k \le x^* \le b_k です。ckc_k はこの区間の中点なので、∣ck−x∗∣≤(bk−ak)/2=(b−a)/2k+1|c_k - x^*| \le (b_k - a_k)/2 = (b-a)/2^{k+1} です。

誤差の上界が ff によらず決まることが二分法の強みです。許容誤差 τ\tau を達成する回数は k≥log⁡2((b−a)/τ)−1k \ge \log_2((b-a)/\tau) - 1 で、[0,1][0,1] から 10−1210^{-12} まで縮めるには約 40 回かかります。1 回で 1 ビットずつ、10 進では約 3.3 回で 1 桁ずつ進みます。Q 次数の言葉では、∣ck−x∗∣|c_k - x^*| そのものは規則正しく減るとは限りませんが、上界 (b−a)/2k+1(b-a)/2^{k+1} が比 1/21/2 で減るので、R 次数 1 で収束します。

例 3.2(二分法と Brent 法の反復回数)

f(x)=cos⁡x−xf(x) = \cos x - x を [0,1][0,1] で考えます。f(0)=1>0f(0) = 1 \gt 0、f(1)=cos⁡1−1<0f(1) = \cos 1 - 1 \lt 0 です。二分法と、SciPy の brentq(二分法・割線法・逆 2 次補間を組み合わせた Brent 法)で、区間幅の許容値ごとに反復回数を比べます。

# Python 3.12.10 / SciPy 1.18.1 / mpmath 1.3.0
import math
from scipy.optimize import brentq
from mpmath import mp, findroot, cos as mcos
def bisect(f, a, b, tol):
"""f(a)f(b) < 0。区間幅が tol 以下になるまで二分する"""
fa, n = f(a), 0
while b - a > tol:
c = (a + b) / 2
fc = f(c)
n += 1
if fc == 0:
return c, n
if (fa < 0) == (fc < 0):
a, fa = c, fc
else:
b = c
return (a + b) / 2, n
f = lambda x: math.cos(x) - x # 根は x = cos x の解
mp.dps = 30
ref = float(findroot(lambda x: mcos(x) - x, 0.7))
print(" tol 二分法: 回数 誤差 brentq: 回数 誤差")
for tol in [1e-4, 1e-8, 1e-12]:
xb, nb = bisect(f, 0.0, 1.0, tol)
xr, res = brentq(f, 0.0, 1.0, xtol=tol, full_output=True)
print(f"{tol:.0e} {nb:3d} {abs(xb - ref):.2e}"
f" {res.iterations:2d} {abs(xr - ref):.2e}")
print("理論の回数 ceil(log2(1/tol)):",
[math.ceil(math.log2(1 / t)) for t in [1e-4, 1e-8, 1e-12]])
tol 二分法: 回数 誤差 brentq: 回数 誤差
1e-04 14 2.01e-05 5 5.93e-06
1e-08 27 9.04e-10 6 1.35e-09
1e-12 40 5.17e-14 7 4.77e-15
理論の回数 ceil(log2(1/tol)): [14, 27, 40]

検算. (a) 参照値は mpmath の findroot で 30 桁まで求めた x=cos⁡xx = \cos x の解 0.739085…0.739085\ldots です。(b) 二分法の誤差はどれも許容値の半分以下で、定理 3.1 の 3 の評価(最後の区間の中点なので幅の半分以下)の内側にあります。Brent 法は別の算法で、同じ解に許容値以下の誤差で達しています。(c) 二分法の回数は幅 2−n≤τ2^{-n} \le \tau となる最小の n=⌈log⁡2(1/τ)⌉n = \lceil \log_2(1/\tau) \rceil に一致し、許容値を 10−410^{-4} 倍するごとに約 13.3 回ずつ増えます。Brent 法は 1 回ずつしか増えず、超 1 次収束していることが読み取れます。

方程式 f(x)=0f(x) = 0 を x=g(x)x = g(x) の形に書き直し、xk+1=g(xk)x_{k+1} = g(x_k) と反復するのが不動点反復です。g(x)=x−f(x)g(x) = x - f(x)、g(x)=x−f(x)/f′(x)g(x) = x - f(x)/f'(x) など書き直し方はいくらでもあり、どれを選ぶかで収束の速さが変わります。

定理 4.1(不動点反復の大域収束と誤差評価)✓

a<ba \lt b とし、g:[a,b]→[a,b]g : [a,b] \to [a,b] が定数 0≤L<10 \le L \lt 1 について

∣g(x)−g(y)∣≤L∣x−y∣(x,y∈[a,b])|g(x) - g(y)| \le L |x - y| \qquad (x, y \in [a,b])

を満たすとする。このとき gg は [a,b][a,b] にただ 1 つの不動点 x∗x^* をもち、任意の x0∈[a,b]x_0 \in [a,b] から xk+1=g(xk)x_{k+1} = g(x_k) で作った列は x∗x^* に収束して、k≥1k \ge 1 について次が成り立つ。

  1. (事前評価)∣xk−x∗∣≤Lk1−L∣x1−x0∣|x_k - x^*| \le \dfrac{L^k}{1-L}|x_1 - x_0|。
  2. (事後評価)∣xk−x∗∣≤L1−L∣xk−xk−1∣|x_k - x^*| \le \dfrac{L}{1-L}|x_k - x_{k-1}|。
  3. ∣xk−x∗∣≤Lk∣x0−x∗∣|x_k - x^*| \le L^k |x_0 - x^*|。

gg が [a,b][a,b] で微分可能で ∣g′(x)∣≤L|g'(x)| \le L なら、上の Lipschitz 条件が成り立つ。

証明(定理 4.1)

[a,b][a,b] は R\mathbb{R} の空でない閉集合で、gg はそれ自身への縮小率 LL の縮小写像です。縮小写像の原理(定理 3.2[逆関数定理・陰関数定理と条件付き極値]。n=1n = 1 の場合)から、不動点がただ 1 つ存在し、列は x∗x^* に収束して 1 が成り立ちます(一般の完備距離空間での形は 定理 6.2[距離空間])。

3 は ∣xj+1−x∗∣=∣g(xj)−g(x∗)∣≤L∣xj−x∗∣|x_{j+1} - x^*| = |g(x_j) - g(x^*)| \le L|x_j - x^*| を j=0,…,k−1j = 0, \ldots, k-1 について繰り返せば得られます。

2 は三角不等式と 3 の 1 段分から

∣xk−x∗∣≤∣xk−xk+1∣+∣xk+1−x∗∣=∣g(xk−1)−g(xk)∣+∣g(xk)−g(x∗)∣≤L∣xk−1−xk∣+L∣xk−x∗∣|x_k - x^*| \le |x_k - x_{k+1}| + |x_{k+1} - x^*| = |g(x_{k-1}) - g(x_k)| + |g(x_k) - g(x^*)| \le L|x_{k-1} - x_k| + L|x_k - x^*|

となり、∣xk−x∗∣|x_k - x^*| について解けば得られます。

最後の主張は、x≠yx \ne y について平均値の定理(定理 3.3[平均値の定理とテイラーの定理])から g(x)−g(y)=g′(c)(x−y)g(x) - g(y) = g'(c)(x - y) となる cc が xx と yy の間にあることから従います。

事前評価は反復を始める前に必要な回数を見積もるのに、事後評価は計算しながら止めどきを決めるのに使います。事後評価の右辺は計算済みの量 ∣xk−xk−1∣|x_k - x_{k-1}| だけで書けているので、「差が τ(1−L)/L\tau(1-L)/L 以下になったら止める」とすれば誤差は τ\tau 以下に保証されます。ただし LL が 11 に近いと係数 L/(1−L)L/(1-L) が大きくなり、差が小さくても誤差は小さいとは限りません。

例 4.2(x = cos x の不動点反復)

g(x)=cos⁡xg(x) = \cos x は [0,1][0,1] を [cos⁡1,1]⊂[0,1][\cos 1, 1] \subset [0,1] に写し、∣g′(x)∣=∣sin⁡x∣≤sin⁡1≈0.8415=L|g'(x)| = |\sin x| \le \sin 1 \approx 0.8415 = L です。定理 4.1 から不動点 x∗≈0.739085x^* \approx 0.739085 がただ 1 つあり、x0=0x_0 = 0 から反復すると x1=1x_1 = 1 です。事前評価で誤差 10−1010^{-10} を保証するには Lk/(1−L)≤10−10L^k/(1-L) \le 10^{-10}、すなわち k≥145k \ge 145 回が必要と見積もられます。実際に倍精度で反復すると、56 回目で誤差が 10−1010^{-10} を下回ります。差が出るのは、x∗x^* の近くでの縮小率 ∣g′(x∗)∣=sin⁡x∗≈0.674|g'(x^*)| = \sin x^* \approx 0.674 が区間全体での最悪値 0.84150.8415 より小さいためです。次の定理はこの「解の近くでの縮小率」を正確に述べます。

定理 4.3(不動点反復の局所収束と次数)✓

II を開区間、g:I→Rg : I \to \mathbb{R} を C1C^1 級、x∗∈Ix^* \in I を g(x∗)=x∗g(x^*) = x^* を満たす点とする。

  1. ∣g′(x∗)∣<1|g'(x^*)| \lt 1 ならば δ>0\delta \gt 0 があって、x0∈[x∗−δ,x∗+δ]x_0 \in [x^* - \delta, x^* + \delta] から作った xk+1=g(xk)x_{k+1} = g(x_k) はこの区間にとどまり、x∗x^* に収束する。さらにすべての kk で ek≠0e_k \ne 0 なら ek+1/ek→g′(x∗)e_{k+1}/e_k \to g'(x^*) である。とくに g′(x∗)≠0g'(x^*) \ne 0 なら 1 次収束で、漸近誤差定数は ∣g′(x∗)∣|g'(x^*)| である。
  2. p≥2p \ge 2 とし、gg が CpC^p 級で g′(x∗)=⋯=g(p−1)(x∗)=0g'(x^*) = \cdots = g^{(p-1)}(x^*) = 0 ならば、1 の δ\delta について、すべての kk で ek≠0e_k \ne 0 のとき ek+1ekp→g(p)(x∗)p!\frac{e_{k+1}}{e_k^p} \to \frac{g^{(p)}(x^*)}{p!} である。とくに g(p)(x∗)≠0g^{(p)}(x^*) \ne 0 なら pp 次収束で、漸近誤差定数は ∣g(p)(x∗)∣/p!|g^{(p)}(x^*)|/p! である。
証明(定理 4.3)

1. ∣g′(x∗)∣<L<1|g'(x^*)| \lt L \lt 1 となる LL を取ります。g′g' は x∗x^* で連続なので、δ>0\delta \gt 0 を小さく取れば J=[x∗−δ,x∗+δ]⊂IJ = [x^* - \delta, x^* + \delta] \subset I かつ JJ 上で ∣g′∣≤L|g'| \le L にできます。x∈Jx \in J、x≠x∗x \ne x^* なら平均値の定理(定理 3.3[平均値の定理とテイラーの定理])から xx と x∗x^* の間の η\eta があって

g(x)−x∗=g(x)−g(x∗)=g′(η)(x−x∗)g(x) - x^* = g(x) - g(x^*) = g'(\eta)(x - x^*)

です。η∈J\eta \in J なので ∣g(x)−x∗∣≤L∣x−x∗∣≤δ|g(x) - x^*| \le L|x - x^*| \le \delta、すなわち g(J)⊂Jg(J) \subset J です。よって列は JJ にとどまり、∣ek∣≤Lk∣e0∣→0|e_k| \le L^k |e_0| \to 0 です。ek≠0e_k \ne 0 のとき上の式で x=xkx = x_k とすると ek+1/ek=g′(ηk)e_{k+1}/e_k = g'(\eta_k) で、ηk\eta_k は xkx_k と x∗x^* の間にあるので ηk→x∗\eta_k \to x^* です。g′g' の連続性から g′(ηk)→g′(x∗)g'(\eta_k) \to g'(x^*) です。

2. g′(x∗)=0g'(x^*) = 0 なので 1 が適用でき、列は JJ にとどまって x∗x^* に収束します。ek≠0e_k \ne 0 のとき、Taylor の定理(定理 5.3[平均値の定理とテイラーの定理]。a=x∗a = x^*、n=p−1n = p - 1)を x=xkx = x_k に使うと、xkx_k と x∗x^* の間の ξk\xi_k があって

g(xk)=g(x∗)+∑j=1p−1g(j)(x∗)j!ekj+g(p)(ξk)p!ekp=x∗+g(p)(ξk)p!ekpg(x_k) = g(x^*) + \sum_{j=1}^{p-1} \frac{g^{(j)}(x^*)}{j!} e_k^j + \frac{g^{(p)}(\xi_k)}{p!} e_k^p = x^* + \frac{g^{(p)}(\xi_k)}{p!} e_k^p

です。よって ek+1/ekp=g(p)(ξk)/p!e_{k+1}/e_k^p = g^{(p)}(\xi_k)/p! で、ξk→x∗\xi_k \to x^* と g(p)g^{(p)} の連続性から主張を得ます。

この定理は「反復関数の導関数が解で消えるほど速い」ことを述べています。Newton 法はちょうど g′(x∗)=0g'(x^*) = 0 となるように gg を設計した反復と見ることができ、次の節の 2 次収束はその具体化です。

同じ方程式 x=cos⁡xx = \cos x を 2 通りに書き直した反復で、∣g′(x∗)∣|g'(x^*)| が 1 より小さいか大きいかの違いを見ます。左は g=cos⁡g = \cos で、折れ線が不動点に巻き込むように近づきます。右は逆関数 g=arccos⁡g = \arccos で、同じ不動点から離れていきます。

y = g(x) と y = x のグラフの上で、反復の折れ線が階段状に進む動画。左の g = cos x では折れ線が不動点のまわりを渦巻きながら近づき、右の g = arccos x では同じ不動点のまわりを渦巻きながら外へ広がっていく。
x = cos x の不動点 x* ≈ 0.7391(星印)に向かう反復。左は |g'(x*)| ≈ 0.674 で収束し、右は |g'(x*)| ≈ 1.484 で、初期値 0.70 が x* のすぐ近くでも離れていく。

xkx_k での ff の接線 y=f(xk)+f′(xk)(x−xk)y = f(x_k) + f'(x_k)(x - x_k) と xx 軸の交点を次の近似とするのが Newton 法

xk+1=xk−f(xk)f′(xk)x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}

です。ff を xkx_k のまわりで 1 次式に置き換えて解く、と言っても同じです。置き換えの誤差は Taylor の定理の 2 次の剰余なので、誤差が 2 乗で効くことが予想できます。

x₀ = 2x₁ = 1.5x₂√2y = x² − 2
凸関数 f(x) = x² − 2 での Newton 法。x₀ = 2 での接線が x 軸と交わる点が x₁ = 1.5、x₁ での接線から x₂ ≈ 1.4167 が得られ、解 √2 ≈ 1.4142 に右から近づく。

次の定理は文字が多いですが、言っていることは 1 つです。解の近くの範囲 JJ で、mm は ∣f′∣|f'| の最小値、MM は ∣f′′∣|f''| の最大値、K=M/(2m)K = M/(2m) はそこから決まる 1 つの定数です。Kδ<1K\delta \lt 1 は「範囲の幅が十分狭い」という条件で、これが成り立てば 1 回の反復で誤差が必ず小さくなり、2 乗の形で減っていきます。

定理 5.1(Newton 法の局所 2 次収束)✓

II を開区間、f:I→Rf : I \to \mathbb{R} を C2C^2 級とし、x∗∈Ix^* \in I で f(x∗)=0f(x^*) = 0、f′(x∗)≠0f'(x^*) \ne 0 とする。

  1. δ>0\delta \gt 0 が J=[x∗−δ,x∗+δ]⊂IJ = [x^* - \delta, x^* + \delta] \subset I と m=min⁡J∣f′∣>0m = \min_J |f'| \gt 0 を満たすとし、M=max⁡J∣f′′∣M = \max_J |f''|、K=M/(2m)K = M/(2m) とおく。Kδ<1K\delta \lt 1 ならば、任意の x0∈Jx_0 \in J について Newton 法の列は定義されて JJ にとどまり、 ∣ek+1∣≤Kek2,K∣ek∣≤(K∣e0∣)2k|e_{k+1}| \le K e_k^2, \qquad K|e_k| \le (K|e_0|)^{2^k} が成り立つ。とくに xk→x∗x_k \to x^* で、R 次数 2 以上で収束する。
  2. 1 の条件を満たす δ>0\delta \gt 0 は存在する。
  3. 1 の状況で、すべての kk で ek≠0e_k \ne 0 ならば ek+1ek2→f′′(x∗)2f′(x∗)\frac{e_{k+1}}{e_k^2} \to \frac{f''(x^*)}{2f'(x^*)} である。とくに f′′(x∗)≠0f''(x^*) \ne 0 なら 2 次収束で、漸近誤差定数は ∣f′′(x∗)∣/(2∣f′(x∗)∣)|f''(x^*)|/(2|f'(x^*)|) である。
証明(定理 5.1)

mm と MM は、連続関数 ∣f′∣|f'|、∣f′′∣|f''| の閉区間 JJ 上の最小値・最大値として存在します(定理 2.3[平均値の定理とテイラーの定理])。

1. xk∈Jx_k \in J と仮定します。∣f′(xk)∣≥m>0|f'(x_k)| \ge m \gt 0 なので xk+1x_{k+1} は定義されます。xk=x∗x_k = x^* なら f(xk)=0f(x_k) = 0 で xk+1=x∗x_{k+1} = x^* となり、以後の主張は自明です。xk≠x∗x_k \ne x^* なら、Taylor の定理(定理 5.3[平均値の定理とテイラーの定理]。a=xka = x_k、n=1n = 1)を x=x∗x = x^* に使うと、xkx_k と x∗x^* の間の ξk\xi_k があって

0=f(x∗)=f(xk)+f′(xk)(x∗−xk)+f′′(ξk)2(x∗−xk)20 = f(x^*) = f(x_k) + f'(x_k)(x^* - x_k) + \frac{f''(\xi_k)}{2}(x^* - x_k)^2

です。両辺を f′(xk)f'(x_k) で割って整理すると

ek+1=xk−f(xk)f′(xk)−x∗=f′′(ξk)2f′(xk) ek2e_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} - x^* = \frac{f''(\xi_k)}{2f'(x_k)}\, e_k^2

を得ます。ξk∈J\xi_k \in J なので ∣ek+1∣≤Kek2≤(Kδ)∣ek∣<∣ek∣≤δ|e_{k+1}| \le K e_k^2 \le (K\delta)|e_k| \lt |e_k| \le \delta となり、xk+1∈Jx_{k+1} \in J です。帰納法で列は JJ にとどまり、∣ek+1∣≤Kek2|e_{k+1}| \le K e_k^2 がすべての kk で成り立ちます。両辺に KK を掛けると K∣ek+1∣≤(K∣ek∣)2K|e_{k+1}| \le (K|e_k|)^2 で、これを繰り返すと K∣ek∣≤(K∣e0∣)2kK|e_k| \le (K|e_0|)^{2^k} です。K∣e0∣≤Kδ<1K|e_0| \le K\delta \lt 1 なので右辺は 00 に収束します(K=0K = 0 なら e1=0e_1 = 0 です)。K>0K \gt 0 なら A=1/KA = 1/K、q=Kδq = K\delta として R 次数 2 以上の定義(定義 2.2)を満たします。

2. f′f' は x∗x^* で連続で f′(x∗)≠0f'(x^*) \ne 0 なので、δ0>0\delta_0 \gt 0 を J0=[x∗−δ0,x∗+δ0]⊂IJ_0 = [x^* - \delta_0, x^* + \delta_0] \subset I かつ J0J_0 上で ∣f′∣≥∣f′(x∗)∣/2|f'| \ge |f'(x^*)|/2 となるように取れます。M0=max⁡J0∣f′′∣M_0 = \max_{J_0} |f''| とします。0<δ≤δ00 \lt \delta \le \delta_0 なら m≥∣f′(x∗)∣/2m \ge |f'(x^*)|/2、M≤M0M \le M_0 なので K≤M0/∣f′(x∗)∣K \le M_0/|f'(x^*)| です。よって δ<∣f′(x∗)∣/M0\delta \lt |f'(x^*)|/M_0(M0=0M_0 = 0 なら任意)かつ δ≤δ0\delta \le \delta_0 と取れば Kδ<1K\delta \lt 1 です。

3. 1 の式から ek+1/ek2=f′′(ξk)/(2f′(xk))e_{k+1}/e_k^2 = f''(\xi_k)/(2f'(x_k)) です。xk→x∗x_k \to x^* で、ξk\xi_k は xkx_k と x∗x^* の間にあるので ξk→x∗\xi_k \to x^* です。f′f'、f′′f'' の連続性と f′(x∗)≠0f'(x^*) \ne 0 から主張が従います。

1 の評価は、解を知らなくても使える形にはなっていません(JJ が解を中心にしているため)。それでも「K∣e0∣<1K|e_0| \lt 1 なら収束し、正しい桁数がほぼ倍々に増える」という見通しを与え、KK が大きい(f′f' が小さい、または f′′f'' が大きい)ほど収束域が狭いことを教えます。解を含む球を計算だけで保証する形は、後の節の Kantorovich 型の定理(定理 7.4)で扱います。

例 5.2(√2 の Newton 法)

f(x)=x2−2f(x) = x^2 - 2 では Newton 法は xk+1=(xk+2/xk)/2x_{k+1} = (x_k + 2/x_k)/2 です。定理 5.1 の 3 が予言する極限は f′′(x∗)/(2f′(x∗))=2/(42)=1/(22)≈0.35355f''(x^*)/(2f'(x^*)) = 2/(4\sqrt2) = 1/(2\sqrt2) \approx 0.35355 です。倍精度では誤差がすぐ丸めの水準(10−1610^{-16})で止まってしまうので、mpmath の 60 桁で反復して 2 次収束を最後まで追います。

# Python 3.12.10 / mpmath 1.3.0
from mpmath import mp, mpf, sqrt
mp.dps = 60 # 10 進 60 桁で反復する
r = sqrt(mpf(2)) # 参照値 √2
x, xd = mpf(2), 2.0 # x0 = 2(xd は倍精度で同じ反復)
ep = None
print(" k e_k(60桁) e_k/e_{k-1}^2 e_k(倍精度)")
for k in range(7):
e, ed = abs(x - r), abs(mpf(xd) - r)
q = f"{float(e / ep**2):.5f}" if ep else " - "
print(f"{k:2d} {float(e):.2e} {q} {float(ed):.2e}")
ep = e
x = x - (x * x - 2) / (2 * x) # Newton 法
xd = xd - (xd * xd - 2) / (2 * xd)
print("1/(2√2) =", f"{float(1 / (2 * r)):.5f}")
k e_k(60桁) e_k/e_{k-1}^2 e_k(倍精度)
0 5.86e-01 - 5.86e-01
1 8.58e-02 0.25000 8.58e-02
2 2.45e-03 0.33333 2.45e-03
3 2.12e-06 0.35294 2.12e-06
4 1.59e-12 0.35355 1.59e-12
5 8.99e-25 0.35355 9.67e-17
6 2.86e-49 0.35355 1.25e-16
1/(2√2) = 0.35355

検算. (a) この例では証明の式が閉じた形になり、ek+1=ek2/(2xk)e_{k+1} = e_k^2/(2x_k) です。表の比 0.25000=1/(2⋅2)0.25000 = 1/(2 \cdot 2)、0.33333=1/(2⋅1.5)0.33333 = 1/(2 \cdot 1.5) はちょうど 1/(2xk−1)1/(2x_{k-1}) に一致します。(b) 60 桁の計算と倍精度の計算は k=4k = 4 まで表示の桁で一致し、倍精度は k=5k = 5 以降、2\sqrt2 に最も近い倍精度数の誤差 9.67×10−179.67 \times 10^{-17} の近くで止まります(この先は丸め誤差の問題で、Newton 法の収束の問題ではありません)。(c) 誤差の指数は −1,−3,−6,−12,−25,−49-1, -3, -6, -12, -25, -49 とほぼ倍々に増え、比 ek/ek−12e_k/e_{k-1}^2 は理論値 1/(22)1/(2\sqrt2) に近づきます。

定理 5.3(凸関数での Newton 法の単調な大域収束)✓

II を区間、f:I→Rf : I \to \mathbb{R} を微分可能で f′f' が広義単調増加(ff が凸)な関数とし、x∗∈Ix^* \in I で f(x∗)=0f(x^*) = 0、f′(x∗)>0f'(x^*) \gt 0 とする。x0∈Ix_0 \in I、x0>x∗x_0 \gt x^* ならば、Newton 法の列は定義されて

x∗≤xk+1≤xk(k≥0)x^* \le x_{k+1} \le x_k \qquad (k \ge 0)

を満たし、x∗x^* に収束する。

証明(定理 5.3)

xk∈Ix_k \in I、xk>x∗x_k \gt x^* と仮定します。f′f' は増加なので f′(xk)≥f′(x∗)>0f'(x_k) \ge f'(x^*) \gt 0 で、xk+1x_{k+1} は定義されます。平均値の定理(定理 3.3[平均値の定理とテイラーの定理])から f(xk)=f(xk)−f(x∗)=f′(c)(xk−x∗)f(x_k) = f(x_k) - f(x^*) = f'(c)(x_k - x^*) となる c∈(x∗,xk)c \in (x^*, x_k) があり、0<f′(x∗)≤f′(c)≤f′(xk)0 \lt f'(x^*) \le f'(c) \le f'(x_k) です。よって

xk−xk+1=f(xk)f′(xk)=f′(c)f′(xk)(xk−x∗)∈(0, xk−x∗]x_k - x_{k+1} = \frac{f(x_k)}{f'(x_k)} = \frac{f'(c)}{f'(x_k)}(x_k - x^*) \in (0,\, x_k - x^*]

で、x∗≤xk+1<xkx^* \le x_{k+1} \lt x_k です。xk+1=x∗x_{k+1} = x^* なら以後 xj=x∗x_j = x^* で主張は自明です。そうでなければ xk+1∈(x∗,x0]⊂Ix_{k+1} \in (x^*, x_0] \subset I で同じ議論を繰り返せます。

列 (−xk)(-x_k) は単調増加で −x∗-x^* で上に有界なので、単調有界定理(定理 2.1[級数と収束判定])から xkx_k はある ℓ≥x∗\ell \ge x^* に収束します。f′(xk)f'(x_k) は [f′(x∗),f′(x0)][f'(x^*), f'(x_0)] に入るので、f(xk)=(xk−xk+1)f′(xk)f(x_k) = (x_k - x_{k+1}) f'(x_k) から ∣f(xk)∣≤f′(x0)∣xk−xk+1∣→0|f(x_k)| \le f'(x_0)|x_k - x_{k+1}| \to 0 です。ff は連続なので(定理 4.5[極限と連続性])f(ℓ)=0f(\ell) = 0 です。一方 x>x∗x \gt x^* なら上と同じく平均値の定理から f(x)=f′(c)(x−x∗)>0f(x) = f'(c)(x - x^*) \gt 0 なので、ℓ=x∗\ell = x^* です。

f(x)=x2−2f(x) = x^2 - 2 は凸で、x0>2x_0 \gt \sqrt2 ならどこから始めても単調に減少して 2\sqrt2 に収束します。0<x0<20 \lt x_0 \lt \sqrt2 から始めても、1 回目で x1=(x0+2/x0)/2≥2x_1 = (x_0 + 2/x_0)/2 \ge \sqrt2(相加相乗平均)となるので、2 回目以降はこの定理の状況に入ります。

f′(x∗)=0f'(x^*) = 0 のとき、定理 5.1 は使えません。重根では Newton 法は収束はしますが、1 次収束に落ちます。

次の命題の hh は、解 x∗x^* を mm 回掛けて取り出した残りの部分です。h(x∗)≠0h(x^*) \ne 0 なので、解の近くでは ff は (x−x∗)m(x - x^*)^m とほぼ同じ形に振る舞います。たとえば (x−1)2(x-1)^2 なら m=2m = 2、h≡1h \equiv 1 です。

命題 5.4(重根での Newton 法と修正 Newton 法)✓

II を開区間、x∗∈Ix^* \in I、m≥2m \ge 2 を整数とし、f(x)=(x−x∗)mh(x)f(x) = (x - x^*)^m h(x)(h:I→Rh : I \to \mathbb{R} は C1C^1 級、h(x∗)≠0h(x^*) \ne 0)とする。

  1. δ>0\delta \gt 0 があって、x0∈[x∗−δ,x∗+δ]x_0 \in [x^* - \delta, x^* + \delta] から Newton 法の列は定義されて x∗x^* に収束し、すべての kk で ek≠0e_k \ne 0 なら ek+1/ek→(m−1)/me_{k+1}/e_k \to (m-1)/m である。すなわち 1 次収束で、漸近誤差定数は (m−1)/m(m-1)/m である。
  2. 重複度 mm を掛けた修正 Newton 法 xk+1=xk−mf(xk)/f′(xk)x_{k+1} = x_k - m f(x_k)/f'(x_k) では、δ\delta を小さく取り直せば列は x∗x^* に収束し、すべての kk で ek≠0e_k \ne 0 なら ek+1ek2→h′(x∗)m h(x∗)\frac{e_{k+1}}{e_k^2} \to \frac{h'(x^*)}{m\, h(x^*)} である。
証明(命題 5.4)

積の微分から f′(x)=(x−x∗)m−1D(x)f'(x) = (x - x^*)^{m-1} D(x)、D(x)=mh(x)+(x−x∗)h′(x)D(x) = m h(x) + (x - x^*) h'(x) です。DD は連続で D(x∗)=mh(x∗)≠0D(x^*) = m h(x^*) \ne 0 なので、δ0>0\delta_0 \gt 0 を J0=[x∗−δ0,x∗+δ0]⊂IJ_0 = [x^* - \delta_0, x^* + \delta_0] \subset I 上で DD も hh も 00 にならないように取れます。このとき x∈J0x \in J_0、x≠x∗x \ne x^* で f′(x)≠0f'(x) \ne 0 であり、

f(x)f′(x)=(x−x∗) φ(x),φ(x)=h(x)D(x)\frac{f(x)}{f'(x)} = (x - x^*)\,\varphi(x), \qquad \varphi(x) = \frac{h(x)}{D(x)}

と書けます。φ\varphi は J0J_0 上連続で φ(x∗)=1/m\varphi(x^*) = 1/m です。

1. Newton 法では ek+1=ek−ekφ(xk)=ek(1−φ(xk))e_{k+1} = e_k - e_k\varphi(x_k) = e_k(1 - \varphi(x_k)) です(ek=0e_k = 0 なら止める)。1−φ1 - \varphi は連続で 1−φ(x∗)=(m−1)/m<11 - \varphi(x^*) = (m-1)/m \lt 1 なので、(m−1)/m<L<1(m-1)/m \lt L \lt 1 を取り、δ≤δ0\delta \le \delta_0 を [x∗−δ,x∗+δ][x^* - \delta, x^* + \delta] 上で ∣1−φ∣≤L|1 - \varphi| \le L となるように取れます。すると ∣ek+1∣≤L∣ek∣|e_{k+1}| \le L|e_k| で列はこの区間にとどまり、ek→0e_k \to 0 です。さらに ek+1/ek=1−φ(xk)→1−1/me_{k+1}/e_k = 1 - \varphi(x_k) \to 1 - 1/m です。

2. 修正 Newton 法では ek+1=ek(1−mφ(xk))e_{k+1} = e_k(1 - m\varphi(x_k)) で、

1−mφ(x)=D(x)−mh(x)D(x)=(x−x∗) h′(x)D(x)1 - m\varphi(x) = \frac{D(x) - m h(x)}{D(x)} = \frac{(x - x^*)\, h'(x)}{D(x)}

です。よって ek+1=ek2 h′(xk)/D(xk)e_{k+1} = e_k^2\, h'(x_k)/D(x_k) です。ψ=h′/D\psi = h'/D は J0J_0 上連続なので ∣ψ∣≤N|\psi| \le N と有界です(定理 2.3[平均値の定理とテイラーの定理])。δ≤δ0\delta \le \delta_0 を Nδ<1N\delta \lt 1 となるように取れば、∣ek+1∣≤Nδ∣ek∣|e_{k+1}| \le N\delta|e_k| で列は区間にとどまり ek→0e_k \to 0 です。最後に ek+1/ek2=ψ(xk)→h′(x∗)/D(x∗)=h′(x∗)/(mh(x∗))e_{k+1}/e_k^2 = \psi(x_k) \to h'(x^*)/D(x^*) = h'(x^*)/(m h(x^*)) です。

例 5.5((x − 1)² での遅い収束)

f(x)=(x−1)2f(x) = (x-1)^2 では Newton 法は xk+1=xk−(xk−1)/2x_{k+1} = x_k - (x_k - 1)/2、すなわち ek+1=ek/2e_{k+1} = e_k/2 です。誤差は 1 回で半分にしかならず、二分法と同じ速さです。命題 5.4 の比 (m−1)/m=1/2(m-1)/m = 1/2 がちょうど実現しています。修正 Newton 法 xk+1=xk−2f(xk)/f′(xk)x_{k+1} = x_k - 2f(x_k)/f'(x_k) では 1 回で x1=1x_1 = 1 に到達します(h≡1h \equiv 1 で h′=0h' = 0)。実際の問題では重複度 mm は未知のことが多く、ek+1/eke_{k+1}/e_k の観測値から mm を推定するか、f/f′f/f'(単根しかもたない)に Newton 法を使う工夫をします。

重根は収束が遅いだけでなく、達成できる精度そのものが低い根でもあります。ff の値の計算に含まれる誤差を ff に加わる小さな摂動 εq\varepsilon q とみなすと、多項式の単根は 1 次で −εq(x∗)/f′(x∗)-\varepsilon q(x^*)/f'(x^*) だけ動きます(命題 3.1[条件数と安定性])。∣f′(x∗)∣|f'(x^*)| が小さい根、つまり重根に近い根ほど、反復法をどう選んでも計算できる桁数は少なくなります。これは算法ではなく問題の性質(条件数)で、条件数と安定性の章で扱います。

例 5.6(周期軌道に入る Newton 法)

f(x)=x3−2x+2f(x) = x^3 - 2x + 2 に x0=0x_0 = 0 から Newton 法を使います。f(0)=2f(0) = 2、f′(0)=−2f'(0) = -2 なので x1=1x_1 = 1、f(1)=1f(1) = 1、f′(1)=1f'(1) = 1 なので x2=0x_2 = 0 です。列は 0,1,0,1,…0, 1, 0, 1, \ldots と周期 2 で回り続け、実数解 x∗≈−1.769x^* \approx -1.769 には近づきません。f′f' は ±2/3\pm\sqrt{2/3} で 00 になり、その間では接線が大きく外れた点に飛ばします。局所収束の定理の「解の十分近くから始めれば」という仮定は飾りではなく、実用では二分法で区間を絞ってから Newton 法に切り替える、減速(xk+1=xk−λkf(xk)/f′(xk)x_{k+1} = x_k - \lambda_k f(x_k)/f'(x_k)、0<λk≤10 \lt \lambda_k \le 1)を入れる、などの対策をとります。

Newton 法は毎回 f′(xk)f'(x_k) を必要とします。導関数の式がない、あるいは計算が高くつくときは、接線の傾きを直前の 2 点を結ぶ直線の傾きで置き換えます。これが割線法

xk+1=xk−f(xk) xk−xk−1f(xk)−f(xk−1)x_{k+1} = x_k - f(x_k)\,\frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})}

です。1 回の反復で新しく計算する関数値は f(xk)f(x_k) の 1 つだけです。解析には差分商を使います。f[a,b]f[a,b] は 2 点 a,ba, b を結ぶ線分の傾き、f[a,b,c]f[a,b,c] はその傾きの変化の割合で、導関数と 2 階導関数の「差分版」です。相異なる点 a,b,ca, b, c について

f[a,b]=f(b)−f(a)b−a,f[a,b,c]=f[a,b]−f[b,c]a−cf[a,b] = \frac{f(b) - f(a)}{b - a}, \qquad f[a,b,c] = \frac{f[a,b] - f[b,c]}{a - c}

と書きます。定義から f[a,b]=f[b,a]f[a,b] = f[b,a] です。割線法は xk+1=xk−f(xk)/f[xk−1,xk]x_{k+1} = x_k - f(x_k)/f[x_{k-1}, x_k] と書けます。

次の補題のねらいは、割線法の誤差の式に出てくる差分商を、導関数の値に置き換えることです。1 つ目は平均値の定理の言い換えで、2 つ目は 3 点を通る 2 次式を作り、Rolle の定理を 2 回使って示します。

補題 6.1(2 階の差分商の平均値定理)✓

II を開区間、f:I→Rf : I \to \mathbb{R} を 2 回微分可能とし、a,b,c∈Ia, b, c \in I を相異なる点とする。このとき aa と bb の間の η\eta と、a,b,ca, b, c の最小値と最大値の間の ξ\xi があって

f[a,b]=f′(η),f[a,b,c]=f′′(ξ)2f[a,b] = f'(\eta), \qquad f[a,b,c] = \frac{f''(\xi)}{2}

が成り立つ。

証明(補題 6.1)

1 つ目は平均値の定理(定理 3.3[平均値の定理とテイラーの定理])そのものです。2 つ目のために 2 次式

p(x)=f(a)+f[a,b](x−a)+f[a,b,c](x−a)(x−b)p(x) = f(a) + f[a,b](x - a) + f[a,b,c](x - a)(x - b)

を考えます。p(a)=f(a)p(a) = f(a)、p(b)=f(a)+f[a,b](b−a)=f(b)p(b) = f(a) + f[a,b](b-a) = f(b) です。また定義から f[a,b,c](c−a)=f[b,c]−f[a,b]f[a,b,c](c - a) = f[b,c] - f[a,b] なので

f[a,b,c](c−a)(c−b)=(f[b,c]−f[a,b])(c−b)=f(c)−f(b)−f[a,b](c−b)f[a,b,c](c-a)(c-b) = \bigl(f[b,c] - f[a,b]\bigr)(c - b) = f(c) - f(b) - f[a,b](c - b)

となり、p(c)=f(a)+f[a,b](c−a)+f(c)−f(b)−f[a,b](c−b)=f(a)+f[a,b](b−a)+f(c)−f(b)=f(c)p(c) = f(a) + f[a,b](c - a) + f(c) - f(b) - f[a,b](c-b) = f(a) + f[a,b](b - a) + f(c) - f(b) = f(c) です。

g=f−pg = f - p は II 上 2 回微分可能で、3 点 a,b,ca, b, c で 00 になります。3 点を小さい順に t1<t2<t3t_1 \lt t_2 \lt t_3 とすると、Rolle の定理(定理 3.1[平均値の定理とテイラーの定理])を [t1,t2][t_1, t_2] と [t2,t3][t_2, t_3] に使って g′(s1)=g′(s2)=0g'(s_1) = g'(s_2) = 0 となる t1<s1<t2<s2<t3t_1 \lt s_1 \lt t_2 \lt s_2 \lt t_3 を得ます。g′g' は微分可能なので連続で、もう一度 Rolle の定理を [s1,s2][s_1, s_2] に使うと g′′(ξ)=0g''(\xi) = 0 となる ξ∈(s1,s2)⊂(t1,t3)\xi \in (s_1, s_2) \subset (t_1, t_3) があります。p′′=2f[a,b,c]p'' = 2f[a,b,c] なので f′′(ξ)=2f[a,b,c]f''(\xi) = 2f[a,b,c] です。

一般の nn 階の差分商についての同じ形の定理は、多項式補間の章で補間多項式を使って示します。ここでは割線法に必要な 2 階の場合だけを使います。

定理 6.2(割線法の局所収束)✓

II を開区間、f:I→Rf : I \to \mathbb{R} を C2C^2 級とし、x∗∈Ix^* \in I で f(x∗)=0f(x^*) = 0、f′(x∗)≠0f'(x^*) \ne 0 とする。δ>0\delta \gt 0 が J=[x∗−δ,x∗+δ]⊂IJ = [x^* - \delta, x^* + \delta] \subset I と m=min⁡J∣f′∣>0m = \min_J |f'| \gt 0 を満たし、M=max⁡J∣f′′∣M = \max_J |f''| に対して K≥M/(2m)K \ge M/(2m) かつ Kδ<1K\delta \lt 1 とする(このような δ\delta と KK は存在する)。相異なる x0,x1∈Jx_0, x_1 \in J から割線法を始め、f(xk)=0f(x_k) = 0 となったら止める。このとき

  1. 列は定義されて JJ にとどまり、∣ek+1∣≤K∣ek∣∣ek−1∣|e_{k+1}| \le K|e_k||e_{k-1}| を満たして x∗x^* に収束する。
  2. 止まらないならば、 ek+1ek ek−1→f′′(x∗)2f′(x∗)\frac{e_{k+1}}{e_k\, e_{k-1}} \to \frac{f''(x^*)}{2f'(x^*)} である。
証明(定理 6.2)

δ\delta と KK の存在は、定理 5.1 の 2 と同じ議論で δ\delta を取り、K=M/(2m)K = M/(2m) とすれば得られます。

1. xk−1,xk∈Jx_{k-1}, x_k \in J が相異なり、f(xk−1)≠0f(x_{k-1}) \ne 0、f(xk)≠0f(x_k) \ne 0 と仮定します(k=1k = 1 では仮定どおりで、止まっていなければ以下で帰納的に保たれます)。補題 6.1 から f[xk−1,xk]=f′(η)f[x_{k-1}, x_k] = f'(\eta)(η∈J\eta \in J)で、∣f′(η)∣≥m>0|f'(\eta)| \ge m \gt 0 なので xk+1x_{k+1} は定義されます。f(xk)=f(xk)−f(x∗)=f[xk,x∗] ekf(x_k) = f(x_k) - f(x^*) = f[x_k, x^*]\, e_k なので

ek+1=ek−f[xk,x∗]f[xk−1,xk] ek=ek f[xk−1,xk]−f[xk,x∗]f[xk−1,xk]=ek ek−1 f[xk−1,xk,x∗]f[xk−1,xk]e_{k+1} = e_k - \frac{f[x_k, x^*]}{f[x_{k-1}, x_k]}\, e_k = e_k\,\frac{f[x_{k-1}, x_k] - f[x_k, x^*]}{f[x_{k-1}, x_k]} = e_k\, e_{k-1}\,\frac{f[x_{k-1}, x_k, x^*]}{f[x_{k-1}, x_k]}

です。最後の等号は 2 階の差分商の定義そのものです(3 点 xk−1,xk,x∗x_{k-1}, x_k, x^* は相異なります。f(xk−1)≠0f(x_{k-1}) \ne 0、f(xk)≠0f(x_k) \ne 0 だからです)。補題 6.1 から分子は f′′(ξk)/2f''(\xi_k)/2(ξk∈J\xi_k \in J)なので、∣ek+1∣≤M2m∣ek∣∣ek−1∣≤K∣ek∣∣ek−1∣|e_{k+1}| \le \frac{M}{2m}|e_k||e_{k-1}| \le K|e_k||e_{k-1}| です。∣ek−1∣≤δ|e_{k-1}| \le \delta なので ∣ek+1∣≤Kδ∣ek∣<∣ek∣|e_{k+1}| \le K\delta|e_k| \lt |e_k| で、xk+1∈Jx_{k+1} \in J です。さらに xk+1=xkx_{k+1} = x_k となるのは f(xk)=0f(x_k) = 0 のときだけなので、xk≠xk+1x_k \ne x_{k+1} です。帰納法ですべての段で成り立ち、∣ek+1∣≤(Kδ)k∣e1∣→0|e_{k+1}| \le (K\delta)^k |e_1| \to 0 です。

2. 1 の式から ek+1/(ekek−1)=f′′(ξk)/(2f′(ηk))e_{k+1}/(e_k e_{k-1}) = f''(\xi_k)/(2f'(\eta_k)) です。ηk\eta_k と ξk\xi_k はいずれも xk−1,xk,x∗x_{k-1}, x_k, x^* の最小値と最大値の間にあるので x∗x^* に収束し、f′f'、f′′f'' の連続性から主張が従います。

系 6.3(割線法の R 次数は黄金比以上)✓

定理 6.2 の状況で割線法が止まらないならば、φ=(1+5)/2\varphi = (1 + \sqrt5)/2 について、列は R 次数 φ\varphi 以上で x∗x^* に収束する。

証明の方針は、dk+1≤dkdk−1d_{k+1} \le d_k d_{k-1} という「前の 2 回の積」の形の不等式から、指数がどう増えるかを調べることです。指数が φk\varphi^k のとき、φk+φk−1=φk+1\varphi^k + \varphi^{k-1} = \varphi^{k+1} となって指数の足し算がぴったり合います。これが黄金比が出てくる理由で、qq の取り方は最初の 2 項でこの形を満たすためのものです。

証明(系 6.3)

K~=max⁡(K,1/(2δ))\tilde K = \max(K, 1/(2\delta)) とおくと K~>0\tilde K \gt 0、K~δ<1\tilde K\delta \lt 1 で、定理 6.2 の 1 の不等式は KK を K~\tilde K に替えても成り立ちます。dk=K~∣ek∣d_k = \tilde K|e_k| とおくと dk+1≤dkdk−1d_{k+1} \le d_k d_{k-1}、0<dk≤K~δ<10 \lt d_k \le \tilde K\delta \lt 1 です。q=max⁡(d0,d11/φ)q = \max(d_0, d_1^{1/\varphi}) とすると 0<q<10 \lt q \lt 1 で、d0≤q=qφ0d_0 \le q = q^{\varphi^0}、d1≤qφd_1 \le q^{\varphi} です。dk−1≤qφk−1d_{k-1} \le q^{\varphi^{k-1}}、dk≤qφkd_k \le q^{\varphi^k} なら、φ2=φ+1\varphi^2 = \varphi + 1 を使って

dk+1≤qφk+φk−1=qφk−1(φ+1)=qφk+1d_{k+1} \le q^{\varphi^k + \varphi^{k-1}} = q^{\varphi^{k-1}(\varphi + 1)} = q^{\varphi^{k+1}}

です。帰納法で ∣ek∣≤K~−1qφk|e_k| \le \tilde K^{-1} q^{\varphi^k} がすべての kk で成り立ち、R 次数 φ\varphi 以上の定義(定義 2.2)を満たします。

1 回あたりの関数値の計算が Newton 法では ff と f′f' の 2 回、割線法では ff の 1 回です。f′f' の計算が ff と同じ手間なら、割線法の 2 回分(次数 φ2≈2.618\varphi^2 \approx 2.618)は Newton 法の 1 回分(次数 2)より速いことになります。

注意 6.4(割線法の Q 次数)

f′′(x∗)≠0f''(x^*) \ne 0 なら、C=∣f′′(x∗)/(2f′(x∗))∣C = |f''(x^*)/(2f'(x^*))| として ∣ek+1∣/∣ek∣φ→Cφ−1|e_{k+1}|/|e_k|^{\varphi} \to C^{\varphi - 1} が成り立ち、割線法は Q 次数の意味でもちょうど φ\varphi 次収束します。定理 6.2 の 2 の漸化式の対数をとって線形差分方程式として解析する証明が知られています(Atkinson『An Introduction to Numerical Analysis』第 2 章の割線法の節、Ortega–Rheinboldt『Iterative Solution of Nonlinear Equations in Several Variables』の収束の速さを扱う章)。本章では R 次数の評価だけを使います。

例 6.5(割線法の次数を誤差から推定する)

f(x)=x2−2f(x) = x^2 - 2、x0=2x_0 = 2、x1=1.5x_1 = 1.5 から割線法を mpmath の 300 桁で回し、誤差の比から次数を pk=log⁡(ek/ek−1)/log⁡(ek−1/ek−2)p_k = \log(e_k/e_{k-1})/\log(e_{k-1}/e_{k-2}) で推定します。

# Python 3.12.10 / mpmath 1.3.0
from mpmath import mp, mpf, sqrt, log
mp.dps = 300 # 誤差が 1e-250 程度まで追えるように
f = lambda x: x * x - 2
r = sqrt(mpf(2))
xs = [mpf(2), mpf("1.5")] # x0, x1
for k in range(10): # 割線法
a, b = xs[-2], xs[-1]
xs.append(b - f(b) * (b - a) / (f(b) - f(a)))
e = [abs(x - r) for x in xs]
C = 2 / (2 * 2 * r) # f''(x*)/(2 f'(x*)) = 1/(2√2)
print(" k e_k e_k/(e_{k-1}e_{k-2}) 推定次数")
for k in range(2, len(e)):
q = e[k] / (e[k - 1] * e[k - 2])
p = log(e[k] / e[k - 1]) / log(e[k - 1] / e[k - 2])
s = f"{float(e[k]):.2e}".ljust(10)
print(f"{k:2d} {s} {float(q):.5f} {float(p):.4f}")
print(f"1/(2√2) = {float(C):.5f}, 黄金比 = {float((1+sqrt(5))/2):.4f}")
k e_k e_k/(e_{k-1}e_{k-2}) 推定次数
2 1.44e-02 0.28571 0.9305
3 4.21e-04 0.34146 1.9750
4 2.12e-06 0.35172 1.4980
5 3.16e-10 0.35350 1.6666
6 2.37e-16 0.35355 1.6000
7 2.65e-26 0.35355 1.6250
8 2.22e-42 0.35355 1.6154
9 2.08e-68 0.35355 1.6190
10 1.63e-110 0.35355 1.6176
11 1.20e-178 0.35355 1.6182
1/(2√2) = 0.35355, 黄金比 = 1.6180

検算. (a) 極限値 1/(22)1/(2\sqrt2) は Newton 法の例と同じ f′′(x∗)/(2f′(x∗))f''(x^*)/(2f'(x^*)) で、表の第 3 列はそこに収束します(定理 6.2 の 2)。(b) x2−2x^2 - 2 では f[a,b]=a+bf[a,b] = a + b、f[a,b,c]=1f[a,b,c] = 1 なので、証明の式は ek+1=ekek−1/(xk−1+xk)e_{k+1} = e_k e_{k-1}/(x_{k-1} + x_k) と閉じた形になり、第 3 列は 1/(xk−2+xk−1)1/(x_{k-2} + x_{k-1}) に等しくなります(k=2k = 2 で 1/3.5=0.285711/3.5 = 0.28571)。(c) 推定次数は 1.5,1.667,1.6,1.625,1.615,…1.5, 1.667, 1.6, 1.625, 1.615, \ldots と黄金比 1.6181.618 のまわりを振動しながら近づきます。これは誤差の指数が Fibonacci 数列に似た漸化式 log⁡ek+1≈log⁡ek+log⁡ek−1\log e_{k+1} \approx \log e_k + \log e_{k-1} に従うためで、連続する Fibonacci 数の比が黄金比に近づく様子と同じです。

注意 6.6(Brent 法)

実用の 1 変数の求根では、符号の変わる区間を保ちながら、割線法や逆 2 次補間(3 点を通る xx の 2 次式で f=0f = 0 となる点を求める)の候補が区間の中に入り十分に縮むときはそれを採用し、そうでなければ二分法の 1 歩に戻る、という組み合わせ(Dekker 法・Brent 法)が標準です。二分法と同じく収束は保証され、解の近くでは超 1 次収束します。SciPy の brentq はこの算法です。正確な規則と、関数値の評価回数が最悪でも二分法の回数の 2 乗程度で済むことの証明は、R. P. Brent『Algorithms for Minimization without Derivatives』(Prentice-Hall, 1973)第 4 章にあります。

7. 多変数の Newton 法と Kantorovich 型の定理

Section titled “7. 多変数の Newton 法と Kantorovich 型の定理”

U⊂RnU \subset \mathbb{R}^n を開集合、F:U→RnF : U \to \mathbb{R}^n を全微分可能な写像とし、ヤコビ行列を J(x)=JF(x)J(x) = J_F(x) と書きます。Rn\mathbb{R}^n にはユークリッドノルム ∥⋅∥\|\cdot\|、行列には作用素ノルム ∥A∥=sup⁡{∥Ax∥:∥x∥≤1}\|A\| = \sup\{\|Ax\| : \|x\| \le 1\} を使います(多変数関数の微分と同じ記法です)。閉球を Bˉ(a,r)={x:∥x−a∥≤r}\bar B(a, r) = \{x : \|x - a\| \le r\}、開球を B(a,r)B(a, r) と書きます。多変数の Newton 法は

xk+1=xk−J(xk)−1F(xk)x_{k+1} = x_k - J(x_k)^{-1} F(x_k)

です。実際の計算では逆行列を作らず、連立 1 次方程式 J(xk) dk=−F(xk)J(x_k)\, d_k = -F(x_k) を解いて xk+1=xk+dkx_{k+1} = x_k + d_k とします。証明の道具を 2 つ用意します。

補題 7.1(逆行列の摂動)✓

A,BA, B を nn 次の実正方行列とし、AA は正則で ∥A−1(B−A)∥≤q<1\|A^{-1}(B - A)\| \le q \lt 1 とする。このとき BB は正則で

∥B−1∥≤∥A−1∥1−q\|B^{-1}\| \le \frac{\|A^{-1}\|}{1 - q}

である。

証明(補題 7.1)

Bx=0Bx = 0 とすると 0=A−1Bx=x+A−1(B−A)x0 = A^{-1}Bx = x + A^{-1}(B - A)x なので ∥x∥=∥A−1(B−A)x∥≤q∥x∥\|x\| = \|A^{-1}(B - A)x\| \le q\|x\| で、q<1q \lt 1 から x=0x = 0 です。斉次方程式が自明な解しかもたないので、正則性の判定(定理 7.1[行列式とその性質] の 4 ⇒ 2)から BB は正則です。任意の yy について x=B−1yx = B^{-1}y とおくと、上と同じ変形で x=A−1y−A−1(B−A)xx = A^{-1}y - A^{-1}(B - A)x なので

∥x∥≤∥A−1∥ ∥y∥+q∥x∥,∥B−1y∥=∥x∥≤∥A−1∥1−q∥y∥\|x\| \le \|A^{-1}\|\,\|y\| + q\|x\|, \qquad \|B^{-1}y\| = \|x\| \le \frac{\|A^{-1}\|}{1 - q}\|y\|

です。∥y∥≤1\|y\| \le 1 について上限をとれば主張を得ます。

これは Banach 空間でのノイマン級数(定理 6.1[有界作用素と共役作用素])による評価の、有限次元での形です。有限次元では単射と全射が同値なので、級数を使わずに示せます。

次の補題は、多変数の場合に「1 次式で近似した誤差は距離の 2 乗程度」と言うためのものです。証明では、ベクトル値の量を vv の向きの単位ベクトル ww との内積で 1 つの実数にして、1 変数の単調性の議論に持ち込みます。

補題 7.2(1 次近似の剰余の評価)✓

U⊂RnU \subset \mathbb{R}^n を開集合、F:U→RnF : U \to \mathbb{R}^n を全微分可能、a,b∈Ua, b \in U を結ぶ線分 [a,b][a, b] が UU に含まれるとする。γ≥0\gamma \ge 0 について、線分上のすべての点 zz で ∥J(z)−J(a)∥≤γ∥z−a∥\|J(z) - J(a)\| \le \gamma\|z - a\| ならば

∥F(b)−F(a)−J(a)(b−a)∥≤γ2∥b−a∥2\|F(b) - F(a) - J(a)(b - a)\| \le \frac{\gamma}{2}\|b - a\|^2

である。

証明(補題 7.2)

v=F(b)−F(a)−J(a)(b−a)v = F(b) - F(a) - J(a)(b - a) とし、v≠0v \ne 0 としてよいとします。w=v/∥v∥w = v/\|v\| とおきます。V={s∈R:a+s(b−a)∈U}V = \{s \in \mathbb{R} : a + s(b - a) \in U\} は開集合 UU の連続写像による逆像なので開集合で、[0,1][0,1] を含みます。VV 上の実数値関数

ϕ(s)=wT(F(a+s(b−a))−s J(a)(b−a))\phi(s) = w^{\mathsf{T}}\bigl(F(a + s(b - a)) - s\,J(a)(b - a)\bigr)

は、連鎖律(定理 9.4[多変数関数の微分と偏微分])から微分可能で、zs=a+s(b−a)z_s = a + s(b-a) とおくと ϕ′(s)=wT(J(zs)−J(a))(b−a)\phi'(s) = w^{\mathsf{T}}\bigl(J(z_s) - J(a)\bigr)(b - a) です。s∈[0,1]s \in [0,1] では、コーシー・シュワルツの不等式(定理 4.2[内積空間とグラム・シュミット直交化])と ∥w∥=1\|w\| = 1、∥zs−a∥=s∥b−a∥\|z_s - a\| = s\|b - a\| から

∣ϕ′(s)∣≤∥J(zs)−J(a)∥ ∥b−a∥≤γs∥b−a∥2|\phi'(s)| \le \|J(z_s) - J(a)\|\,\|b - a\| \le \gamma s \|b - a\|^2

です。c=γ∥b−a∥2c = \gamma\|b - a\|^2 とし、ψ(s)=ϕ(s)−cs2/2\psi(s) = \phi(s) - c s^2/2 とおくと [0,1][0,1] 上で ψ′(s)=ϕ′(s)−cs≤0\psi'(s) = \phi'(s) - cs \le 0 なので、導関数の符号と単調性(系 3.4[平均値の定理とテイラーの定理]。−ψ-\psi に適用)から ψ(1)≤ψ(0)\psi(1) \le \psi(0)、すなわち ϕ(1)−ϕ(0)≤c/2\phi(1) - \phi(0) \le c/2 です。一方 ϕ(1)−ϕ(0)=wTv=∥v∥\phi(1) - \phi(0) = w^{\mathsf{T}} v = \|v\| なので、∥v∥≤c/2\|v\| \le c/2 です。

平均値不等式(定理 9.7[多変数関数の微分と偏微分])をそのまま使うと係数は γ\gamma になり、1/21/2 が出ません。この 1/21/2 が、Kantorovich 型の定理の条件 h<1/2h \lt 1/2 の 1/21/2 の出どころです。

定理 7.3(多変数の Newton 法の局所 2 次収束)✓

U⊂RnU \subset \mathbb{R}^n を開集合、F:U→RnF : U \to \mathbb{R}^n を全微分可能、x∗∈Ux^* \in U で F(x∗)=0F(x^*) = 0、J(x∗)J(x^*) は正則とし、β=∥J(x∗)−1∥\beta = \|J(x^*)^{-1}\| とおく。r>0r \gt 0 と γ>0\gamma \gt 0 が Bˉ(x∗,r)⊂U\bar B(x^*, r) \subset U と

∥J(x)−J(y)∥≤γ∥x−y∥(x,y∈Bˉ(x∗,r))\|J(x) - J(y)\| \le \gamma\|x - y\| \qquad (x, y \in \bar B(x^*, r))

を満たすとし、ρ=min⁡(r,1/(2βγ))\rho = \min(r, 1/(2\beta\gamma)) とする。このとき ∥x0−x∗∥≤ρ\|x_0 - x^*\| \le \rho なら Newton 法の列は定義されて Bˉ(x∗,ρ)\bar B(x^*, \rho) にとどまり、ek=xk−x∗e_k = x_k - x^* について

∥ek+1∥≤βγ∥ek∥2,βγ∥ek∥≤(βγ∥e0∥)2k≤2−2k\|e_{k+1}\| \le \beta\gamma\|e_k\|^2, \qquad \beta\gamma\|e_k\| \le (\beta\gamma\|e_0\|)^{2^k} \le 2^{-2^k}

が成り立つ。とくに xk→x∗x_k \to x^* である。

証明(定理 7.3)

∥ek∥≤ρ\|e_k\| \le \rho と仮定します。∥J(x∗)−1(J(xk)−J(x∗))∥≤βγ∥ek∥≤1/2\|J(x^*)^{-1}(J(x_k) - J(x^*))\| \le \beta\gamma\|e_k\| \le 1/2 なので、補題 7.1(A=J(x∗)A = J(x^*)、B=J(xk)B = J(x_k)、q=1/2q = 1/2)から J(xk)J(x_k) は正則で ∥J(xk)−1∥≤2β\|J(x_k)^{-1}\| \le 2\beta です。F(x∗)=0F(x^*) = 0 を使うと

ek+1=xk−x∗−J(xk)−1F(xk)=J(xk)−1[F(x∗)−F(xk)−J(xk)(x∗−xk)]e_{k+1} = x_k - x^* - J(x_k)^{-1}F(x_k) = J(x_k)^{-1}\bigl[F(x^*) - F(x_k) - J(x_k)(x^* - x_k)\bigr]

です。閉球は凸なので線分 [xk,x∗][x_k, x^*] は Bˉ(x∗,r)\bar B(x^*, r) に含まれ、その上で ∥J(z)−J(xk)∥≤γ∥z−xk∥\|J(z) - J(x_k)\| \le \gamma\|z - x_k\| です。補題 7.2(a=xka = x_k、b=x∗b = x^*)から括弧の中のノルムは (γ/2)∥ek∥2(\gamma/2)\|e_k\|^2 以下なので、

∥ek+1∥≤2β⋅γ2∥ek∥2=βγ∥ek∥2≤βγρ ∥ek∥≤12∥ek∥\|e_{k+1}\| \le 2\beta \cdot \frac{\gamma}{2}\|e_k\|^2 = \beta\gamma\|e_k\|^2 \le \beta\gamma\rho\,\|e_k\| \le \frac12\|e_k\|

です。よって xk+1∈Bˉ(x∗,ρ)x_{k+1} \in \bar B(x^*, \rho) で、帰納法ですべての kk で成り立ちます。βγ∥ek+1∥≤(βγ∥ek∥)2\beta\gamma\|e_{k+1}\| \le (\beta\gamma\|e_k\|)^2 を繰り返せば 2 つ目の不等式を得ます。

この定理も解 x∗x^* を中心に述べているので、解が存在することは前提です。次の定理は、解を知らないまま、初期値 x0x_0 で計算できる 3 つの量 β,η,γ\beta, \eta, \gamma から、解の存在と一意性と、解を含む球の半径を結論します。

定理 7.4(Kantorovich 型の存在と一意性(簡易 Newton 法))✓

U⊂RnU \subset \mathbb{R}^n を開集合、F:U→RnF : U \to \mathbb{R}^n を全微分可能、x0∈Ux_0 \in U で J0=J(x0)J_0 = J(x_0) は正則とする。正の数 β,η,γ,R\beta, \eta, \gamma, R が

∥J0−1∥≤β,∥J0−1F(x0)∥≤η,∥J(x)−J0∥≤γ∥x−x0∥(x∈Bˉ(x0,R))\|J_0^{-1}\| \le \beta, \qquad \|J_0^{-1}F(x_0)\| \le \eta, \qquad \|J(x) - J_0\| \le \gamma\|x - x_0\| \quad (x \in \bar B(x_0, R))

と Bˉ(x0,R)⊂U\bar B(x_0, R) \subset U を満たすとし、L=βγL = \beta\gamma、h=Lηh = L\eta とおく。h<1/2h \lt 1/2 かつ

t−=1−1−2hL≤Rt_- = \frac{1 - \sqrt{1 - 2h}}{L} \le R

ならば次が成り立つ。

  1. FF は Bˉ(x0,t−)\bar B(x_0, t_-) に零点 x∗x^* をもち、Bˉ(x0,t−)\bar B(x_0, t_-) の中の零点はそれだけである。
  2. 簡易 Newton 法 xk+1=xk−J0−1F(xk)x_{k+1} = x_k - J_0^{-1}F(x_k) の列は Bˉ(x0,t−)\bar B(x_0, t_-) にとどまって x∗x^* に収束し、q=Lt−=1−1−2h<1q = L t_- = 1 - \sqrt{1 - 2h} \lt 1 について ∥xk−x∗∥≤qkη/(1−q)\|x_k - x^*\| \le q^k\eta/(1 - q) である。
  3. FF は開球 B(x0,min⁡(R,1/L))B(x_0, \min(R, 1/L)) の中に x∗x^* 以外の零点をもたない。
証明(定理 7.4)

T(x)=x−J0−1F(x)T(x) = x - J_0^{-1}F(x) とおきます。J0−1J_0^{-1} は正則なので、T(x)=xT(x) = x と F(x)=0F(x) = 0 は同値です。

縮小性. 0<t≤R0 \lt t \le R とします。TT は全微分可能で JT(z)=I−J0−1J(z)=J0−1(J0−J(z))J_T(z) = I - J_0^{-1}J(z) = J_0^{-1}(J_0 - J(z)) なので、z∈Bˉ(x0,t)z \in \bar B(x_0, t) で ∥JT(z)∥≤βγ∥z−x0∥≤Lt\|J_T(z)\| \le \beta\gamma\|z - x_0\| \le Lt です。x,y∈Bˉ(x0,t)x, y \in \bar B(x_0, t) を結ぶ線分は閉球に含まれるので、平均値不等式(定理 9.7[多変数関数の微分と偏微分])から

∥T(x)−T(y)∥≤Lt ∥x−y∥\|T(x) - T(y)\| \le Lt\,\|x - y\|

です。

閉球を自分自身に写すこと. x∈Bˉ(x0,t)x \in \bar B(x_0, t) について

T(x)−x0=−J0−1[F(x)−F(x0)−J0(x−x0)]−J0−1F(x0)T(x) - x_0 = -J_0^{-1}\bigl[F(x) - F(x_0) - J_0(x - x_0)\bigr] - J_0^{-1}F(x_0)

です。仮定から線分 [x0,x][x_0, x] 上で ∥J(z)−J0∥≤γ∥z−x0∥\|J(z) - J_0\| \le \gamma\|z - x_0\| なので、補題 7.2(a=x0a = x_0、b=xb = x)から

∥T(x)−x0∥≤β⋅γ2∥x−x0∥2+η≤L2t2+η\|T(x) - x_0\| \le \beta\cdot\frac{\gamma}{2}\|x - x_0\|^2 + \eta \le \frac{L}{2}t^2 + \eta

です。t−t_- は 2 次方程式 L2t2−t+η=0\frac{L}{2}t^2 - t + \eta = 0 の小さいほうの根(h<1/2h \lt 1/2 なので実根)なので、t=t−t = t_- のとき右辺はちょうど t−t_- です。よって TT は Bˉ(x0,t−)\bar B(x_0, t_-) を自分自身に写します。

1・2. Bˉ(x0,t−)\bar B(x_0, t_-) は Rn\mathbb{R}^n の空でない閉集合で、TT はそこで縮小率 q=Lt−=1−1−2hq = Lt_- = 1 - \sqrt{1 - 2h} の縮小写像です。h<1/2h \lt 1/2 から q<1q \lt 1 です。縮小写像の原理(定理 3.2[逆関数定理・陰関数定理と条件付き極値])から TT は Bˉ(x0,t−)\bar B(x_0, t_-) にただ 1 つの不動点 x∗x^* をもち、簡易 Newton 法の列(TT の反復そのもの)はそこに収束して ∥xk−x∗∥≤qk∥x1−x0∥/(1−q)\|x_k - x^*\| \le q^k\|x_1 - x_0\|/(1-q) です。∥x1−x0∥=∥J0−1F(x0)∥≤η\|x_1 - x_0\| = \|J_0^{-1}F(x_0)\| \le \eta から 2 の評価を得ます。不動点の一意性は、この閉球の中の零点の一意性です。

3. yy を ∥y−x0∥<min⁡(R,1/L)\|y - x_0\| \lt \min(R, 1/L) を満たす零点とし、t=max⁡(∥y−x0∥,t−)t = \max(\|y - x_0\|, t_-) とおきます。t−≤Rt_- \le R と Lt−=q<1Lt_- = q \lt 1 から t≤Rt \le R、Lt<1Lt \lt 1 です。x∗x^* と yy はともに Bˉ(x0,t)\bar B(x_0, t) の中の TT の不動点なので、縮小性から ∥x∗−y∥=∥T(x∗)−T(y)∥≤Lt∥x∗−y∥\|x^* - y\| = \|T(x^*) - T(y)\| \le Lt\|x^* - y\| となり、x∗=yx^* = y です。

条件 h=βγη<1/2h = \beta\gamma\eta \lt 1/2 は、「初期値での Newton 補正の大きさ η\eta」が「ヤコビ行列の変化の速さ γ\gamma」と「ヤコビ行列の逆の大きさ β\beta」に比べて十分小さい、という意味です。左辺の 3 つの量はすべて x0x_0 と FF の情報から計算でき、とくに区間演算で上界を計算すれば、計算機の出力がそのまま「この球の中に解がただ 1 つある」ことの証明になります。

x₀x*(零点)t₋破線の球の中の零点は x* だけ
Kantorovich 型の定理の状況(模式図)。初期値 x₀ を中心とする半径 t₋ の閉球の中に零点 x* があり、より大きい開球(半径 min(R, 1/L))の中では零点は x* だけである。

注意 7.5(Kantorovich の定理の完全な形)

本章で証明したのは簡易 Newton 法(ヤコビ行列を J0J_0 に固定)が縮小写像になる形です。L. V. Kantorovich による元の定理は、同じ β,η\beta, \eta と、x0x_0 を中心とする条件より強い球の中での Lipschitz 条件 ∥J(x)−J(y)∥≤γ∥x−y∥\|J(x) - J(y)\| \le \gamma\|x - y\|(x,y∈Bˉ(x0,R)x, y \in \bar B(x_0, R))を満たす γ\gamma について、h≤1/2h \le 1/2 のもとで、Newton 法そのものの列が定義されて Bˉ(x0,t−)\bar B(x_0, t_-) の中の零点に収束し、h<1/2h \lt 1/2 なら 2 次収束することを述べます。証明は Newton 法の列を 1 変数の多項式 L2t2−t+η\frac{L}{2}t^2 - t + \eta に Newton 法を使った列で上から抑える「優級数」の議論で、不等式の帰納法が 2〜3 頁続きます。J. M. Ortega, W. C. Rheinboldt『Iterative Solution of Nonlinear Equations in Several Variables』(Academic Press, 1970)§12.6、P. Deuflhard『Newton Methods for Nonlinear Problems』(Springer, 2004)第 2 章を参照してください。ここでは定理環境にせず、後の章の証明でも使いません。

演習 8.1易

ff が [1,2][1, 2] で連続、f(1)f(2)<0f(1) f(2) \lt 0 とします。二分法の中点 ckc_k の誤差が 10−610^{-6} 以下であることが定理 3.1から保証される最小の kk を求めてください。

解答

誤差の上界は (2−1)/2k+1=2−(k+1)(2 - 1)/2^{k+1} = 2^{-(k+1)} です。2−(k+1)≤10−62^{-(k+1)} \le 10^{-6} は k+1≥6log⁡210≈19.93k + 1 \ge 6\log_2 10 \approx 19.93、すなわち k≥19k \ge 19 と同値です。よって k=19k = 19 で、c0,…,c19c_0, \ldots, c_{19} の 20 回の関数値の計算で保証されます。

演習 8.2標準

a>0a \gt 0 とします。f(x)=1/x−af(x) = 1/x - a に Newton 法を使うと xk+1=xk(2−axk)x_{k+1} = x_k(2 - a x_k) となり、割り算なしで 1/a1/a を計算できることを確かめてください。さらに rk=1−axkr_k = 1 - a x_k とおいて rk+1=rk2r_{k+1} = r_k^2 を示し、列が 1/a1/a に収束するための x0x_0 の必要十分条件と、ek+1/ek2e_{k+1}/e_k^2 の値を求めてください。

解答

f′(x)=−1/x2f'(x) = -1/x^2 なので xk+1=xk+xk2(1/xk−a)=xk(2−axk)x_{k+1} = x_k + x_k^2(1/x_k - a) = x_k(2 - a x_k) です。rk+1=1−axk(2−axk)=1−2axk+a2xk2=(1−axk)2=rk2r_{k+1} = 1 - a x_k(2 - a x_k) = 1 - 2a x_k + a^2 x_k^2 = (1 - a x_k)^2 = r_k^2 なので rk=r02kr_k = r_0^{2^k} です。xk→1/ax_k \to 1/a は rk→0r_k \to 0 と同値で、それは ∣r0∣<1|r_0| \lt 1、すなわち 0<x0<2/a0 \lt x_0 \lt 2/a と同値です(∣r0∣=1|r_0| = 1 なら rk=1r_k = 1 で xk=0x_k = 0 にとどまり、∣r0∣>1|r_0| \gt 1 なら ∣rk∣→∞|r_k| \to \infty です)。ek=xk−1/a=−rk/ae_k = x_k - 1/a = -r_k/a なので ek+1=−rk2/a=−aek2e_{k+1} = -r_k^2/a = -a e_k^2、すなわち ek+1/ek2=−ae_{k+1}/e_k^2 = -a です。これは 定理 5.1 の極限 f′′(x∗)/(2f′(x∗))=(2a3)/(2⋅(−a2))=−af''(x^*)/(2f'(x^*)) = (2a^3)/(2 \cdot (-a^2)) = -a と一致します。

演習 8.3標準

a>0a \gt 0 とし、g(x)=x(x2+3a)/(3x2+a)g(x) = x(x^2 + 3a)/(3x^2 + a) とします。g(a)=ag(\sqrt a) = \sqrt a と

g(x)−a=(x−a)33x2+ag(x) - \sqrt a = \frac{(x - \sqrt a)^3}{3x^2 + a}

を示し、x0x_0 が a\sqrt a に十分近いとき不動点反復が 3 次収束することと、漸近誤差定数を求めてください。

解答

s=as = \sqrt a とおくと (x−s)3=x3−3sx2+3ax−as(x - s)^3 = x^3 - 3sx^2 + 3ax - as です。一方 g(x)−s=(x3+3ax−s(3x2+a))/(3x2+a)=(x3−3sx2+3ax−as)/(3x2+a)g(x) - s = \bigl(x^3 + 3ax - s(3x^2 + a)\bigr)/(3x^2 + a) = (x^3 - 3sx^2 + 3ax - as)/(3x^2 + a) なので等式が成り立ち、とくに g(s)=sg(s) = s です。ek≠0e_k \ne 0 なら ek+1/ek3=1/(3xk2+a)e_{k+1}/e_k^3 = 1/(3x_k^2 + a) です。gg は x>0x \gt 0 で C1C^1 級で、等式から g′(s)=0g'(s) = 0 なので、定理 4.3 の 1 により ss の近くから始めた列は ss に収束します。よって ek+1/ek3→1/(3a+a)=1/(4a)e_{k+1}/e_k^3 \to 1/(3a + a) = 1/(4a) で、3 次収束、漸近誤差定数は 1/(4a)1/(4a) です(定理 4.3 の 2 で p=3p = 3 としても同じ結論で、g′′′(s)/3!=1/(4a)g'''(s)/3! = 1/(4a) です)。

演習 8.4難

F(x,y)=(x2+y2−1, x−y)F(x, y) = (x^2 + y^2 - 1,\ x - y) とします。(i) (1,0)(1, 0) から多変数の Newton 法を 2 回行い、(3/4,3/4)(3/4, 3/4) を得ることを確かめてください。(ii) x0=(3/4,3/4)x_0 = (3/4, 3/4) で定理 7.4の β,η,γ\beta, \eta, \gamma(RR は任意)を求め、hh と t−t_- を計算して、零点 (1/2,1/2)(1/\sqrt2, 1/\sqrt2) が Bˉ(x0,t−)\bar B(x_0, t_-) に入っていることを確かめてください。

解答

(i) J(x,y)=(2x2y1−1)J(x, y) = \begin{pmatrix} 2x & 2y \\ 1 & -1 \end{pmatrix} です。(1,0)(1,0) では J=(201−1)J = \begin{pmatrix} 2 & 0 \\ 1 & -1\end{pmatrix}、F=(0,1)F = (0, 1) で、Jd=−FJ d = -F を解くと d=(0,1)d = (0, 1)、次の点は (1,1)(1, 1) です。そこでは J=(221−1)J = \begin{pmatrix} 2 & 2 \\ 1 & -1\end{pmatrix}、F=(1,0)F = (1, 0) で、d=(−1/4,−1/4)d = (-1/4, -1/4)、次の点は (3/4,3/4)(3/4, 3/4) です。

(ii) J(p)−J(p′)=(2(p1−p1′)2(p2−p2′)00)J(p) - J(p') = \begin{pmatrix} 2(p_1 - p'_1) & 2(p_2 - p'_2) \\ 0 & 0\end{pmatrix} の作用素ノルムは第 1 行のベクトルのノルム 2∥p−p′∥2\|p - p'\| なので、γ=2\gamma = 2(どこでも成り立つので RR は任意)です。J0=(3/23/21−1)J_0 = \begin{pmatrix} 3/2 & 3/2 \\ 1 & -1\end{pmatrix}、det⁡J0=−3\det J_0 = -3 で、J0−1=(1/31/21/3−1/2)J_0^{-1} = \begin{pmatrix} 1/3 & 1/2 \\ 1/3 & -1/2\end{pmatrix} です。この行列の 2 つの列は直交し、ノルムは 2/3\sqrt2/3 と 2/2\sqrt2/2 なので、∥J0−1v∥2=v12⋅2/9+v22⋅1/2\|J_0^{-1}v\|^2 = v_1^2 \cdot 2/9 + v_2^2 \cdot 1/2 から β=∥J0−1∥=2/2\beta = \|J_0^{-1}\| = \sqrt2/2 です。F(x0)=(1/8,0)F(x_0) = (1/8, 0)、J0−1F(x0)=(1/24,1/24)J_0^{-1}F(x_0) = (1/24, 1/24) なので η=2/24\eta = \sqrt2/24 です。L=βγ=2L = \beta\gamma = \sqrt2、h=Lη=1/12<1/2h = L\eta = 1/12 \lt 1/2 で、

t−=1−5/62≈0.06161t_- = \frac{1 - \sqrt{5/6}}{\sqrt2} \approx 0.06161

です。零点との距離は 2 (3/4−1/2)≈0.06066<t−\sqrt2\,(3/4 - 1/\sqrt2) \approx 0.06066 \lt t_- で、たしかに球に入っています。さらに 3 から、半径 1/L≈0.7071/L \approx 0.707 の開球には零点がこれ以外にありません(もう 1 つの零点 (−1/2,−1/2)(-1/\sqrt2, -1/\sqrt2) との距離は約 2.062.06 です)。

  • 山本哲朗『数値解析入門』(増訂版)サイエンス社、2003 — 非線形方程式の章(二分法・Newton 法・Kantorovich の定理)。
  • 杉原正顯・室田一雄『数値計算法の数理』岩波書店、1994 — 非線形方程式の章。
  • K. E. Atkinson, An Introduction to Numerical Analysis, 2nd ed., Wiley, 1989 — Chapter 2(求根法、割線法の次数)。
  • J. M. Ortega, W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, 1970(SIAM Classics in Applied Mathematics 30, 2000)— 収束の次数と Newton–Kantorovich の定理。
  • P. Deuflhard, Newton Methods for Nonlinear Problems: Affine Invariance and Adaptive Algorithms, Springer, 2004 — Chapter 2。
  • R. P. Brent, Algorithms for Minimization without Derivatives, Prentice-Hall, 1973 — Chapter 4(Brent 法)。

このテーマの 100 本ノック(50 問)

執筆: Claude Code(Opus 5.5)、Claude Code(Sonnet 5.5)・監査 2 回(Claude Code)・最終 2026-10-12

この記事の誤りを報告する ・運営: 夢現技研合同会社 ・料金 ・利用条件 ・特定商取引法に基づく表記

© 2026 夢現技研合同会社 ・コード例は MIT ライセンスです。