お話：数値解析 第 10 回
非線型方程式 (前編)長田直樹

1

4. q = −0.0054+r とおき、f (2) (q) の q 3 を無視した
f¯(2) (q) = 6.3q 2 + 11.23q + 0.061 = 0 に代入し、
f (3) (r) = 6.3r2 + 11.16196r + 0.000541708 =はじめに
非線型方程式は線型方程式 (連立 1 次方程式) に対する概念で、単独方程式0 を得る。f (3) (r) = 0 の 6.3r2 を無視して
r = −0.000541708/11.16196 = −0.00004853を得る。

f (x) = 0と連立方程式

5. y = 2+0.1−0.0054−0.00004853 = 2.09455147
が解の近似値である。

 f1 (x1 , . . . , xn ) = 0···
fn (x1 , . . . , xn ) = 0

現代の記号を用いて説明すると次のようになる。pは初期値 2 の補正であるので、初期値が解に十分近いときには p は十分小さい。f (1) (p) = f (2 + p) とおがある。

非線型方程式の数値解は反復解法で求める。最も くと f (y) は 3 次多項式なのでテイラー展開は有名な反復解法がニュートン法 (ニュートン・ラフソ

f (1) (p) = f (2) + pf 0 (2) +ン法とも呼ばれている) である。今回は非線型単独

p2 00
p3f (2) + f 000 (2)26となる。p3 と p2 の項を無視して f (1) (p) = 0 を解

方程式に対するニュートン法の話をする。

くと2p=−

ニュートンとラフソン

f (2),
f 0 (2)y =2−f (2)f 0 (2)となる。次のステップでは f (1) (p) に対し初期値を万有引力の法則などで名高い I. ニュートン [5,pp.268-269] は 1669 年、3 次方程式

0.1、補正を q として繰り返しているので
q=−f (y) = f (0) (y) = y 3 − 2y − 5 = 0f (1) (0.1)f(1) 0(0.1),

p = 0.1 −f (1) (0.1)
0f (1) (0.1)

である。の実数解を次のようにして求めた [1]。

王立協会会員の J. ラフソンは 1690 年小冊子を出版した。その中で 3 次方程式 f (a) = a3 − ba − c = 0

1. 2 を初期値に取る。

に対し、解の推定値 g に対するよりよい推定値

2. y = 2 + p を方程式 f (0) (y) = 0 に代入し、p の方程式 f (1) (p) = p3 + 6p2 + 10p − 1 = 0 を得g+る。高次の項を無視した 10p − 1 = 0 を解いてc + bg − g 3f (g)=g− 0
3g 2 − bf (g)を与えた [11]。

p = 0.1 を得る。

ニュートンの方法はステップ毎に代入する関数が(1)3. p = 0.1 + q を方程式 f (p) = 0 に代入し、 異なるのに対し、ラフソンのは毎回同じ関数に代入
f (2) (q) = q 3 + 6.3q 2 + 11.23q + 0.061 = 0 を得 するという違いがある。したがって、反復る。f (2) (q) = 0 の高次の項を無視した 11.23q +f (x(ν) )x(ν+1) = x(ν) − 0 (ν)
0.061 = 0 を解いて q = −0.0054 を得る。
f (x )
1

4. 定数項【実】−0.044424 は 0 にならないので、【実】の符号を変えた 0.044424 を 1 次の係はラフソン法あるいはニュートン・ラフソン法と呼ぶのが適切であるが、本連載では習慣に従いニュー数【方】12.8028 で割り、0.044424/12.8028 =0.00346 強を求める。

トン法と呼ぶ。
ニュートンもラフソンも微分及び組み立て除法は用いていないことを注意しておく。

5. x = 1 + 0.2 + 0.06 + 0.00346 強 = 1.26346 強を解【定商】とする。

3関孝和

1-3 は組み立て除法により解を 1 桁ずつ求めてい関孝和は 1685 年に改訂または執筆した『解隠題之る。4-5 は (1) の高次の項を無視して r を求めたこと法』[3] (かいいんだいのほう) のなかで代数方程式のになり、ニュートンの方法と同一である。これによ数値解法を扱っている。関が最後に取り上げた 3 次り、解の有効桁数を 3 桁から 6 桁と倍増させている。

方程式

1683 年に関孝和と弟子の建部賢明 (たけべかたあきら)、賢弘 (かたひろ) 兄弟の 3 人により算法の集大f (x) = −9 + 3x + 2x2 + x3 = 0

成が始められた。1695 年頃に全 12 巻の『算法大成』を用いて、前節と対応させながら解法を説明する。 が一応完成したが、その後も編纂が続けられ関の没関は多項式を縦書きの昇べき順で表したので、昇べ後 1711 年頃建部賢明の手により『大成算経』[7](たき順を用いている。
【】内は『解隠題之法』の用語でいせいさんけい) 全 20 巻が完成した。
『大成算経』巻三は代数方程式の数値解法を扱ってある。

いる。その中で述べられている「窮商」(きゅうしょう)

1. 【商】1 を立て、組み立て除法を行う。
1)−9633211−9 + 3x + 2x2 + x3−3631右側から計算する

410141とは、代数方程式の数値解を高精度で求める方法である。窮商の一般的説明に引き続き 2 次方程式と『解隠題之法』と同一の 3 次方程式 −9 + 3x + 2x2 + x3 = 0の例が取り上げられている。解法を説明している原文を図 1 に示す。

151

前節の記法に合わせると、x = 1 + p を方程式f (x) = 0 に代入し、p の方程式 −3+10p+5p2 +p3 = 0 を得る。
2. 【商】0.2 を立て、組み立て除法を行う。
0.2) −3102.208 1.04

50.2

1−0.792 11.04
1.08

5.20.2112.12

5.4

10.2
5.6

1−0.792 + 12.12q + 5.6q 2 + q 3 = 0 が得られる。
3. 【商】0.06 を立て、組み立て除法を行うと

図 1 :『大成算経』巻三 (東北大学)−0.044424 + 12.8028r + 5.78r2 + r3 = 0 (1)図 1 は小さくて見づらいので、が得られる。

2

東北大学和算ポータルとなる。ラフソンの方法と同一で、今日のニュート(http://www2.library.tohoku.ac.jp/wasan/)ン法そのものである。
→ 江戸後期刊本「窮商」で用いられたニュートン法は、関と建部→ 大成算経 (狩野 7.0820.20) → 三冊兄弟の三人のうち誰がいつ発見したか明らかになっていない。関が『解隠題之法』で示した補正についをたどり、番号 123 を開き拡大画像で見るとよい。てのニュートン法を建部賢明が発展させ整理したの
図 1 の 3 行目から 7 行目までは『解隠題之法』のではないかと想像している。
上記 5 までの要約である。8 行目から 12 行目までは次のようにしている。以下【】内は『大成算経』の用語である。

ニュートン法の収束

4

6. 1 + 0.2 + 0.06 + 0.00346 = 1.26346 4.1 縮小写像の原理を【 次 商 】と し 、組 み 立 て 除 法 を 行 う。
集合 D ⊂ R で定義された関数 g(x) に対し α =(誌 面 の 都 合 で 小 数 第 9 位 ま で 表 す。)
1.26346) −98.999942925−0.000057074

34.123251171
7.123251171
5.719582343
12.842833514

2
1.26346
3.26346
1.26346
4.52692
1.263465.79038

1 g(α) となる点 α ∈ D を g(x) の不動点という。
1ある定数 L(0 < L < 1) に対して|g(x1 ) − g(x2 )| 5 L|x1 − x2 | (x1 , x2 ∈ D)

11を満たすとき g(x) は縮小写像と呼ばれる。

ニュートン法などの反復が収束するための十分条= d + 1.26346 件を与える定理に縮小写像の原理がある。
を元の方程式 f (x) = 0 に代入し補題 1 (縮小写像の原理) D を R の閉集合とする。
縮小写像 g : D → D は D にただ一つの不動点を持−0.000057074730264 + 12.8428335148d
前節の記法を用いると xつ。任意の初期値 x(0) ∈ D に対し、反復+5.79038d2 + d3 = 0x(ν+1) = g(x(ν) )を得ている。【実】の符号を変えたものを【方】

は α に収束する。

で割り

d=0.000057074730264= 0.00000444409
12.8428335148

証明 数値解析あるいは関数解析のテキスト (たとえば、[8, pp.64-65][10, pp.69-70]) を見よ。

を得る。

¤

補題 1 は R における縮小写像の原理であるが、Rn 、

7. 1.26346 + 0.00000444409 = 1.26346444409 が Cn 、あるいはバナッハ空間 (完備ノルム空間) においても同様の定理が成り立つ。n > 1 の場合は、絶次の近似値【三商】である。
対値の代りにノルムが用いられる。

8. 繰り返すと詳しい値が求められる。
現代的に表すと4.2

1. x(1) = 1.26

α が f (x) の単純零点とは、f (α) = 0, f 0 (α) 6= 0のときいう。

2. 【次商】
x(2) = x(1) −定理 1 f (x) が単純零点 α を含む区間 I˜ で C 2 級で、
f 0 (x) は I˜ に零点を持たないものとする。このとき、

f (x(1) )= 1.26346
f 0 (x(1) )閉区間 I ⊂ I˜ と正の数 M > 0 が存在して、反復

3. 【三商】
x(3) = x(2) −単純零点x(0) ∈ I(2)f (x )= 1.26346444409
f 0 (x(2) )x(ν+1) = x(ν) −3f (x(ν) )f 0 (x(ν) )(2)に対し、1 回の反復毎に有効桁数が 2 倍になっており、4回の反復で倍精度での限界の値が得られる。真値は|x(ν+1) − α| 5 M |x(ν) − α|22.09455148154232659 である。本例では (3) の δ は1.78 × 10−15 より大きく取る必要がある。

が成立する。
証明 まず、g(x) = x − f (x)/f 0 (x) が補題 1 の条件

f (x) を n 次多項式とする。解の絶対値は大きくなを満たすことを示す。

く、初期値の絶対値が大きいときは減少率 1 − 1/nf (x)f 00 (x)g (x) =(f 0 (x))20で 0 近づく。

1limx→∞ xより g 0 (α) = 0 である。g 0 (x) は連続だから、0 < L <1 に対し、δ > 0 が存在して I = [α − δ, α + δ] ⊂ I˜µ¶f (x)1x− 0=1−f (x)nにおいて、|g 0 (x)| < L となる。α は g(x) の不動点でからいえる。

ありので、補題 1 より反復 (2) は α に収束する。

例 2 例 1 の関数 f (x) = x3 − 2x − 5 に初期値 x(0) =次に、x ∈ I をとる。テイラーの定理により
1f (α) = f (x) + f 0 (x)(α − x) + f 00 (ξx )(α − x)22となる ξx ∈ I が存在する。f (α) = 0 より

1000 を用いてニュートン法を適用する。
νx(ν)|f (x(ν) )|01000.000000000000000

1.00 × 109

12666.667112778075193
444.445412269271287

2.96 × 108
8.78 × 107

1718

中略
2.104403822434531
2.094605697648467

1.11 × 10−1
6.05 × 10−4

1920

2.094551483197066
2.094551481542327

1.85 × 10−8
1.78 × 10−15

f 00 (ξx )f (x)=(x − α)2x−α− 0f (x)2f 0 (x)となる。仮定より、A 5 |f 0 (x)|, |f 00 (x)| 5 B(x ∈ I)となる定数 A, B > 0 が存在する。M = B/(2A), x =x(ν) とおくと|x(ν+1) − α| 5 M |x(ν) − α|2¤が得られる。

ν = 1, 2 のときは減少率 2/3 で 0 に近づき、ν =17, . . . , 20 のときは唯一の実数解に 2 次収束してい定理 1 より、初期値を解の近くに取ればニュートる。

ン法は 2 次収束する。倍精度計算のとき通常 4 ∼ 5回の反復で解が得られる。
反復は f (x) の計算の際の丸め誤差の上界 (のでき4.3るだけ小さい値) を δ とするとき、|f (x(ν) )| < δ

多重零点

α が f (x) の m 重零点とは、α で有界な関数 h(x)が存在し、(3)となったときに終了すればよい [4]。δ を小さく取り

f (x) = (x − α)m h(x),すぎると反復は終了しないことがある。このようなh(α) 6= 0事態を避けるため、プログラミングに際して反復回となるときいう。f (x) が α で C m 級のときは、数の上限を設けておく。

f (α) = f 0 (α) = · · · = f (m−1) (α) = 0, f (m) (α) 6= 0

例 1 ニュートンが扱った関数 f (x) = x3 − 2x − 5 とと同値になる。(証明は容易である。) 2 重零点、3 重初期値 x(0) = 2 を用いてニュートン法を適用する。

零点、· · · をまとめて多重零点という。

|f (x(ν) )|νx(ν)

02.000000000000000

1.00

定理 2 f (x) が m 重零点 α を含む開区間 I において122.100000000000000
2.094568121104185

6.10 × 10−2
1.86 × 10−4C m 級で、f 0 (x) 6= 0, (x ∈ I \ {α}) とする。x(0) ∈2.094551481698199
2.094551481542327−934I \ {α} をとると、反復

1.74 × 10
1.78 × 10−15x(ν+1) = x(ν) −4

f (x(ν) )f 0 (x(ν) )はとなる。m 重零点の有効数字の桁数は、計算桁数の1x(ν+1) − αlim
=1−ν→∞ x(ν) − α
m1/m 程度であることが分かる。たとえば、2 重解をδ = 10−15 で求めたとき、|h(x(ν) )| ; 1 とすると予√想される誤差は 10−15 = 3.2 × 10−8 となる。このを満たす。

ことは、ニュートン法にもシュレーダー法にも当て

証明 定義よりはまる。
0f (x) = m(x − α)m−1m 0h(x) + (x − α) h (x)例 3 関数となる。

f (x) = x3 − 5x2 + 8x − 4f (x)(x − α)m h(x)
=(4)
0
m−1
f (x)
m(x − α)h(x) + (x − α)m h0 (x)= (x − 1)(x − 2)2に初期値 x(0) = 3 を用いてニュートン法と m = 2 とより

f (x)x− 0−α=x−α−f (x)

してシュレーダー法を適用する。|f (x(ν) )| < 10−15x−αm + (x − α)h0 (x)
h(x)となったとき反復を停止する。

f (x)
−α
f 0 (x)=1−x−α1m + (x − α)0h (x)h(x)x = x(ν) として ν → ∞ により得られる。

¤E. シュレーダー [6] が 1870 年に提案したシュレーx(ν)

x(ν+1) = x(ν) − mf (x(ν) )f 0 (x(ν) )(5)は m 重零点に 2 次収束する。

は

03.0

3.0

1
23

2.60
2.35
2.19

2.20
2.015
2.0001242.10
中略2.0000000067

で反復が終了する。ニュートン法は縮小率 1/2 の線定理 3 定理 2 と同じ条件の下で、反復x(ν+1) = x(ν) − mx

24 2.00000011
25 2.000000056
26 2.0000000024
ニュートン法では f (x(26) ) = 8.88 × 10−16 で反復が終了し、シュレーダー法は f (x(4) ) = 8.88 × 10−16ダー法

型収束、シュレーダー法は 2 次収束であるが、得らf (x(ν) )f 0 (x(ν) )れた解の誤差はともに 10−8 程度である。

x(ν+1) − α
1 h0 (α)=(ν)
2ν→∞ (xm h(α)
− α)

5lim

正則関数に対するニュートン法
ここまで実数解のみを扱ってきたが、複素数計算を満たす。

ができるシステムでは複素数解を扱うこともできる。

証明 (4) を用いると定理 2 と同様に証明できる。 ¤必要最小限の事項をまとめておく。

初期値 x(0) を f (x) = (x − α)m h(x) の零点 α の近くに取り、|f (x(ν)

複素平面上の集合 D の任意の 2 点が D 内の折れ線で結べるとき連結という。連結開集合を領域、連)| < δ で反復が終了したとする。

結閉集合を閉領域という。領域 D の上で定義された|(x(ν) − α)m h(x(ν) )| < δより、シュレーダー法

(ν)

ν

よって、

x−ニュートン法

関数 f (z) が z ∈ D で微分可能であるとは、複素数

h が 0 以外の値を取りながら 0 に近づくとき¯¯¯ δ ¯1/m¯|x(ν) − α| < ¯¯h(x(ν) ) ¯

f 0 (z) = lim

h→05f (z + h) − f (z)h(6)

が存在して有限のときいう。(6) は実関数の微分と同6じ形をしているが h は 0 の任意の近傍の点を取り得

ニュートン法の周期点複素平面上の写像 g : C → C の反復合成をるので、かなり強い条件になっている。f (z) が D の各点で微分可能なとき正則という。

g 0 (z) = z

補題 2 f (z) が開円盤領域 D = {|z − a| < R} で正g ν+1 (z) = g(g ν (z)),

則のとき、f (z) は級数f (z) =により定義する。ある自然数 p に対し、∞Xf (k) (a)k=0

ν = 0, 1, . . .

k!

(z − a)kg(α) 6= α, · · · , g p−1 (α) 6= α, g p (α) = α(7)となるとき、α ∈ C を g の周期pの周期点という。に展開できる。級数 (7) はすべての z ∈ D に対し収多項式 f (z) に対するニュートン法の反復関数束する。

g(z) = z −¤級数 (7) を a を中心とするテイラー展開という。

証明 関数論のテキストを見よ。

f (z)
f 0 (z)により周期点に収束する初期値の集合は、空でなければたいていの場合フラクタルと呼ばれる複雑な形定理 4 f (z) が単純零点 α を含む領域 D̃ で正則で、 をする。フラクタルとはおおざっぱにいうと、図形
f 0 (z) は D̃ に零点を持たないものとする。閉領域 D ⊂ のどんな一部分も全体と同じ複雑さを持つという広D̃ と正の数 M > 0 が存在して、初期値 z (0) ∈ D をとった反復z(ν+1)=z(ν)

い意味での自己相似性を持つ集合である。ニュートン法の反復から生じるフラクタルについては複素力学系のテキスト (たとえば [2, 9]) を見よ。

f (z (ν) )− 0 (ν)f (z )

f (z) = z 3 − 2z + 2 に対するニュートン法z (ν+1) = z (ν) −に対し、(z (ν) )3 − 2z (ν) + 23(z (ν) )2 − 2(8)を考える。z (ν) = 0 または 1 となると、ν から先は 0|z (ν+1) − α| 5 M |z (ν) − α|2と 1 を繰り返し収束しない。0 と 1 は反復 (8) に対すが成立する。

る周期 2 の周期点である。周期点 0 と 1 に収束する
複素平面上 −100 5 Re z 5 100, −80 5 Im z 5 80¤

証明 定理 1 とほぼ同様である。の範囲の初期値の集合を図 2 に示す。実軸と虚軸に= 平行な直線で 250000(= 500 × 500) の長方形に分割(ν)−1.0 + 1.0i を用いてニュートン法を適用する。計算 し頂点を初期値としてニュートン法を適用し z がは gcc 4.0.1 の double complex(実部虚部とも 10 進 周期点に近づいたら初期値をプロットしている。
例 4 例 1 の関数 f (z) = z − 2z − 5 に初期値 z3(0)16 桁弱) を用いた。

80(ν)

νz

0−1.00000000 + 1.00000000i|f (z(ν)

)|

60

1.00

40−1123−1.05000000 + 1.15000000i−1.04726405 + 1.13606317i−1.04727573 + 1.13593990i

1.09 × 10
9.40 × 10−4
7.11 × 10−84−1.04727574 + 1.13593989i

1.78 × 10−1520

0−20−40−1.047275740771163 + 1.135939889088928i に 2次収束している。

−60−80−100−80−60−40−20020406080100

図 2: 周期点に収束する初期値−100 5 Re z 5 100, −80 5 Im z 5 806図 3 には −10 5 Re z 5 10, −8 5 Im z 5 8 の[7] 関孝和・建部賢明・建部賢弘、大成算経、巻三、東北大学、狩野 7.20820.20

範囲を示す。図 2,3 はほぼ相似であることが見て取れる。

http://www2.library.tohoku.ac.jp/wasan/[8] 杉原正顯・室田一雄、数値計算法の数理、岩波書店、1994864[9] 上田哲生・谷口雅彦・諸澤俊介、複素力学系序説、
培風館、1995

20[10] 山本哲朗、数値解析入門 [増訂版]、サイエンス社、2003−2−4−6[11] T.J. Ypma, Development of the Newton-−8−10−8−6−4−202468Raphson Method, SIAM Review, 37(1995), 531551.

10

図 3: 周期点に収束する初期値−10 5 Re z 5 10, −8 5 Im z 5 8

おさだ なおき (東京女子大学)

丸め誤差がないと仮定したとき、収束しない初期値の集合が (有限集合、可算集合、代数曲線の部分集合、測度 0 の集合、など) 無視できるとき、反復法は大域的収束性を持つという。図 2 から、ニュートン法は大域的収束性を持たないことが分かる。なお、フラクタルや大域的収束性は多様な定義がある。

参考文献[1] F. Cajori, Historical Note on the NewtonRaphson Method of Approximation, Amer.
Math. Monthly, 18(1911), 29-32.
[2] K. ファルコナー、服部・村井訳、フラクタル幾
何学、共立出版、2006[3] 平山諦・下平和夫・広瀬秀雄、関孝和全集、大阪
教育図書、1974[4] 伊理正夫・藤野和建、数値計算の常識、共立出
版、1985[5] I. Newton, De analysi per æquationes numero
terminorum inﬁnitas, 257-282, Opera quae exstant omnia, Bd 3, F. Fromann, 1964
[6] E. Schröder, Ueber unendlich viele Algorithmen
zur Auﬂösung der Gleichungen, Math. Annal.,2(1870), 317-365.

7