非線型方程式、特に代数方程式の解法

桂田祐史

2007 年 3 月 24 日

目次

第 1 章非線型方程式概説 4 1.1 不動点定理に基づく反復法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Newton 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.2.1 基礎事項 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2 Newton 法を用いた初等的な計算 . . . . . . . . . . . . . . . . . . . . . 9 1.2.3 停止則について . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.2.4 Kantorovich の定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2.5 減速 Newton 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2.6 実例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.3 連続変形法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4 その他 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

第 2 章代数方程式の解法 14

2.1 序 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.1 落書き . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.2 1 次元代数方程式の解法概観 . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1.3 連立代数方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.2 連立法, 特に Durand-Kerner 法 . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.2.1 Durand-Kerner 法 — 名前の由来 . . . . . . . . . . . . . . . . . . . . . 18 2.2.2 DK 法 — Durand の解釈 . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2.3 復習: 根と係数の関係 . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.2.4 DK 法 — Kerner の解釈 . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.2.5 解の精度の検証 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2.6 DK 法の長所・短所 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2.7 初期値の取り方の重要性 . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2.8 Aberth の初期値, DKA 法 . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2.9 Ehrlich-Aberth 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2.10 根の大きさの限界 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.11 伊理の初期値 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

付録 A おもちゃ箱 (その他の方法) 26 A.1 平野法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 A.2 Horner 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 A.3 Graffe 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 A.4 Bairstow Hitchcock の方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 A.5 Strum の方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

A.5.1 スツルムの定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 A.5.2 ユークリッドの互除法による Strum 列の生成 . . . . . . . . . . . . . . 32 A.5.3 3 重対角行列の固有多項式と Strum 列 . . . . . . . . . . . . . . . . . . 34 A.5.4 直交多項式の作る Strum 列 . . . . . . . . . . . . . . . . . . . . . . . . 36 A.5.5 一般化された Strum 列 . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

A.6 Bernoulli 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 A.7 商差法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 A.8 次数低下 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 A.9 Jenkins and Traub 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 A.10 その他 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

A.10.1 Euler 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 A.10.2 ハレー法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 A.10.3 ポメンタール法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

付録 B 応用解析 IV の数値実験 43 B.1 6701 号室のワークステーションの利用 . . . . . . . . . . . . . . . . . . . . . . 43 B.2 応用解析 IV ホームページ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 B.3 C++ について . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

B.3.1 まずはサンプル・プログラム . . . . . . . . . . . . . . . . . . . . . . . 44 B.3.2 C++ の (とても簡潔な) 紹介 . . . . . . . . . . . . . . . . . . . . . . . 45 B.3.3 サンプル・プログラムを読むための C++ の知識 . . . . . . . . . . . . 45 B.3.4 g++ — 実験に用いる C++ 処理系 . . . . . . . . . . . . . . . . . . . . 46

B.4 グラフィックス・ライブラリィ GLSC について . . . . . . . . . . . . . . . . . 46 B.5 素朴な Newton 法の実行例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 B.6 DK 法のサンプル・プログラム . . . . . . . . . . . . . . . . . . . . . . . . . . 49

B.6.1 ソースプログラム DK3.C . . . . . . . . . . . . . . . . . . . . . . . . . . 49 B.6.2 コンパイル＆実行例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

付録 C 情報処理 II から 54 C.1 レポート課題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 C.2 例題を解くプログラム例 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

C.2.1 Newton 法の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 C.2.2 二分法の場合 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

付録 D 多次元の Newton 法の例 57 D.1 ターゲット問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 D.2 復習 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 D.3 ターゲット問題の差分近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 D.4 Newton 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 D.5 サンプル・プログラム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

メモ

まずは Kantorovich の定理, それから山本 [15] を勉強することだ。

• 2003 年度 Gröbner 基底がらみで面白い本 (齋藤・竹島・平野 [3]) が出た。誰か卒研でやらないかな。

• 連立法が Newton 法であることの証明は、山本 [15] にある。これを読解して解説するようにしたい。

• Kantorovich の定理の証明を消化して書いておこう。

• 篠原 [4] にある占部^{1}の定理。

1篠原能材氏の恩師故占部実氏のこと。

第1章非線型方程式概説

有限次元の線型方程式については、問題の規模がさほど大きくなければ十分に満足できる解法があるが、非線型方程式一般に適用できる解法はない。

• 不動点定理に基づく反復法

• Newton 法

1.1 不動点定理に基づく反復法

与えられた方程式をそれと同値な

(1.1)  x = F (x)

に変換してそれを解くと解釈できる方法が多い。

方程式 (1.1) の解 x のことを F の不動点と呼び、方程式 (1.1) の解の存在を保証する定理を不動点定理と呼ぶ。

ここでは (1.1) のタイプの方程式を不動点型の方程式と呼ぶことにする。

適当に選んだ X の要素 x_{0} から漸化式

(1.2)  x_{k+1} = F (x_{k}) (k = 0, 1, 2, ⋯)

で定めた列 {x_{k}}k ∈ N が極限 x_{∞} ∈ X を持つことがある。 F が連続であれば x_{∞} は F の不動点である (これは (1.2) において k → ∞ の極限を取れば良い)。このことを背景に列 {x_{k}}_{k∈N} を実際に計算して、十分大きな番号 k に対する x_{k} を方程式の近似解として採用する方法を反復法と呼ぶ。

[figure]

定理 1.1.1 (Banach の不動点定理 (縮小写像に関する不動点定理))(X, d) は完備距離空間で、 F: X → X を縮小写像とする。すなわち

d(F (x), F (y)) ≤ kd(x, y) (x, y ∈ X)

を満たす k ∈ [0, 1) が存在すると仮定する。このとき F の不動点が一意的に存在し、それは X の任意の要素 x_{0} を初期値とした反復法

x_{k+1} = F (x_{k}) (k = 0, 1, 2, ⋯)

で生成される列 {x_{k}} の極限として得られる。

[figure]

証明有名なので省略する。例えば Schwartz [23] を参照せよ。注意 1.1.1 Lipschitz 定数 < 1 の Lipschitz 条件を満たす写像は縮小写像である。ゆえに f が ∥f′(x)∥ ≤ ∃L < 1 を満たすならば縮小写像である。

[figure]

注意 1.1.2 (Banach の不動点定理の簡単な拡張)F そのものでなく、適当な自然数 p に対して F^{p} が縮小写像になる場合にも、 F の不動点は一意的に存在し、それは反復法で得られる。ただし F^{p} とは

F^{1} =^{def.} F, F^{m+1} = F^{m} ∘ F (m ∈ N)

で帰納的に定義される写像である。

[figure]
[figure]

定理 1.1.2 (Brouwer の不動点定理)D を R^{n} の有界閉凸集合とし、 f: D → D が連続とするとき、 f は少なくとも一つの不動点を持つ。

[figure]

証明 Zeidler [24], 増田 [10] 等を参照せよ。

1.2 Newton 法

山本哲朗先生の一連の著作、特に山本 [15] を見よ。

1.2.1 基礎事項

X を Banach 空間、 U を X の開集合、 f: U → X を Fréchet 微分可能な写像とするとき方程式

(1.3)  f(x) = 0

を考える。適当な x_{0} ∈ U を取って漸化式

(1.4)  x_{k+1} = x_{k} − (f′(x_{k}))^{−1} f(x_{k})

で列 {x_{k}}_{k∈N} を生成したとき、これが極限 x_{∞} を持つことがある。このとき x_{∞} は (1.3) の解となる。

以下 (f′(x_{k}))^{−1} のことを f′(x_{k})^{−1} と書く。

そこで十分大きな番号 k を取って、 x_{k} を (1.4) の近似解として採用する方法が考えられるが、それを Newton 法と呼ぶ。形式的には

F (x) := x − f′(x)^{−1}f(x)

で定められる F の不動点を反復法で求めていることになる^{1}。注意 1.2.1 (Newton 法という名前について)Newton が 3 次方程式の近似解法に利用した

1F の定義域を考えると、 Banach の不動点定理をすぐに適用できるような状況にはないことが分かる。しかし次のことに注意しよう。簡単のため、 1 次元とする。 F′(x) = f(x)f′′(x)/f′(x)^{2} なので、解の十分近くでは |F′(x)| ≤ ∃L < 1. これから解を中心とする十分小さな閉球を X とすると、 F |_{X} は縮小写像である。

ことに基づくそうだが、慣用の形に一般化したのは Raphson なので、 Newton-Raphson^{ラフソン}法と呼ばれることが多い。

注意 1.2.2 (準 Newton 法) f の値の計算コストが高い場合は、 (1.4) でなく

x_{k+1} = x_{k} − f′(x_{0})^{−1}f(x_{k})

で列 {x_{k}}_{k∈N} を生成することがある。この方法を準 Newton 法, 簡易 Newton 法, von Mises^{フォン・ミーゼ} 法、あるいは修正 Newton 法などと呼ぶ。これは線形収束となる。

[figure]

定義 1.2.1 (p 次収束) p ≥ 1 に対して、 {x^{(k)}}_{k∈N} が a に p 次収束するとは、正定数 L が存在して (p = 1 のときは L < 1 を要請する)、十分大きい任意の k に対して

∥x^{(k+1)} − a∥ ≤ L∥x^{(k)} − a∥^{p}

が成り立つことをいう。特に p = 1 のとき、線形収束ともいい、 L を収束率と呼ぶ。なお、

lim_{k→∞} \frac{∥x^{(k+1)} − a∥}{∥x^{(k)} − a∥} = 0

のとき、超 1 次収束 (superlinear convergence) するという。

[figure]

線形収束の場合、

∥x^{(k)} − a∥ ≤ L^{k}∥x^{(0)} − a∥

が成り立つ。つまり等比数列的 (あるいは指数関数的) に収束するわけである。

[figure]

定理 1.2.1 (Newton) Ω を R^{n} の開集合、 f: Ω → R^{n} を C^{2}-級の写像、 a ∈ Ω, そして f(a) = 0, f′(a) が正則とするとき、 a に十分近い x_{0} から始めて、 Newton 法で作った列 {x_{k}}_{k∈N} は a に 2 乗収束する。すなわち正定数 C が存在して

∥a − x_{k+1}∥ ≤ C∥a − x_{k}∥^{2} (k ∈ N). [figure]

証明 (概略) f は C^{2}-級で、 f(a) = 0 と仮定すると、十分 a に近い x_{k} に対して

x_{k+1} − a = [x_{k} − f′(x_{k})^{−1}f(x_{k})] − a
= x_{k} − a − f′(x_{k})^{−1}f(x_{k})
= x_{k} − a − f′(x_{k})^{−1} (f(x_{k}) − f(a))
= x_{k} − a − f′(x_{k})^{−1} [f′(x_{k})(x_{k} − a) − O ((a − x_{k})^{2})]
= O ((a − x_{k})^{2}) .

詳しくは例えば山本 [15] を見よ。

重根の場合は収束の速さが線形収束になり (つまり「収束が遅い」 )、実際の数値計算で得られる最終的な精度も低い。注意 1.2.3 上の定理では舞台を R^{n} としたが、 Banach 空間 X でも大丈夫である。それには Banach 空間 X の開集合 Ω で定義された写像 f: Ω → X の 2 回微分可能性を定義する必要があるが、 f が x^{∗} で Fréchet 微分可能で、線型写像

T: X → L(X, X) := {φ; φ: X → X 連続線形}

が存在して、

f(x) = f(x^{∗}) + f′(x^{∗})(x − x^{∗}) + T (x − x^{∗})(x − x^{∗}) + S(x; x^{∗})

によって S(x; x^{∗}) を定義するとき

lim_{x→x^{∗}} \frac{∥S(x; x^{∗})∥}{∥x − x^{∗}∥^{2}} = 0

が成り立つとき、 f は x^{∗} で 2 回強微分可能、 T を 2 階Fréchet 微分という。 T は T: X×X → X という有界双線型写像とも見なせる。 — ここの説明は一松によるものだが、 Schwartz [23] も見ておこう。

微積分の演習問題から

次の命題は微積分の演習書によく登場する。

[figure]

命題 1.2.1 (凸関数の Newton 法) f: [a, b] → R が 2 回微分可能で、

f′′(x) > 0 (x ∈ [a, b]), f(a)f(b) < 0

を満すならば、方程式 f(x) = 0 の解 s は (a, b) 内にただ一つだけ存在し、

x_{0} = \{ \begin{matrix}b & (f(b) > 0 のとき) \\ a & (f(a) > 0 のとき) \end{matrix}

を初期値とする Newton 法で定まる列 {x_{n}}_{n∈N} の極限となる:

lim_{n→∞} x_{n} = s. [figure]

証明 (それなりに長い)

1. (存在) f(s) = 0 となる s ∈ (a, b) が存在することは中間値の定理から分かる。

2. (一意性)

f(s_{1}) = f(s_{2}) = 0, a < s_{1} < s_{2} < b

を満たす s_{1}, s_{2} が存在したと仮定すると、 Rolle の定理から

∃z ∈ (s_{1}, s_{2}) s.t. f′(z) = 0.

仮定 f′′ > 0 から、

f′ < 0 ((a, z) 内で) ∴ f(a) > f(s_{1}) = 0,
f′ > 0 ((z, b) 内で) ∴ f(b) > f(s_{2}) = 0.

これは f(a)f(b) < 0 に反する。ゆえに f(s) = 0 となる s は unique である。

3. (Newton 法の収束) 以下、 f(b) > 0 の場合に x_{0} = b を初期値とする Newton 法が収束することを示す (f(a) > 0 の場合も同様にできる)。 f(x) = 0 の unique な解を s とすると、

(1.5)  f(x) > 0 (x ∈ (s, b]),
(1.6)  f′(x) > 0 (x ∈ [s, b])

が成り立つ。

(1.5) の証明もしも f(x^{∗}) ≤ 0, x^{∗} ∈ (s, b] となる x^{∗} があれば、中間値の定理より f(x^{∗∗}) = 0, x^{∗∗} ∈ [x^{∗}, b) となる x^{∗∗} が存在することになるが、これは f(x) = 0 の解の一意性に反する。

(1.6) の証明まず f(s) = 0, h > 0 のとき f(s + h) > 0 であるから、

f′(s) = lim_{h↓0} \frac{f(s + h) − f(s)}{h} ≥ 0.

ところが、もし f′(s) = 0 ならば f′′ > 0 より、 f′(x) < 0 (x ∈ (a, s)). ゆえに f(a) > f(s) = 0. これは矛盾であるから f′(s) ≠ 0. ゆえに f′(s) > 0. 仮定 f′′ > 0 から (1.6) が導かれる。

さて、

x_{0} = b, x_{n+1} = x_{n} − \frac{f(x_{n})}{f′(x_{n})}

で数列 {x_{n}}_{n∈N} を決めようとするわけであるが、次の不等式が成り立つことが分かれば、 {x_{n}} は well-defined で、 s を下界とする単調減少数列であることが分かる。

(1.7)  x ∈ (s, b] ⟹ s < x′ < x, where x′ := x − \frac{f(x)}{f′(x)}.

(1.7) の証明既に示したように f(x) > 0, f′(x) > 0 より

\frac{f(x)}{f′(x)} > 0 ∴ x′ = x − \frac{f(x)}{f′(x)} < x.

さらに

x′ − s = x − \frac{f(x)}{f′(x)} − (s − \frac{f(s)}{f′(x)}) (∵ f(s) = 0)
= (x − s) − \frac{f(x) − f(s)}{f′(x)}
= \frac{f′(x)(x − s) − (f(x) − f(s))}{f′(x)}
= \frac{f(s) − [f(x) + f′(x)(s − x)]}{f′(x)}
= \frac{f′′(x + θ(s − x))(s − x)^{2}/2}{f′(x)}, ∃θ ∈ (0, 1) (∵ Taylor の定理)
> 0.

ラスト {x_{n}}_{n∈N} は s を一つの下界とする単調減少数列であるから、収束する。その極限を x_{∞} とすると、

x_{n+1} = x_{n} − \frac{f(x_{n})}{f′(x_{n})}

で n → ∞ として

x_{∞} = x_{∞} − \frac{f(x_{∞})}{f′(x_{∞})}.

これから f(x_{∞}) = 0. 解の一意性から s = x_{∞}. すなわち

lim_{n→∞} x_{n} = s.

T 先生の講義の思い出 Newton 法の説明をしていただいたときのこと。 Newton 法はいつも収束するわけではないと言って、反例の存在を図で示そうとし始めたのだが、なかなかうまく行かない。見ると関数の二階導関数の符号が一定のグラフばかり書いている。微積分の勉強で上の命題を知っていた私は、少しどきどきしながら「変曲点がないと…」とつぶやいた。無事に T 先生の耳に届き、すぐに了解して反例のグラフを描かれました。めでたし、めでたし。

1.2.2 Newton 法を用いた初等的な計算

例 1.2.1 (平方根) a > 0 に対して、\sqrt{a} を f(x) = x^{2} − a = 0 の根と考えて、 x_{0} > 0 から始めて

(1.8)  x_{k+1} = x_{k} − \frac{x_{k}^{2} − a}{2x_{k}}

により列 {x_{k}}_{k∈N} を生成すると、

x_{k+1} − \sqrt{a} = \frac{(x_{k} − \sqrt{a})^{2}}{2x_{k}} .

なお (1.8) は、

x_{k+1} = \frac{1}{2} (x_{k} + \frac{a}{x_{k}} )

のように整理することもできるが、数値計算においては、丸め誤差の観点から不利である。例 1.2.2 (平方根の逆数) 1/\sqrt{a} を f(x) = 1/x^{2} − a = 0 の根と考えて、

x_{k+1} = x_{k} + \frac{(1/x_{k})^{2} − a}{2/(x_{k})^{3}} = x_{k} (\frac{3}{2} − \frac{a}{2}(x_{k})^{2})

あらかじめ a/2 などを計算しておけば、除算なしで、乗算 3 回だけで計算可能である。例 1.2.3 (立方根) \sqrt[3]{a} を f(x) = x^{3} − a = 0 の根と考えると、

x_{k+1} = \frac{2}{3} (x_{k} + \frac{a}{(x_{k})^{2}} )

となるが、 f(x) = x^{2} − \frac{a}{x} = 0 の根と考えると

x_{k+1} = x_{k} \frac{(x_{k})^{3} + 2a}{2(x_{k})^{3} + a}

となる。例 1.2.4 (逆数) 1/a を f(x) = a − 1/x = 0 の根と見なすと

x_{k+1} = x_{k}(2 − ax_{k})

となる。

 newton.c

[figure]

例 1.2.5 (プログラム例) /*

 * newton.c -- Newton 法で方程式 f(x)=0 を解く
 * コンパイル: cc -o newton newton.c -lm
 * 実行: ./newton
 */
#include <stdio.h>
#include <math.h>
int main()
{
    int i, maxitr = 100;
    double f(double), dfdx(double), x, dx, eps;
    printf(" 初期値 x0, 許容精度ε=");
    scanf("%lf %lf", &x, &eps);
    for (i = 0; i < maxitr; i++) {
        dx = -f(x) / dfdx(x);
        x += dx;
        printf("f(%24.15e)=%9.2e\n", x, f(x));
        if (fabs(dx) <= eps)
            break;
    }
    return 0;
}
double f(double x)
{
    return cos(x) - x;
}

/* 関数 f の導関数 (df/dx のつもりで名前をつけた) */

double dfdx(double x)
{
    return -sin(x) - 1.0;
}

[figure]

1.2.3 停止則について

小さな正の数 ε や δ を与えて

|Δx_{k}| < ε (ただし Δx_{k} := x_{k+1} − x_{k})

あるいは

|f(x_{k})| < δ

のとき、 x_{k} が十分良い近似解となったと判断して反復を停止するという考えで書かれているプログラムが多いが、 ε や δ をどう取るべきかまで考えておかないと不十分であろう。

例えば代数方程式

∑_{i=0}^{n} a_{i}x^{i} = a_{n}x^{n} + a_{n−1}x^{n−1} + ⋯ + a_{1}x + a_{0} = 0

を Newton 法で解く場合には、 ε_{M} を計算機イプシロンとして、

|f(x_{k})| ≤ (∑_{i=0}^{n} |a_{i}||x_{k}|^{i}) ε_{M}

が良いと言われている (杉原・室田 [5] p.64 を見よ)。

1.2.4 Kantorovich の定理

Newton 法の収束域、収束速度を与えるだけでなく、解の存在そのものまで保証できる優れた定理である。

h < 1/2 が単根の場合、 h = 1/2 が重根の場合に相当する。

またこの定理の仮定は、精度保証付きの計算で確認可能な可能性が高いものであることに注意しよう。

[figure]

定理 1.2.2 (Newton-Kantorovich) D を R^{n} の開凸集合とし、 f: D → R^{n} は微分可能で、 x^{(0)} ∈ D で、次の条件が満たされるとする。

(a) Jacobi 行列 f′ は D で Lipschitz 連続、すなわち

∃L > 0 s.t. ∥f′(x) − f′(y)∥ ≤ L∥x − y∥ (x, y ∈ D)

が成り立つ。

(b) f′(x^{(0)}) は正則行列で、

∥f′(x^{(0)})^{−1}∥ ≤ a, ∥f′(x^{(0)})^{−1}f(x^{(0)})∥ ≤ b, h := abL ≤ \frac{1}{2}.

(c) t^{∗} := (1 − \sqrt{1 − 2h})/(aL), U := B(x^{(0)}; t^{∗}) とするとき、

\overline{U} ⊂ D.

このとき次の (1), (2) が成り立つ。

(1) \overline{U} の中に f(x) = 0 の解 x^{∗} がただ一つ存在する。

(2) x^{(0)} を初期値とする Newton 法による列 \{x_{k}\}_{k∈N} が定義され、 x^{(k)} ∈ \overline{U} であり、

∥x^{(k)} − x^{∗}∥ ≤ \frac{(1 − \sqrt{1 − 2h})^{2^{k}}}{2^{k}aL} (k = 0, 1, 2, ⋯). [figure]

証明杉原・室田 [5] pp.71–75 を見よ (これは Zeidler [24] の証明を改良したものだとか)。他にも山本 [16] や、 Rall [22] に別証があるとか。

この証明を読んでみよう。それから、この定理を利用して精度保証するプログラムを書いてみよう。

1.2.5 減速 Newton 法

杉原・室田 [5] を見よ。

1.2.6 実例

杉原・室田 [5] の例 4.1, 4.2 など。

1.3 連続変形法

杉原・室田 [5] の §4.4 を見よ。

1.4 その他

• 歴史

• Kantorovich の定理の証明

• Aitken 加速, Steffensen 反復

第2章代数方程式の解法

2.1 序

2.1.1 落書き

この章に書いてあるのは、あくまでも桂田個人の学習ノートであり、書きかけのものも多く、単なる写経も多い。

まじめにやる必要がある場合 (実際に問題を解く、数値解析の講義をする場合) は、一松 [8], 杉原・室田 [5], 伊理・藤野 [2] などを一読すること。

杉原・室田 [5] はきちんとまとまっていて面白い。読み物としては一松 [8] が面白い。山本・北川 [17] は数値例が豊富で楽しい。

一松先生曰く

• 古くて新しい難問である。 4000 年以上昔のバビロニアの粘土板にもあるくらい歴史があるが、コンピューターの時代になってかえって難しくなったとも言える (扱える数の精度が低く、範囲が狭い)。

• すべての問題をたった一つの方法で解くのは困難である。すくなくとも「序盤向き」の解法と「終盤向き」の解法を使い分けるのがよい。

序盤向き精度は低くても、根の位置をつきとめ、他の根と分離する。

二分法 (Strum 列を利用した二分法)、商差法終盤向き大域収束性はなくとも、十分よい近似値から始めれば速く収束する。

Newton 法、 Bairstow 法、 Durand-Kerner 法

伊理・藤野 [2] 「数値計算の常識」では

Newton 法だ、と。

杉原・室田 [5] 「数値計算法の数理」では

• 平野法、連立法 (Durand-Kerner 法) について論じている。

2.1.2 1 次元代数方程式の解法概観

(今のところ、この項はガラクタ箱になってしまっている。 )

2 次方程式

2 次方程式 ax^{2} + bx + c = 0 は根の公式

x = \frac{−b ± \sqrt{b^{2} − 4ac}}{2a}

がある。高校数学では係数 a, b, c は実数であるが、複素数にしてもこの式のままで大丈夫である (なお、複素数の平方根は、実数の範囲で計算できる — Ahlfors の有名な複素解析の教科書に説明がある)。

この根の公式のまま計算しようとすると、 b^{2} ≫ 4|ac| の場合に桁落ち (cancellation) が起きてしまい、精度が低下してしまう。

例えば b > 0 の場合、

x_{1} = \frac{−b + \sqrt{b^{2} − 4ac}}{2a}

の計算で桁落ちが起こるので、

x_{1} = \frac{−b + \sqrt{b^{2} − 4ac}}{2a} = \frac{(−b + \sqrt{b^{2} − 4ac})(−b − \sqrt{b^{2} − 4ac})}{2a(−b − \sqrt{b^{2} − 4ac})} = \frac{−2c}{b + \sqrt{b^{2} − 4ac}}

のように分子の有理化を行うか、桁落ちの起こらない方の根

x_{2} = \frac{−b − \sqrt{b^{2} − 4ac}}{2a}

を先に計算してから、根と係数の関係 x_{1}x_{2} = c/a から

x_{2} = \frac{c}{ax_{1}}

と計算する (この式の右辺を整理すると上と同じになる)。

3 次方程式

いわゆる Cardano の方法がある^{1}。かつてはこの方法は数値計算には向かないと考えられていたが (例えば森口 [14])、最近では見直されているそうである ( 「特に 1 実根 2 複素根を持つ実 3 次方程式の標準算法として見直されてきている」 — 一松 [8])。

まず

x^{3} + a_{2}x^{2} + a_{1}x + a_{0} = 0

を

x + \frac{a_{2}}{3} = X

と置き換えて、

(2.1)  X^{3} − 3px + 2q = 0, p = (\frac{a_{2}}{3} )^{2} − \frac{a_{1}}{3} , q = \frac{a_{0}}{2} − \frac{a_{1}a_{2}}{6} + (\frac{a_{2}}{3} )^{3}

1Cardano (1501–1576). 『偉大な学術すなわり代数学の法則について』の中で\sqrt{−1} という記号を使った。「10 をその積が 40 になる 2 つの数に分解すると 5 ± \sqrt{−15} になる」。

に変換する。次に

X = u + v, uv = p

とすると、 (2.1) は

u^{3} + v^{3} = −2q, u^{3}v^{3} = p^{3}

となる。ゆえに

(2.2)  t^{2} + 2qt + p^{3} = 0

を解けば u^{3}, v^{3} を得る。実行すると

u^{3}, v^{3} = −q ± \sqrt{q^{2} − p^{3}}.

ここで場合わけをする。

(I) もし q^{2} ≥ p^{3} ならば、 (2.2) は 2 実根

t_{1} = −q + \sqrt{q^{2} − p^{3}}, t_{2} = −q − \sqrt{q^{2} − p^{3}}, ω = e^{2πi/3} = \frac{−1 + \sqrt{3}i}{2}

を持つ。その三乗根 (実数！ ) から

x_{1} = t_{1}^{1/3}ω + t_{2}^{1/3}ω^{2} − \frac{a_{2}}{3} ,
x_{2} = t_{1}^{1/3}ω^{2} + t_{2}^{1/3}ω − \frac{a_{2}}{3} ,
x_{3} = t_{1}^{1/3} + t_{2}^{1/3} − \frac{a_{2}}{3} .

(II) q^{2} < p^{3} ならば、 t^{1/3} の計算には、複素数の 3 乗根を要する。注意すべきは t_{1}^{1/3}t_{2}^{1/3} = uv = p であるから、 t_{1}^{1/3} を求めるごとに対応する t_{2}^{1/3} を

t_{2}^{1/3} = \frac{p}{t_{1}^{1/3}}

で定める。これらを数値計算するには極座標に直して計算すればよい。この場合は、代数方程式は 3 実根を持つが、それは実数の四則と累乗根だけでは求められないことが証明されている (不還元の場合と呼ばれる)。

3 次方程式の数値計算の際の注意事項一松 [8] の §21 を見よ。なお、渡辺・名取・小国 [18] には Fortran サブルーチンが載っている。

実根のみを持つ実係数代数方程式の解きにくさ

例 2.1.1 (Wilkinson の例)

p(x) = ∏_{j=1}^{20}(x − j) = (x − 1)(x − 2) ⋯ (x − 19)(x − 20) = x^{20} − 210x^{19} + ⋯ + 20! = 0

を考える。係数が変化したとする。これは q(x) を定まった方程式、 ε を小さな正数として

p(x) + εq(x) = 0

を解くことになる。その結果、根 x_{0} が x_{0} + δ に変わったとする。

0 = p(x_{0} + δ) + εq(x_{0} + δ) = p(x_{0}) + δp′(x_{0}) + ⋯ + εq(x_{0}) + εδq′(x_{0}) + · · · .

大ざっぱに

δ ≒ εq(x_{0})/p′(x_{0}).

q(x) = x^{19} とすると、

δ ≒ Cε, C = ±\frac{x_{0}^{19}}{(x_{0} − 1)!(20 − x_{0})!}.

C は x_{0} = 16 に対して最大で、そのとき

|C| ≒ 2.4 × 10^{9}.

Wilkinson は x^{19} の係数 210 が 210 + ε に変わったとき、 x_{0} = 16 に ε のほぼ 24 億倍の誤差が入ることを確かめている。 ε = 2^{−23} ≒ 10^{−7} とすると、もはや ε の 2 乗以上の項が無視できない。実際にはこのときには、 20 個の実根のうち 10 個が, 5 対の複素根に化けてしまい、その虚数部は最大 2.8 にもなっている。

2.1.3 連立代数方程式

一松 [9] 第 3 章 6 節「連立代数方程式」の写経。

計算機代数の話題はまだ数多くあるが、近年有名になった「グレブナー基底」について略説する。多変数多項式に関する連立代数方程式を解く問題は良く現れるが、一般的な解法はほとんどない。理

論上ではつぎつぎに消去法によって変数を消去してゆけば解けるはずだが、それを正直に実行すると、たちまち大きな係数を持つ高次方程式が現れて、手におえなくなる。

これまでは与えられた方程式を適当に組み合わせて、因数分解できる形に直す方法が主力だった。しかしこれを下手に実行すると、考察すべき場合の個数が激増し、解をうまくまとめるのが厄介いになる。

変数に対称性がある場合には、それを生かしてうまく解ける場合もあった。ただしそれにこだわるのは得策でなく、対称性を無視して消去法を使うほうがよいことも多い。

この問題は、一般化すると、多項式イデアルに関する所属判定問題になる。イデアルとは、歴史的な由来のある述語だが、いくつかの基底要素 p_{1}, ⋯, p_{m} に多項式 u_{1}, ⋯, u_{m} をかけて加えて表わされる

I = {u_{1}p_{1} + ⋯ + u_{m}p_{m}}

という形の多項式の族と考えてよい。

一変数の場合には、このようなイデアルが、 p_{1}, ⋯, p_{m} の最大公因子 d の倍式全体として表される (単項イデアル)。この d は「よい基底」である。

多変数の場合には、単項でないイデアルがふつうであり、「よい基底」を求めることが重要である。連立代数方程式 p_{1} = 0, ⋯, p_{m} = 0 を解くには、原理的には p_{1}, ⋯, p_{m} の生成するイデアルの任

意の基底の共通零点を求めればよい。しかし基底が悪いと、この操作が実行困難である。

イデアルに関する所属判定問題とは、与えられた多項式 f が、そのイデアル I に含まれるかどうかを判定する算法を問う問題である。その応用例は後述する。もし具体的に u_{1}, ⋯, u_{m} を求めて

f = u_{1}p_{1} + ⋯ + u_{m}p_{m}

と書くことができれば f は I に含まれる、他方もし変数 x_{1}, ⋯, x_{n} にたいして適当な値 a_{1}, ⋯, a_{n} で

p_{1}(a_{1}, ⋯, a_{n}) = 0, ⋯, p_{m}(a_{1}, ⋯, a_{n}) = 0, f(a_{1}, ⋯, a_{n}) ≠ 0

である組がみつかれば、 f はイデアル I には含まれない。

その昔、何日も所属問題がうまく解けずに苦しみ、試しに数値を代入したら、前期のような組がみつかり、 f が I に属さないことがわかった、そして結果的にそれまでの計算はむだだった、という実例がある。現在でもこのような予備テストは重要である。しかしこのような a_{1}, ⋯, a_{n} の組がみつからない (いつでも f(a_{1}, ⋯, a_{n}) = 0 となる) とき、 f が I に含まれそうだと予想されても証明にはならない。

この判定には I の基底自体が重要である。じっさいに基底のとり方が悪いと、奇妙な現象が起こる。例えば変数の個数を 3 とし、 x_{1}, x_{2}, x_{3} を x, y, z と書くことにして、次の基底を考える。

p_{1} = x^{3}yz − xz^{2}, p_{2} = xy^{2}z − xyz, p_{3} = x^{2}y^{2} − z^{2}

これにたいして f = x^{2}yz − z^{3} の所属問題を考えてみよう。

原則として最も高い次数の項を合わせて消去する。しかし f の主項 x^{2}yz は、 p_{1}, p_{2}, p_{3} のどれを使っても、そのままでは簡約できない。この f は^{x}p_{2} + zp_{3} と表現できるので、イデアルに属する。このことは x^{2}y^{2}z というおなじ項を加減して変形すればわかるが、これを機械的に求めることはできない。

この例は f が悪いというよりも、イデアルの基底が不足しているのである。これにさらに新しい基底を補充して、基底 (base) を固める必要がある。

多項式イデアルの基底を補充して、標準的な基底を整備する理論は、かなり昔から考えられていた。たとえば広中平祐氏が、フィールズ賞の対象となった「特異点の還元」の研究中に、その種の「標準基底」を扱っている。したがって、これを「広中基底」とよぶべきかもしれない。

しかしその具体的な構成算法を示したのはブーフバーガーであり、彼は恩師の名をとって、それをグレブナー基底とよんだ。今日その名で広く使われているので、以下でもそうよんで論じることにする。

それはある意味では、連立一次方程式を解くガウスの消去法の一般化である。対称式については基本対称式が「グレブナー基底」に相当する。一般のグレブナー基底は、これらの概念の一般化に相当する。

2.2 連立法, 特に Durand-Kerner 法

n 次代数方程式を n 個の初期値 {x_{i}^{(0)}}_{i=1}^{n} から始めて、

x_{i}^{(k+1)} = x_{i}^{(k)} + Δx_{i}^{(k)} (i = 1, …, n)

の形の漸化式で計算していく方法を連立法と呼ぶ。 2 次の連立法 (Durand-Kerner 法)、 3 次の連立法 (Ehrlich-Aberth 法) を解説する。

2.2.1 Durand-Kerner 法 — 名前の由来

実際には古くからあるとも言われているが (山本 [15])、フランスの Durand が古く提唱し、ドイツの Kerner が別の立場から解釈を与えた (1966 年) ため、今では Durand-Kerner 法と呼ばれるのが普通である、と一松先生は書いているが、杉原・室田 [5] では「2 次の連立法」と呼んでいる。

2.2.2 DK 法 — Durand の解釈

代数方程式

p(x) ≡ x^{n} + a_{n−1}x^{n−1} + ⋯ + a_{1}x + a_{0} = 0

の根を ξ_{1}, ⋯, ξ_{n} とすると

p(x) = (x − ξ_{1}) ⋯ (x − ξ_{n}),
(2.3)  p′(ξ_{i}) = ∏_{j≠i}(ξ_{i} − ξ_{j}) = (ξ_{i} − ξ_{1}) ⋯ (ξ_{i} − ξ_{i−1})(ξ_{i} − ξ_{i+1}) ⋯ (ξ_{i} − ξ_{n}).

Newton 法では現在の近似値 x^{(k)} から次の近似値 x^{(k+1)} を求めるには

(2.4)  x^{(k+1)} = x^{(k)} − \frac{p(x^{(k})}{p′(x^{(k)})}

とする。ここで導関数の計算をしないで済ませるために、この分母を (2.3) で代用することを考える。つまり x^{(k)} が ξ_{i} の近似値だと考えて p′(x^{(k)}) を p′(ξ_{i}) で置き換えるわけである。各根 ξ_{i} に対する近似値 x_{i}^{(k)} が i = 1, ⋯, n について全部そろっているとして、 (2.4) の代わりに

(2.5)  x_{i}^{(k+1)} = x_{i}^{(k)} − \frac{p(x_{i}^{(k)})}{∏_{j≠i} (x_{i}^{(k)} − x_{j}^{(k)})}.

を用いる。 Durand は (2.5) の収束性を示し、最終的に 2 乗収束になることを証明した。 ξ_{1}, ⋯, ξ_{n} がすべて相異なり、 x_{i}^{(k)} が ξ_{i} に十分近いとして、 ε_{i}^{(k)} := x_{i}^{(k)} − ξ_{i} とおけば、

ε_{i}^{(k+1)} = ε_{i}^{(k)} [1 − \frac{∏_{j≠i}(x_{i}^{(k)} − ξ_{j})}{∏_{j≠i}(x_{i}^{(k)} − x_{j}^{(k)})}]

= ε_{i}^{(k)} [1 − ∏_{j≠i} (1 + \frac{ε_{j}^{(k)}}{x_{i}^{(k)} − x_{j}^{(k)}} )]

2 次の項 ε_{j}^{(k)} · ε_{ℓ}^{(k)} を無視すると

ε_{i}^{(k)} ≒ −ε_{i}^{(k)} ∑_{j≠i} \frac{ε_{j}^{(k)}}{x_{i}^{(k)} − x_{j}^{(k)}} .

|x_{i}^{(k)} − x_{j}^{(k)}| ≥ \frac{1}{C} (C は i, j, k によらない)

とすれば

|ε_{i}^{(k+1)}| ≤ C|ε_{i}^{(k)} · ∑_{j≠i} |ε_{j}^{(k)}|

これは 2 乗収束を意味するが、さらに詳しくみると

[figure]

{定理 2.2.1 (DK 法の誤差) \atop | 次の誤差| ≤ C × | 自分自身の誤差| × | 他の誤差の和|}

[figure]
[figure]

定理 2.2.2 (重心の不変性)任意の初期値に対して、 (2.5) によって作られる x_{i}^{(k+1)} (i = 1, 2, ⋯, n) は

(2.6)  ∑_{i=1}^{n} x_{i}^{(k+1)} = −a_{n−1} = ∑_{i=1}^{n} ξ_{i}

を満たす。すなわち、近似根 x_{1}^{(k)}, ⋯, x_{n}^{(k)} の重心は一定で、根 ξ_{1}, ⋯, ξ_{n} の重心に等しい。ゆえに誤差 ε_{i}^{(k+1)} = x_{i}^{(k+1)} − ξ_{i} の和は 0 である。

∑_{i=1}^{n} ε_{i}^{(k+1)} = 0.

[figure]

証明 (そんなに簡単ではない。一松 [8] の §20 や、山本・北川 [17] を参照。 )

2.2.3 復習: 根と係数の関係

この項の内容は代数学の入門書に解説されていることが多い。色々あるが、例えば山崎 [19] をあげておく。

2 次方程式の場合

a_{2}x^{2} + a_{1}x + a_{0} = a_{2}(x − ξ_{1})(x − ξ_{2})

より

−a_{2}(ξ_{1} + ξ_{2}) = a_{1}, a_{2}(ξ_{1}ξ_{2}) = a_{0}.

3 次方程式の場合

a_{3}x^{3} + a_{2}x^{2} + a_{1}x + a_{0} = a_{3}(x − ξ_{1})(x − ξ_{2})(x − ξ_{3})

より

−a_{3}(ξ_{1} + ξ_{2} + ξ_{3}) = a_{2}, a_{3}(ξ_{2}ξ_{3} + ξ_{3}ξ_{1} + ξ_{1}ξ_{2}) = a_{1}, −a_{3}(ξ_{1}ξ_{2}ξ_{3}) = a_{0}.

これを一般化して

a_{n}x^{n} + a_{n−1}x^{n−1} + ⋯ + a_{1}x + a_{0} = a_{n}(x − ξ_{1})(x − ξ_{2}) ⋯ (x − ξ_{n})

について、

(−1)^{s}a_{n} × (ξ_{1}, ⋯, ξ_{n} から s 個取って作った積の和 ) = a_{n−s} (s = 1, 2, ⋯, n).

[figure]

定義 2.2.1 (基本対称式) 不定元 X_{1}, X_{2}, ⋯, X_{n} の整式

S_{s}(X_{1}, X_{2}, ⋯, X_{n}) = ∑_{1≤i_{1}<i_{2}<···<i_{s}≤n} X_{i_{1}}X_{i_{2}} ⋯ X_{i_{s}} (s = 1, 2, ⋯, n)

(i_{1}, ⋯, i_{s} は 1, ⋯, n から取った長さ s の増加列で_{n}C_{s} 個ある) を X_{1}, X_{2}, ⋯, X_{n} の s 次基本対称式と呼ぶ。

[figure]
[figure]

定理 2.2.3 (根と係数の関係)可換体 K 上の n 次の整式

f(X) = ∑_{j=0}^{n} a_{j}X^{j} = a_{n}X^{n} + a_{n−1}X^{n−1} + ⋯ + a_{1}X + a_{0} (a_{n} ≠ 0)

のすべての根 ξ_{1}, ξ_{2}, ⋯, ξ_{n} について

(−1)^{s} \frac{a_{n−s}}{a_{n}} = S_{s}(ξ_{1}, ξ_{2}, ⋯, ξ_{n}). [figure]

[figure]

定義 2.2.2 (対称式) 整域 K における不定元 X_{1}, X_{2}, ⋯, X_{n} の整式で、不定元をどのように入れ換えても等しいものがあったとき、その整式を X_{1}, X_{2}, ⋯, X_{n} の対称式であるという。

[figure]
[figure]

定理 2.2.4 不定元 X_{1}, X_{2}, ⋯, X_{n} の任意の対称式 f は、ある n 不定元の整式 g によって

f = g(S_{1}, S_{2}, ⋯, S_{n})

と表される。

[figure]

その他、判別式、終結式などが勉強しておくべき内容である。

2.2.4 DK 法 — Kerner の解釈

X_{1}, ⋯, X_{n} の s 次基本対称式 (s 個の積全体の和) を S_{s}(X_{1}, ⋯, X_{n}) とおく:

S_{s}(X_{1}, ⋯, X_{n}) = ∑_{i_{1}<i_{2}<···<i_{s}} X_{i_{1}}X_{i_{1}} ⋯ X_{i_{s}} = ∑_{i_{1}<i_{2}<···<i_{s}} ∏_{j=1}^{s} X_{i_{j}}.

根と係数の関係から

(2.7)  f_{s}(ξ_{1}, ⋯, ξ_{n}) := (−1)^{s}S_{s}(ξ_{1}, ⋯, ξ_{n}) − a_{n−s} = 0 (s = 1, 2, ⋯, n).
[figure]

定理 2.2.5 (Kerner 1966) Durand の反復 (2.5) は、 (2.7) を ξ_{1}, ⋯, ξ_{n} についての連立方程式とみなして、それを Newton 法で解く反復計算に他ならない。

[figure]

証明 (書きかけ) 与えられた z_{1}, ⋯, z_{n} に対して

∏_{j=1}^{n} (x − z_{j}) = Q(x) = x^{n} + b_{n−1}x^{n−1} + ⋯ + b_{1}x + b_{0}

とおくと、

\frac{∂Q}{∂z_{j}} = (\frac{∂b_{n−1}}{∂z_{j}} ) x^{n−1} + ⋯ + (\frac{∂b_{1}}{∂z_{j}} ) x + (\frac{∂b_{0}}{∂z_{j}} )
= −(x − z_{1}) ⋯ (x − z_{j−1})(x − z_{j+1}) ⋯ (x − z_{n}).

また根と係数の関係

b_{n−k} = (−1)^{k}S_{k}(z_{1}, ⋯, z_{n})

が成り立つ。ここで x = z_{i} とおくと、 −Q′(z_{i})δ_{ij} となるから、

(−\frac{z_{i}^{n−1}}{Q′(z_{i})}, ⋯, −\frac{z_{i}}{Q′(z_{i})}, −\frac{1}{Q′(z_{i})})

は J(z)^{−1} = (∂b_{n−i}/∂z_{j})^{−1} の第 i 行と一致する。ゆえに J(z)^{−1}f(z) の第 i 成分は

\frac{−z_{i}^{n−1}}{Q′(z_{i})} (b_{n−1} − a_{n−1}) − ⋯ − \frac{1}{Q′(z_{i})}(b_{0} − a_{0}) = −\frac{z_{i}^{n}}{Q′(z_{i})}(b_{n} − a_{n}) − \frac{z_{i}^{n−1}}{Q′(z_{i})}(b_{n−1} − a_{n−1}) − \frac{1}{Q′(z_{i})}(b_{0} − a_{0})
= \frac{−Q(z_{i}) + P (z_{i})}{Q′(z_{i})} = \frac{P (z_{i})}{Q′(z_{i})}.

ゆえに ... の各成分を比較して ... を得る。

[figure]

{系 2.2.1 代数方程式の根がすべて相異なるとき、 DK 法は根に十分近いところで 2 次の \atop 収束となる。}

[figure]

証明 Newton 法の一般的な性質である。

2.2.5 解の精度の検証

次の定理は Durand-Kerner 法に限らず、連立法一般に利用できる便利な定理である。

[figure]

定理 2.2.6 (Smith (1970)) z_{1}, ⋯, z_{n} を相異なる複素数とするとき、

p(z) = z^{n} + a_{n−1}z^{n−1} + ⋯ + a_{1}z + a_{0} = 0

のすべての根は閉円板

Γ_{ℓ} = \{z ∈ C; |z − z_{ℓ}| ≤ r_{ℓ}\}, r_{ℓ} = \frac{n|p(z_{ℓ})|}{∏_{j≠ℓ}(z_{ℓ} − z_{j})} (ℓ = 1, 2, ⋯, n)

の合併に含まれる。その連結成分の一つが m 個の閉円板からなれば、その中にちょうど

m 個の根がある。 [figure]

証明山本 [15] に載っていた定理だが証明は書いていない。杉原・室田 [5] ではステートメントも少し詳しく、証明つきである。

2.2.6 DK 法の長所・短所

DK 法の長所

[figure]

• すべての根を同時に求める反復法。並列演算向き。

• 導関数の計算が不要 (多項式を計算する手続きさえあれば計算できるので便利)。

• 各根の収束の速さがほぼ一様。

• Newton 法と異なり暴走の危険性が少ない。

[figure]

短所

[figure]

• 能率が今一つと言われている。ただし、初期値次第とも言える。

• 代数方程式専用であって、他への応用がない。

• 近接根が求めにくい (一般的な困難と言えるが)。

• 一つの根だけが欲しい場合には無駄が大きい。

[figure]

2.2.7 初期値の取り方の重要性

例 2.2.1 (DK 法が収束しない初期値)p(x) ∈ R[x], x_{i}^{(0)} ∈ R (i = 1, 2, …, n) とすると、

x_{i}^{(k)} ∈ R (i = 1, 2, …, n, k ∈ N)

となるので、虚根には収束し得ないことが分る。実際の数値例については、山本・北川 [17] p.61 を見よ。

2.2.8 Aberth の初期値, DKA 法

以下に示す連立法の初期値は Aberth (1973) により提唱された。アルゴリズム全体を DKA 法と呼ぶ。

p(z) = z^{n} + a_{n−1}z^{n−1} + ⋯ + a_{0} = 0

の根を求める DK 法の初期値 z_{j}^{(0)} (j = 1, 2, ⋯, n) として、

(2.8)  z_{j}^{(0)} = −\frac{a_{n−1}}{n} + r_{0} exp [\sqrt{−1} (\frac{2π(j − 1)}{n} + \frac{π}{2n})] (j = 1, 2, ⋯, n)

ただし r_{0} は閉円板

\overline{B}(−a_{n−1}/n; r_{0}) := \{z ∈ C; |z + a_{n−1}/n| ≤ r_{0}\}

がすべての根を含むように選ぶ。例えば

r_{0} = \frac{|a_{n−1}|}{n} + 1 + max_{0≤i≤n−1} |a_{i}|.

p(z) = 0 の根は行列

( \begin{matrix}0 & & & & −a_{0} \\ 1 & 0 & & & −a_{1} \\ & ... & ... & & ⋮ \\ & & ... & 0 & −a_{n−2} \\ & & & 1 & −a_{n−1} \end{matrix} )

の固有値である。 Gerschgorin の定理から、それらはすべて

|z| ≤ 1 + max_{0≤i≤n−1} |a_{i}|

に含まれる。

また

x_{j}^{(1)} + \frac{a_{1}}{n} = (1 − \frac{1}{n}) (x_{j}^{(0)} + \frac{a_{1}}{n} ) + O ((x_{j}^{(0)} + \frac{a_{1}}{n} )^{−1})

が成り立つので、 r_{0} が十分大きいいとき、 x_{j}^{(k)} (j = 1, 2, ⋯, n) は、円の中心 (根の重心) −a_{1}/n に向かって「直進」する。このことを「最初は 1 次収束」、「終盤は 2 次収束」と言われることもあるが、中盤がどうなるかは分からず、収束証明になっているわけではない。

2.2.9 Ehrlich-Aberth 法

x_{i}^{(k+1)} = x_{i}^{(k)} − \frac{p(x_{i}^{(k)})}{p′(x_{i}^{(k)}) − p(x_{i}^{(k)}) ∑_{j≠i} \frac{1}{x_{i}^{(k)} − x_{j}^{(k)}}}

という漸化式に基づく方法は 3 次収束する。この式の導出は杉原・室田 [5] にある。

2.2.10 根の大きさの限界

[figure]

定理 2.2.7 (根の大きさの評価)

p(z) = ∑_{k=0}^{n} a_{k}z^{k} = a_{n}z^{n} + a_{n−1}z^{n−1} + ⋯ + a_{1}z + a_{0} (a_{n} ≠ 0)

に対して、 K = {k; 0 ≤ k ≤ n − 1, a_{k} ≠ 0} の要素数を κ とし、 κ ≥ 1 と仮定する。このとき、実係数多項式

q(x) = |a_{n}|x^{n} − ∑_{k=0}^{n−1} |a_{k}|x^{k}

は unique な正の実根 r をもち、 p(z) の任意の根 α に対して

|α| ≤ r ≤ r^{∗} := max_{k∈K}(κ|a_{k}/a_{n}|)^{1/(n−k)}

が成り立つ。

[figure]

q(x) = 0 の根 r は、 x = r^{∗} を初期値とする Newton 法により求められる。 (杉原・室田 [5] を見よ。 )

2.2.11 伊理の初期値

(2.9)  z_{j}^{(0)} = −\frac{a_{n−1}}{na_{n}} + r_{0} exp [\sqrt{−1} (\frac{2π(j − 1)}{n} + \frac{3}{2n})] (j = 1, 2, ⋯, n)

伊理 [1] p.152 による。森口 [14] の解説によると「特に 3/(2n) だけ位相をずらせて、実軸や虚軸と特別な関係にならないようにしてあるのが特色である。」とある。

付録A おもちゃ箱 (その他の方法)

A.1 平野法

(この節の内容は全面的に杉原・室田 [5] のつまみ食いである。 )

平野法は Newton 法を少し修正したものであるが、大域的収束性を保証できるという大きな長所がある。 1960 年代後半に

• 平野菅保, 代数方程式の解法および誤差, 第 8 回プログラミングシンポジウム報告集, 情報処理学会, 1967.

• B.Dejon and K.Nickel, A never failing fast convergent root-finding algorithm, Constructive Aspects of the Fundamental Theorems of Algebra (B.Dejon and P.Henrici, eds.), John Wiley, New York, 1969, pp.1–35.

によって独立に提案された算法である。アルゴリズムの解説と収束を保証する定理については、

室田一雄, 平野の変形 Newton 法の大域的収束性, 情報処理学会論文誌, Vol.21 (1980), pp.469–474.

で初めて証明された。

p(z) = ∑_{k=0}^{n} a_{n−k}z^{k}

z^{(ν+1)} = z^{(ν)} + Δz^{(ν)} (ν = 0, 1, ⋯)

Δz^{(ν)} は次のように考えて決める。 p(z) を z^{(ν)} の回りで

(A.1)  p(z^{ν} + ζ) = ∑_{k=0}^{n} c_{n−k}ζ^{k}

と展開する。もちろん

c_{n−k} = \frac{p^{(k)}(z^{(ν)})}{k!} (k = 0, 1, ⋯, n),

特に

c_{n} = p(z^{(ν)}), c_{n−1} = p′(z^{(ν)})

通常の Newton 方では、 z^{(ν)} が根に十分近いことを前提にして、 (A.1) の右辺から k = 0 と k = 1 に対応する項だけを取り出して、

(A.2)  p(z^{(ν)+ζ}) ≒ c_{n−1}ζ + c_{n}

と近似して、これを 0 にする値

ζ = −\frac{c_{n}}{c_{n−1}} = −\frac{p(z^{ν})}{p′(z^{(ν)})}

を修正量としている。平野法では、

p(z^{(ν)} + ζ) ≒ c_{n−k}ζ^{k} + c_{n} (k = 1, 2, ⋯, n)

という 2 項近似を考え、これから定まる

ζ_{k} := (−\frac{c_{n}}{c_{n−k}} )^{1/k} (c_{n−k} = 0 のときは ζ_{k} = ∞ とおく)

を修正量の候補とし、それらの中で絶対値が最小のもの (これを ζ_{m} と書く) を修正量に採用する。

[figure]

補題 A.1.1 (|ζ_{m}| が余裕を持って最小ならば関数値は減少する) 0 < β < 1 とする。

|ζ_{m}| ≤ \frac{β}{1 + β} |ζ_{k}| (1 ≤ k ≤ n, k ≠ m)

ならば

|p(z^{(ν)} + ζ_{m})| ≤ β |p(z^{(ν)})|

である。

[figure]

証明杉原・室田 [5] を見よ。

z^{(ν)} が根 α の十分良い近似値ならば m = 1 となる (すなわち Newton 法と一致)。さらに α が単根であることを仮定すると、補題の条件が成り立つ。

0 < β < 1, λ > 1 は減速の仕方を決めるパラメーターである。例えば β = 3/4, λ = 2 とする。

平野法の第 ν 段の修正量の決定

[figure]

(1) c_{k} (k = 0, 1, 2, ⋯, n) を組み立て除法で計算する。 μ := 1 とおく。

(2) ζ_{k} := (−μc_{n}/c_{n−k})^{1/k} (k = 1, 2, ⋯, n) とおく。

(3) |ζ_{k}| が最小である k (1 ≤ k ≤ n) を m とおく。

(4) |p(z^{(ν)+ζ_{m}})| ≤ (1 − (1 − β)μ)|c_{n}| ならば、 Δz^{(ν)} := ζ_{m} として第 ν 段を終了し、そう

でなければ μ := μ/λ として、 (2) に戻る。

[figure]

(2) における k 乗根の分枝は、 |z^{(ν)} + ζ_{k}| が最も小さくなるなものを選ぶ。具体的には ϕ, ψ_{k} ∈ [0, 1) を

ϕ = \frac{arg z^{(ν)}}{2π} , ψ_{k} = \frac{arg(−c_{n}/c_{n−k})}{2π},

で定めて、 k(1/2 − ϕ) − ψ_{k} に最も近い整数 j_{k} を求め、

(A.3)  ζ_{k} := |μ\frac{c_{n}}{c_{n−k}} |^{1/k} exp [\frac{2πi(ψ_{k} + j_{k})}{k}]

とする。

[figure]

定理 A.1.1 (室田 1980 — 平野法の大域的収束性)任意の初期値 z^{(0)} ∈ C に対して、平野法の近似列 {z^{(ν)}} は、

|p(z^{(ν+1)})| ≤ [1 − (1 − β)θ] |p(z^{(ν)})|

を満たす。ここで、

θ = (\frac{β}{1 + β} )^{2n^{3}} λ^{−n}.

ゆえに {z^{(ν)}} は p(z) のある根に収束する。

[figure]

証明元々は室田 [11] による。杉原・室田 [5] 第 5 章を見よ。

A.2 Horner 法

古くは天元術と言われたとか。一松 [8] にちょこっと説明があるが…。

A.3 Graffe 法

根の自乗を根とする方程式を次々に作る。 URR のように広い範囲を数値を表現できるようなシステムでは有効性が高いかもしれないという指摘がある。

A.4 Bairstow Hitchcock の方法

n 次代数方程式

p_{n}(x) = a_{0}x^{n} + a_{1}x^{n−1} + ⋯ + a_{n}, a_{0} ≠ 0, n ≥ 2

が与えられたとき、 2 次の因数 x^{2} − ux − v を求めることを考える。

p_{n}(x) = (x^{2} − ux − v)(b_{0}x^{n−2} + b_{1}x^{n−3} + ⋯ + b_{n−2}) + b_{n−1}(x − u) + b_{n}

とおいて、係数を比較すると

(A.4)  \{ \begin{matrix}b_{0} = a_{0} \\ b_{1} = a_{1} + ub_{0} \\ b_{k} = a_{k} + ub_{k−1} + vb_{k−2} & (k = 2, 3, ⋯, n) \end{matrix}

が得られる。 b_{k} は (u, v) の関数と考えられるが、我々が求めたいものは b_{n−1} = b_{n} = 0 となるような (u, v) である。そこで u, v に関する連立 1 次方程式

b_{n−1}(u, v) = 0, b_{n}(u, v) = 0

であると考えることができる。これを Newton 法により解くのがベアストウ・ヒッチコックの方法である。詳しい計算手順は省略する。例えば森 [12] に載っている。一松による説明

x^{n} + a_{1}x^{n−1} + ⋯ + a_{n−1}x + a_{n} = 0

の二次の因数 x^{2} − px − q を求めよう。商を x^{n−2} + b_{1}x^{n−3} + ⋯ + b_{n−2}, 剰余を b_{n−1}x + b_{n} とおくと、

b_{1} = a_{1} + p, b_{2} = a_{2} + pb_{1} + q, b_{k} = a_{k} + pb_{k−1} + qb_{k−2} (k = 3, 4, ⋯, n − 1, n)

という漸化式で {b_{k}} が計算できる。次に {c_{k}} を

c_{1} = b_{1} + p, c_{2} = b_{2} + pc_{1} + q, c_{k} = b_{k} + pc_{k−1} + qc_{k−2} (k = 3, 4, ⋯, n − 1, n)

という漸化式で計算すると、

\frac{∂b_{k}}{∂p} = c_{k−1}, \frac{∂b_{k}}{∂q} = c_{k−2}

となることが示される。解くべき方程式は

b_{n}(p, q) = 0, b_{n−1}(p, q) = 0

であるから、

b_{n} + c_{n−1}Δp + c_{n−2}Δq = 0
b_{n−1} + c_{n−2}Δp + c_{n−3}Δq = 0

p_{ν+1} = p_{ν} + \frac{b_{n}c_{n−3} − b_{n−1}c_{n−2}}{(c_{n−2})^{2} − c_{n−1}c_{n−3}} , q_{ν+1} = q_{ν} + \frac{b_{n−1}c_{n−1} − b_{n}c_{n−2}}{(c_{n−2})^{2} − c_{n−1}c_{n−3}}

という反復を行う。初期値には p_{0} = q_{0} = 0 が使われるが、これを Bairstow 法とよぶ。除法を次数の低い方から進める変形を Hitchcock 法と呼ぶ。

一松 [8] から

[figure]

「Bairstow 法は、これまで、実係数の代数方程式の複素根を、実数計算だけで求めるのに便利な方法として、広く使われてきた。実係数なら複素根は必ず共役複素数の対になって現れ、実 2 次因子の根になるからである。しかし機械的にこれを使って失敗した例が多い。その多くは、 Newton 法には大域的収束性が保証されないのに、適当な初期値からいい加減に始めた場合である。特に虚数部が小さい根 α ± iβ があると、あたかも α に実根があるかのように見え、他の実根 γ と合わせた 2 次因子 (x − α)(x − γ) の近似が計算されかねない。ある段階に進むと、 x = α の付近に実根がないため、分母が 0 近くなって、収束しかかった反復列が大きく飛ぶという現象を生ずる。これを防ぐための手法もいろいろ研究されている。筆者のささやかな経験では、根の分布を調べて十分よい初期値から始めることにし、複素数のよい計算体系が利用できる場合には、複素数の 1 変数の Newton 法を使うようにした方が無難と思う。」

[figure]

A.5 Strum の方法

(Fourier 解析で有名な Fourier による先駆的な研究がある ([] を見よ)。行列の固有値問題でよく利用されている。一般の代数方程式の数値解法としては、数値的不安定性が心配になる。 )

Strum 列は数値計算法のあちこちで顔を出す。ここで簡単にまとめておく。なお、以下の記述は森 [12], 一松 [8] から適当に抜き書きしたものである。高木 [6] を読んで必要なところを書いた方がよいか。

固有値問題のノートにも Strum 列のことを書いておいた。時々見比べて、マージすること。

A.5.1 スツルムの定理

実軸上の区間 [a, b] における多項式 f(x) の実根の個数に関する有名な定理。

[figure]

定義 A.5.1 (Strum 列) f(x) ∈ R[x], [a, b] を R の区間とする。このとき次の 4 つの条件を満足する実係数多項式列

f_{0}(x) = f(x), f_{1}(x), f_{2}(x), ⋯, f_{ℓ}(x)

は区間 [a, b] において Strum 列をなすという。

(1) ∀x ∈ [a, b] に対して、隣り合う二つの多項式 f_{k}(x), f_{k+1}(x) は同時に 0 にはならない。

(2) x_{0} ∈ [a, b], k ∈ {1, 2, ⋯, ℓ−1} において f_{k}(x_{0}) = 0 となれば、 f_{k−1}(x_{0})f_{k+1}(x_{0}) < 0.

(3) 列の最後の多項式 f_{ℓ}(x) は [a, b] において一定の符号を持つ。

(4) x_{0} ∈ [a, b] において f(x_{0}) = 0 ならば f′(x_{0})f_{1}(x_{0}) > 0. [figure]

この定理によって、区間 [a, b] 内に存在する f(x) = 0 の解を数を知ることができるが、二分探索の技法と組み合わせることで、解が存在する区間を好きなだけ小さくすることが出来、その区間内の適当な点を解の近似値として採用するという近似解法ができる。これを Strum の方法と呼ぶ。

[figure]

定義 A.5.2 (符号変化の回数)実係数多項式の列 f(x) = f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x) が区間 [a, b] で Strum 列をなすとき、 f(x) の根でない x ∈ [a, b] に対して、関数値の列

f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x)

を左から右に見ていったときの符号の変化の回数 N(x) が定義できる。

• x において、すべての k につき f_{k}(x) ≠ 0 が成り立つ場合の N(x) の定義は明らか。

• ある k について f_{k}(x) = 0 となった場合は、仮定 (Strum 列の条件 (3) と x が f(x) の根でないという仮定) から k = 0, ℓ はありえず 1 ≤ k ≤ ℓ − 1 で、その場合も Strum 列の条件 (2) によって、 f_{k−1}(x)f_{k+1}(x) < 0 となることから、向う三軒両隣では

\begin{matrix}f_{k−1}(x) & f_{k}(x) & f_{k+1}(x) \\ + & 0 & − \end{matrix} または\begin{matrix}f_{k−1}(x) & f_{k}(x) & f_{k+1}(x) \\ − & 0 & + \end{matrix}

のいずれかのパターンしかありえない。そこで、この 3 項で 1 回符号変化したと数えることにする。

[figure]

例 A.5.1 +, +, −, 0, +, +, − では 3 回と数える。

[figure]

定理 A.5.1 (Strum) 実係数多項式の列 f(x) = f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x) は区間 [a, b] で Strum 列をなし、 f(a)f(b) ≠ 0 であるとする。このとき、 f(x) ≠ 0 なる x ∈ [a, b] に対して、関数値の列

f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x)

を左から右に見ていったときの符号の変化の回数を N(x) とすると、方程式 f_{0}(x) = 0 の区間 [a, b] における解の個数は N(a) − N(b) である。

[figure]

証明 f(x) = 0 の解の個数は有限個である。そのうち [a, b] 内にあるものを大きさの順に並べて

x_{1} < x_{2} < ⋯ < x_{n}

とする。 n + 1 個の区間

I_{0} = [a, x_{1}), I_{1} = (x_{1}, x_{2}), I_{2} = (x_{2}, x_{3}), I_{n−1} = (x_{n−1}, x_{n}), I_{n} = (x_{n}, b]

の合併において N(x) は定義できる。以下では

(1) 各区間 I_{j} において N(x) は定数である:

∃{n_{j}}_{j=0}^{n} s.t. N(x) = n_{j} (x ∈ I_{j}, j = 0, 1, ⋯, n)

(2) n_{j+1} − n_{j} = −1 (j = 0, 1, ⋯, n − 1)

であることを証明する (この二つから容易に N(b) − N(a) = n が導かれる)。 (1) の証明 I_{j} 内のある部分区間 J において、すべての多項式 f_{k}(x_{0}) ≠ 0 とすると、 f_{k}(x) の符号は J で一定であるので (符号が変化するには 0 にならないといけない)、 N(x) の値は変化しないことが分かる。

そこで、ある多項式 f_{k}(x) の解 x_{∗} ∈ I_{j} を通過した場合を考える。 Strum 列の条件から k ≠ ℓ, I_{j} の取り方から k ≠ 0 であるから、 1 ≤ k ≤ ℓ − 1 である。 Strum 列の条件 (2) から f_{k−1}(x_{∗})f_{k+1}(x_{∗}) < 0 であるが、連続性から x_{∗} の十分小さな近傍 V を取れば、そこで f_{k−1}(x) と f_{k+1}(x) は定符号となり、 V 上で f_{k−1}(x)f_{k+1}(x) < 0 がなりたつ。それゆえ、列

f_{0}(x), f_{1}(x), ⋯, f_{k−1}(x), f_{k}(x), f_{k+1}(x), ⋯, f_{ℓ}(x)

の部分列

f_{k−1}(x), f_{k}(x), f_{k+1}(x)

の符号については次の二つのいずれか一方だけしか起こらない。

(a) \begin{matrix}f_{k−1}(x) & f_{k}(x) & f_{k+1}(x) \\ V 上 − & & V 上 + \end{matrix} (b) \begin{matrix}f_{k−1}(x) & f_{k}(x) & f_{k+1}(x) \\ V 上 + & & V 上 − \end{matrix}

それゆえ、 (f_{k}(x) の符号が何であっても) 任意の x ∈ V について、部分列 f_{k−1}(x), f_{k}(x), f_{k+1}(x) における符号の変化は 1 として N(x) の値に算入される。これから十分小さな x_{∗} の近傍で N(x) の値に変化はないことが分かる。ゆえに N(x) は I_{j} 上で定数であることが分かる。

(2) の証明 f(x) = 0 の解 x_{j} において考える。 Strum 列の条件 (4) は f′(x_{j}) と f_{1}(x_{j}) が同符号であることを示している。例えば f′(x_{j}) > 0 の場合、十分小さな ε > 0 を取ると、

• f(x) < 0 (x ∈ [x_{j} − ε, x_{j})), f(x_{j}) = 0, f(x) > 0 (x ∈ (x_{j}, x_{j} + ε]).

• f_{1}(x) > 0 (x ∈ (x_{j} − ε, x_{j} + ε)).

これから

n_{j+1} − n_{j} = N(x_{j} + 0) − N(x_{j} − 0) = −1

であることが分かる。 f′(x_{j}) < 0 の場合も同様の議論で n_{j+1} − n_{j} = −1 であることが分かる。

A.5.2 ユークリッドの互除法による Strum 列の生成

多項式 f(x) ∈ R[x] が与えられたとき、 f_{0}(x) = f(x) と f_{1}(x) = f′(x) から Euclid の互除法を行い、関数列 f_{0}(x), f_{1}(x), ⋯, f_{ℓ}(x) を作る:

(A.5)  f_{0}(x) := f(x), f_{1}(x) := f′(x),

f(x) = q_{1}(x)f_{1}(x) − f_{2}(x), deg f_{2}(x) < deg f_{1}(x),
f_{1}(x) = q_{2}(x)f_{2}(x) − f_{3}(x), deg f_{3}(x) < deg f_{2}(x),
f_{2}(x) = q_{3}(x)f_{3}(x) − f_{4}(x), deg f_{4}(x) < deg f_{3}(x),
⋮
(A.6)  f_{k−1}(x) = q_{k−1}(x)f_{k}(x) − f_{k+1}(x), deg f_{k+1}(x) < deg f_{k}(x),
⋮
f_{ℓ−2}(x) = q_{ℓ−1}(x)f_{ℓ−1}(x) − f_{ℓ}(x), deg f_{ℓ}(x) < deg f_{ℓ−1}(x),
f_{ℓ−1}(x) = q_{ℓ}(x)f_{ℓ}(x).

(普通の互除法と異なり、 f_{k−1}(x) を f_{k}(x) で割ったときの普通の剰余 ×(−1) を f_{k+1}(x) とする。 )

よく知られているように f_{ℓ}(x) は f(x) と f′(x) の最大公約多項式であるから、 f(x) ∈ R[x] が重根を持たない場合、 f_{ℓ}(x) ≡ 定数(≠ 0) となることに注意しよう。以下この場合に f(x) = 0 の解 (根) を求めることを考える。 f(x) が重根を持つ場合は f(x) の代わりに g(x) = f(x)/f_{ℓ}(x) を考えることで同様の議論ができる。

[figure]

定理 A.5.2 f(x) ∈ R[x] が重根を持たないとき、 f(a)f(b) ≠ 0 を満たす任意の区間 [a, b] において、 f(x) と f′(x) から互除法で作った剰余列 (詳しくは (A.5), (A.6) 参照)

f(x) = f_{0}(x), f′(x) = f_{1}(x), f_{2}(x), ⋯, f_{ℓ}(x)

は Strum 列をなす。

[figure]

証明まず f_{ℓ}(x) ≡ 定数であるから、 Strum 列の条件 (3) は満たされている。次にある x_{0} ∈ [a, b], ある k ∈ {0, 1, ⋯, ℓ − 1} に対して

f_{k}(x_{0}) = f_{k+1}(x_{0}) = 0

となったとすると、式 (A.6) から f_{k+2}(x_{0}) 以降の f_{j}(x_{0}) もすべて 0 になり、特に f_{ℓ}(x_{0}) = 0. これは f_{ℓ}(x) が定数関数 (≠ 0) であることに矛盾する。ゆえに Strum 列の条件 (1) が満たされる。次にある点 x_{0}, ある k ∈ {1, 2, ⋯, ℓ − 1} に対して f_{k}(x_{0}) = 0 となったとすると、式 (A.6) から

f_{k−1}(x_{0}) = −f_{k+1}(x_{0}).

ゆえに Strum 列の条件 (2) も満たされる (条件 (3) から上式の値は 0 にならないことに注意)。最後に x_{0} が f(x) の根であるとき、 f(x) が重根を持たないという仮定から f′(x_{0}) ≠ 0 で、

f′(x_{0})f_{1}(x_{0}) = f′(x_{0})^{2} > 0

となり、条件 (4) も満たされる。

A.5.3 3 重対角行列の固有多項式と Strum 列

T = ( \begin{matrix}a_{1} & b_{1} & 0 & & 0 \\ b_{1} & a_{2} & b_{2} \\ & & ... & ... & ... \\ & & b_{N−2} & a_{N−1} & b_{N−1} \\ 0 & & 0 & b_{N−1} & a_{N} \end{matrix} )

を実対称三重対角行列とする。

注意 A.5.1 (ちょっとした)実は以下の議論で T の対称性はあまり本質的でなく、対角線をはさんで対称の位置にある成分が同符号、すなわち

t_{k,k−1}t_{k−1,k} > 0 k = 2, ⋯, N

が成り立つことが本質的であるという指摘が一松先生の本にあった。なるほど。行列を Hessenberg 形にする算法で、特に与えられた行列が実対称であった場合には、結果が実対称三重対角になってしまうということで、 (以下 Strum の方法を用いて固有値を求める段階になると) 実対称性はある意味では必要性はあまりないのだが、最初に実対称性がないと三重対角化できないわけで、まあ必ずついてくる「おまけ」ですな。

b_{k} ≠ 0 (k = 1, 2, ⋯, N − 1) と仮定する。もしある k に対して b_{k} ≠ 0 ならば

T = ( \begin{matrix}T′ & O \\ O & T′′ \end{matrix} )

となり、 T′, T′′ の固有値を求める問題に帰着できるから、一般性は失われない。 p_{k}(λ) を λI − T の第 k 主座行列式とする (k = 0, 1, ⋯, N)。すなわち

(A.7)  p_{k}(λ) := \{ \begin{matrix}det(λI_{k} − T_{k}) & (k = 1, 2, ⋯, N) \\ 1 & (k = 0) \end{matrix} .

ただし、

I_{k} = k 次の単位行列, T_{k} = ( \begin{matrix}a_{1} & b_{1} & 0 & & 0 \\ b_{2} & a_{2} & b_{1} \\ & & b_{k−2} & a_{k−1} & b_{k−1} \\ 0 & & 0 & b_{k−1} & a_{k} \end{matrix} ) .

すぐ分かる命題を二つ。

[figure]

\begin{matrix}補題 A.5.1 (漸化式) (A.7) で定義された \{p_{j}(λ)\}_{j=0}^{N} について、以下が成り立つ。 \\ \{ & \begin{aligned}p_{0}(λ) & = & 1, \\ p_{1}(λ) & = & λ − a_{1}, \\ p_{k+1}(λ) & = & (λ − a_{k+1})p_{k}(λ) − b_{k}^{2}p_{k−1}(λ) (k = 1, 2, ⋯, N − 1), \\ p_{N}(λ) & = & det(λI − T). \end{aligned} \end{matrix}

[figure]

証明行列式の行に関する展開定理を用いる。

[figure]

補題 A.5.2 (Strum 列であること)(A.7) で定義された {p_{j}(λ)}_{j=0}^{N} を逆順に並べた多項式列

p_{N}(λ), p_{N−1}(λ), ⋯, p_{1}(λ), p_{0}(λ)

は任意の閉区間 [a, b] において Strum 列である。すなわち

(1) 隣り合う二つの多項式 p_{k}(λ), p_{k+1}(λ) は共通根を持たない。

(2) ある λ_{0} ∈ R において p_{k}(λ_{0}) = 0 ならば p_{k−1}(λ_{0})p_{k+1}(λ_{0}) < 0 (k = 1, 2, ⋯, N −1).

(3) 列の最後の多項式 p_{0}(λ) は R において定符号である。

(4) ある λ_{0} ∈ R において p_{N}(λ_{0}) = 0 ならば p′_{N}(λ_{0})p_{N−1}(λ_{0}) > 0.

[figure]

証明

(1) もしも p_{k}(λ_{0}) = p_{k+1}(λ_{0}) = 0 とすると、漸化式から b_{k}^{2}p_{k−1}(λ_{0}) = 0. b_{k} = 0 と仮定したから p_{k−1}(λ_{0}) = 0. これを繰り返すと

0 = p_{k+1}(λ_{0}) = p_{k}(λ_{0}) = p_{k−1}(λ_{0}) = ⋯ = p_{2}(λ_{0}) = p_{1}(λ_{0}) = p_{0}(λ_{0}).

これから

p_{0}(λ_{0}) = 0

これは p_{0} ≡ 1 に矛盾する。

(2) p_{k}(λ_{0}) = 0 を漸化式に代入すると p_{k+1}(λ_{0}) = −b_{k}^{2}p_{k−1}(λ_{0}). 前項より左辺 ≠ 0. これから p_{k+1}(λ_{0}), p_{k−1}(λ_{0}) は異符号である。

(3) p_{0}(λ) ≡ 1 であるから明らか。

(4) 漸化式

(A.8)  p_{k}(λ) = (λ − a_{k})p_{k−1}(λ) − b_{k−1}^{2}p_{k−2}(λ)

を微分すると、

(A.9)  p′_{k}(λ) = p_{k−1}(λ) + (λ − a_{k})p′_{k−1}(λ) − b_{k−1}^{2}p′_{k−2}(λ).

(A.8) と (A.9) から (A.10)

p′_{k}(λ)p_{k−1}(λ) − p_{k}(λ)p′_{k−1}(λ) = β_{k−1}^{2} (p′_{k−1}(λ)p_{k−2}(λ) − p_{k−1}(λ)p′_{k−2}(λ)) + p_{k−1}^{2}(λ)

が得られる。ここで

q_{k}(λ) := p′_{k}(λ)p′_{k−1}(λ) − p_{k}(λ)p_{k−1}(λ)

とおくと (A.10) は

q_{k}(λ) = p_{k−1}(λ)^{2} + β_{k−1}(λ)^{2}q_{k−1}(λ) (k = 2, 3, ⋯, N).

となる。ところで

q_{1}(λ) = p′_{1}(λ)p_{0}(λ) − p_{1}(λ)p′_{0}(λ) = p′_{1}(λ) = 1 · 1 − (λ − α_{1}) · 0 = 1 > 0

であるから、以下帰納的に

q_{k}(λ) > 0 (k = 2, 3, ⋯, N).

特に

q_{N}(λ) = p′_{N}(λ)p_{N−1}(λ) − p_{N}(λ)p′_{N−1}(λ) > 0

であるが、 p_{N}(λ) = 0 であるから

p′_{N}(λ)p_{N−1}(λ) > 0.

符号の変化数の計算に関する注意 p_{N}(a) ≠ 0 なる a に対して

N(a) ≡^{def.} “{p_{0}(a), p_{1}(a), ⋯, p_{N}(a)}” の符号の変化数

とおく (逆順であっても符号の変化数は同じになる)。例えば

\begin{matrix}p_{0}(a) & p_{1}(a) & p_{2}(a) & p_{3}(a) & p_{4}(a) & p_{5}(a) & p_{6}(a) & p_{7}(a) & p_{8}(a) & p_{9}(a) & p_{10}(a) \\ + & − & − & − & 0 & + & + & + & + & − & − \end{matrix}

では N(a) = 3.

この符号の変化数は (特別な注意をせずに) 数値的に安定して計算できる。つまり絶対値が非常に小さくて、符号の判別がつきにくい場合も、「符号の変化数」そのものは疑いがなく計算できる。例えば

\begin{matrix}p_{k−1}(a) & p_{k}(a) & p_{k+1}(a) \\ + & 絶対値小 & − \end{matrix} または\begin{matrix}p_{k−1}(a) & p_{k}(a) & p_{k+1}(a) \\ − & 絶対値小 & + \end{matrix}

において p_{k}(a) が正であっても負であっても 0 であっても符号の変化数の計算にとっては影響がない。注意すべきは Strum 列の条件 (iv) から

\begin{matrix}p_{k−1}(a) & p_{k}(a) & p_{k+1}(a) \\ + & 絶対値小 & + \end{matrix} または\begin{matrix}p_{k−1}(a) & p_{k}(a) & p_{k+1}(a) \\ − & 絶対値小 & − \end{matrix}

のような場合 (もしこうなったら符号の変化数の計算がむつかしい) が起こり得ないことである。

A.5.4 直交多項式の作る Strum 列

(まあここは単なる覚え書き。後で肉付けするかもしれない。 )

w: [a, b] → R は連続で、有限個の点で 0 になる他は正で、

sup_{k∈N} ∫_{a}^{b} x^{k}w(x) dx < ∞

を満たすような関数とする。このとき [a, b] 上の実数値連続関数全体の集合に

(u, v)_{w} := ∫_{a}^{b} u(x)v(x)w(x) dx

で定義される内積を導入して、内積空間としたものを H_{w}(a, b) とする。

関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とする。

[figure]

命題 A.5.1 H_{w}(a, b) において関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とし、 p_{n}(x) の最高次の係数を μ_{n} とすると、

inf {∥q∥_{w}; q(x) ∈ R[x], deg q(x) = n, q(x) の最高次係数 = 1} = ∥p_{n}/μ_{n}∥_{w}. [figure]

[figure]

命題 A.5.2 H_{w}(a, b) において関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とすると、 p_{n}(x) の根はすべて単根で、区間 [a, b] の内部にある。

[figure]
[figure]

命題 A.5.3 H_{w}(a, b) において関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とすると、次のような 3 項漸化式が存在する。

p_{k}(x) = (α_{k}x + β_{k})p_{k−1}(x) − γ_{k}p_{k−2}(x) (k = 1, 2, ⋯).

ただし p_{−1}(x) ≡ 0 とし、

α_{k} = \frac{μ_{k}}{μ_{k−1}} , β_{k} = −\frac{α_{k}(xp_{k−1}, p_{k−1})_{w}}{λ_{k−1}} ,
γ_{k} = \frac{α_{k}(xp_{k−1}, p_{k−2})_{w}}{λ_{k−2}} = \frac{μ_{k}μ_{k−2}λ_{k−1}}{μ_{k−1}^{2}λ_{k−2}} .

ただし、 λ_{k} = (p_{k}, p_{k})_{w}, μ_{k} = p_{k}(x) の最高次係数。

[figure]
[figure]

命題 A.5.4 H_{w}(a, b) において関数列 {1, x, ⋯, x^{n}} から Gram-Schimidt の直交化法によって得られる直交多項式系を {p_{0}(x), p_{1}(x), ⋯, p_{n}(x)} とすると (最高次の係数が正であるようにしておくと)、

p_{n}(x), p_{n−1}(x), ⋯, p_{1}(x), p_{0}(x)

は区間 [a, b] において Strum 列をなす。

[figure]

A.5.5 一般化された Strum 列

一松先生の本にあった Strum 列の定義の条件は森先生の本よりも少し緩い条件であった (つまり一般化してあると言える)。明示していないけれど p_{k}(λ) の連続性は仮定されている、と思う。

次の性質を持つ関数列 {p_{k}(λ)}_{k=0}^{n} を Strum 系と言う。

(1) 各 p_{k}(λ) の解は有限個; p_{k}(λ_{0}) = p_{k+1}(λ_{0}) = 0 はありえない。

(2) p_{k}(λ_{0}) = 0 のとき p_{k−1}(λ_{0})p_{k+1}(λ_{0}) < 0.

(3) p_{0}(λ) は定符号。

もちろん対応する Strum の定理の結論は少し複雑になる。

[figure]

定理 A.5.3 (一般化された Strum の定理){p_{k}(λ)}_{k=0}^{n} を Strum 系とする。 p_{n}(a)p_{n}(b) ≠ 0 (a < b) のとき、

N(a) − N(b) = ∑_{p_{n}(λ_{0})=0} χ(λ_{0}).

ただし χ(λ_{0}) は、 λ を λ_{0}− から λ_{0}+ に変化させるときの p_{n}(λ)/p_{n−1}(λ) の符号変化を表す量で

χ(λ_{0}) := \{ \begin{matrix}1 & (− から + に変化) \\ 0 & (符号の変化なし) \\ −1 & (+ から − に変化). \end{matrix}
[figure]

証明略

[figure]

系 A.5.1 N(a) ≠ N(B) ならば [a, b] 内に p_{n}(λ) = 0 の解がある。

[figure]
[figure]

系 A.5.2 (実対称三重対角行列の固有値問題)n 次実対称三重対角行列

T = ( \begin{matrix}a_{1} & b_{1} & 0 & & 0 \\ b_{1} & a_{2} & b_{2} \\ & ... & ... & ... \\ & & b_{n−2} & a_{n−1} & b_{n−1} \\ 0 & & 0 & b_{n−1} & a_{n} \end{matrix} ) , b_{k} ≠ 0 (k = 1, 2, ⋯, n − 1)

から p_{k}(λ) = det(λI_{k} − T_{k}) (T_{k} は T の第 k 次首座行列) として作った Strum 列 p_{n}(λ), p_{n−1}(λ), ⋯, p_{1}(λ), p_{0}(λ) ≡ 1 について、 n 個の固有値はすべて実単根であり、各固有値 λ_{0} において、つねに χ(λ_{0}) = 1 である。ゆえに p_{n}(a)p_{n}(b) ≠ 0 となる a < b に対して

N(a) − N(b) = 区間 [a, b] 内の f(x) = 0 の根の個数.

そして p_{k}(λ) の隣接した根の中間に、必ず p_{k−1}(λ) = 0 の根がある (k = 1, 2, ⋯, n)。

[figure]

証明略。

A.6 Bernoulli 法

与えられた解析関数の零点のうちで原点に最も近いものを求める方法と、その多項式への応用 (ベルヌイの方法) を述べる。

[figure]

補題 A.6.1 (ケーニッヒ)

h(z) = c_{0} + c_{1}z + c_{2}z^{2} + ⋯, c_{0} ≠ 0

を D_{R} := |z| < R における有理型関数とする。またこの関数は D_{R} には z = z_{r} にただ一つの単純な極を持つものとする。このとき

\frac{c_{k}}{c_{k+1}} = z_{r} [1 + O (|\frac{z_{r}}{ρ} |^{k+1})] .

ただし ρ は |z_{r}| < ρ < R を満たす数である。

[figure]

複素平面上の円板 D_{R} = {z ∈ C; |z| < R} で定義された有理型関数

f(z) = a_{0} + a_{1}z + a_{2}z^{2} + ⋯, a_{0} ≠ 0

が、 D_{R} 内に唯一の零点 z_{a} を持つとする。また

g(z) = b_{0} + b_{1}z + b_{2}z^{2} + ⋯

を D_{R} で正則で、 g(z_{a}) ≠ 0 を満たす関数とする。このとき

h(z) = \frac{g(z)}{f(z)} = c_{0} + c_{1}z + c_{2}z^{2} + ⋯

は補題の条件を満足する。ゆえに

lim_{k→∞} \frac{c_{k+1}}{c_{k}} = \frac{1}{z_{a}} .

これから f(z) = 0 の絶対値最小の零点が得られる。 g(z) = h(z)f(z) を

b_{0} + b_{1}z + b_{2}z^{2} + ⋯ = (a_{0} + a_{1}z + a_{2}z^{2} + ⋯)(c_{0} + c_{1}z + c_{2}z^{2} + ⋯)

から

(A.11)  \begin{matrix}a_{0}c_{0} = b_{0} \\ a_{0}c_{1} + a_{1}c_{0} = b_{1} \\ a_{0}c_{2} + a_{1}c_{1} + a_{2}c_{0} = b_{2} \\ ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ & ⋯ \end{matrix}

ベルヌイ法上述の方法を

p_{n}(z) = a_{0}z^{n} + a_{1}z^{n−1} + ⋯ + a_{n−1}z + a_{n} = 0
f(z) := z^{n}p_{n}(z^{−1})

とおくと f(z) も n 次多項式である。そして f(z) の絶対値最小の零点 z_{a} と p_{n}(z) の絶対値最大の零点 z_{p} との間には

z_{p} = \frac{1}{z_{a}}

の関係がある。 (この項書きかけ) 一松先生曰く、

(代数方程式の) 最大根を求めるには古くから (Daniel) Bernoulli 法が知られている (本質的には (固有値問題を解くための) 累乗法と同じ)。

A.7 商差法

商差法 (quotient-difference algorithm, QD 法) は、昔から知られていたが、 1950 年代初めの Rutishauser が導入し、後に Henrici が発展させた方法である。

一松 [8] (1982) に解説がある。「現在でも低次 (10 次以下) で実根の多い代数方程式の、序盤ないし中盤むきの手法としては有用であり、多くの関連話題とともに、理論的には重要である。しかし、汎用とするには難点が多く、一つの補助算法と考えた方がよい。 (中略) 以上のような理由で、現在では実用価値が低下しているので、算法の解説のみにとどめ、定理の証明を省略した。」

A.8 次数低下

Durand-Kerner 法のように代数方程式 p(x) = 0 の全根を同時に探索する方法を採用した場合には必要ないが、例えば Newton 法のような方法で、一つの根の近似根 ξ を求めた後、 p(x) を (x − ξ) で割って、 1 次低い代数方程式に還元することを次数低下あるいは減次 (deflation) と呼ぶ。

次数低下をする際に、除法を「いつでも最高次の項から行う」のは危ない (一松 [8] の §21 を見よ)。

この点については、櫻井 [7] に詳しい。

n 次多項式 f(x) を x − α で割った商を q(x), 剰余を r とおくと、

f(x) = q(x)(x − α) + r

であり、 f(α) = r であるから、 α が f(x) の根のときは r = 0 となる。

また Q(x), R を

f(x) = Q(x)(x − α) + Rx^{n}

を満たすように決めた場合にも f(α) = 0 のときには R = 0 であり、従って q(x) = Q(x) であるが、数値計算で求める場合、丸め誤差の観点からは二つのやり方に差がある。前者を高次項からの減次、後者を低次項からの減次という。

高次項からの減次 f(x) = a_{0} + a_{1}x + ⋯ + a_{n}x^{n} を因子 x − α によって高次項から減次した多項式 q(x) = q_{0} + q_{1}x + ⋯ + q_{n−1}x^{n−1} を求める。

1. q_{n−1} := a_{n} とおく。

2. i = n − 1, n − 2, ⋯, 1 の順に、以下を繰り返す。 q_{i−1} := q_{i}α + a_{i}.

3. q(x) := q_{0} + q_{1}x + ⋯ + q_{n−1}x^{n−1} とおく。

低次項からの減次 f(x) = a_{0} + a_{1}x + ⋯ + a_{n}x^{n} を因子 x − α によって高次項から減次した多項式 Q(x) = Q_{0} + Q_{1}x + ⋯ + Q_{n−1}x^{n−1} を求める。

1. Q_{0} := −a_{0}/α とおく。

2. i = 1, 2, ⋯, n − 1 の順に、以下を繰り返す。 Q_{i} := (Q_{i−1} − a_{i})/α.

3. Q(x) := Q_{0} + Q_{1}x + ⋯ + Q_{n−1}x^{n−1} とおく。

絶対値の小さな根については高次項からの減次、絶対値の大きな根については低次項からの減次をした方が誤差の影響が少なくなる。絶対値が中くらいの場合、 |a_{k}α^{k}| が最大となる k = k_{max} を探し、

q̃(x) = q_{0} + q_{1}x + ⋯ + q_{k_{max}}x^{k_{max}} + Q_{k_{max}+1}x^{k_{max}+1} + ⋯ + Q_{n−1}x^{n−1}.

を採用すると良い結果が得られることがある (そうだ)。

A.9 Jenkins and Traub 法

Jenkins–Traub 法については、 Jenkins and Traub [20], Ralston and Rabinowitz [21]

A.10 その他

この節の記述は櫻井 [7] による。

A.10.1 Euler 法

f を x_{0} の近傍で 2 次式で近似すると

f(x) ≒ f(x_{0}) + f′(x_{0})(x − x_{0}) + \frac{f′′(x_{0})}{2!}(x − x_{0})^{2}

そこで

f(x_{0}) + f′(x_{0})(x − x_{0}) + \frac{f′′(x_{0})}{2!}(x − x_{0})^{2} = 0

の 2 根のうち x_{0} に近い方を新たな近似解とする。

A.10.2 ハレー法 A.10.3 ポメンタール法

付録B 応用解析IVの数値実験

2003 年度の応用解析 IV の内容から。

本日は、 C++ で書かれたサンプル・プログラムを応用解析 IV ホームページから入手して、 DK 法の簡単な数値実験を行なう。

B.1 6701 号室のワークステーションの利用

数学科の学生は 6701 号室にあるワークステーション、パソコンを利用することができる。ユーザー名は情報科学センターのワークステーションと (とりあえず) 同じユーザー名となるが、独立に管理している (だから、例えばパスワードは一応別物^{1}となる)。

6701 号室の利用については、以下の WWW ページを読むこと。

                     http://www.math.meiji.ac.jp/howto/

B.2 応用解析 IV ホームページ

この講義のホームページ

           http://www.math.meiji.ac.jp/~mk/lecture/ouyoukaiseki4/

が用意されていて、サンプル・プログラムはそこから入手できる。

なお、講義ノートなどは、受講者のみに公開している (一般には公開していない)。ユーザー名は ouyoukaiseki4 で、パスワードはである。

サンプル・プログラムを手元のコンピューターにコピーする方法 Netscape を使う場合、サンプル・プログラムをアクセスしてから、 File メニューの「名前を付けて保存 (S)」を選択する。 IE を使う場合は…

1自分で同じものを選ぶことは可能である。

B.3 C++ について

B.3.1 まずはサンプル・プログラム

 quadratic-equation.C

[figure]

 1  // quadratic-equation.C --- 複素係数の 2 次方程式を解く。
 2  //

3 // コンパイルは (6701 号室では例えば)

 4  //      g++ -o quadratic-equation quadratic-equation.C
 5  //

6 // 注意: 素朴な (桁落ちの対策などしていない) アルゴリズムを使っている。

 7
 8  #include <iostream.h>
 9  #include <complex.h>
10
11  int main()
12  {
13    complex<double> a, b, c, D, x1, x2;
14
15    cout << "複素係数の 2 次方程式 a x^2+b x+c=0 (a ≠ 0) を解きます。 " << endl;

16 cout << " 複素数は ( ) でくくり、実部と虚部をカンマ , で区切って表す。 "

17         << endl;
18    cout << " 例えば 1+2i は (1,2) と表わします。 " << endl;
19    cout << "入力してください。 " << endl;
20
21    cout << "a="; cin >> a;
22    cout << "b="; cin >> b;
23    cout << "c="; cin >> c;
24
25    cout << "a=" << a << ", b=" << b << ", c=" << c << endl;
26
27    D = sqrt(b * b - 4 * a * c);
28    x1 = (-b + D) / (2 * a);
29    x2 = (-b - D) / (2 * a);
30    cout << "x1=" << x1 << endl;
31    cout << "x2=" << x2 << endl;
32
33    return 0;
34  }

[figure]

実行例

[figure]

oyabun% g++ -o quadratic-equation quadratic-equation.C
oyabun% ./quadratic-equation
複素係数の 2 次方程式 a x^2+b x+c=0 (a ≠ 0) を解きます。

複素数は ( ) でくくり、実部と虚部をカンマ , で区切って表す。例えば 1+2i は (1,2) と表わします。入力してください。

a=1
b=1
c=1
a=(1,0), b=(1,0), c=(1,0)
x1=(-0.5,0.866025)
x2=(-0.5,-0.866025)
oyabun%

[figure]

B.3.2 C++ の (とても簡潔な) 紹介

C++ は C 言語を基礎にして作られたプログラミング言語である。

大抵の C プログラムは C++ プログラムでもある

という言葉が示すように、 C++ は C 言語を拡張したものと考えることができる。主な拡張点は、オブジェクト指向プログラミング^{2}を実現しやすくする機能である。

今回は、 C 言語よりも複素数を扱うプログラミングがしやすいという理由から、 C++ を採用した (上のプログラムを C 言語で書くとどうなるかを想像してみよう)。

B.3.3 サンプル・プログラムを読むための C++ の知識

入出力に C 言語の関数 (printf() 等) を使うことも可能だが、普通は >>, << 演算子を使うことが推奨されている。これを使う場合、プログラムの先頭付近で\underline{#include <iostream.h>} としておく必要がある。

標準出力への出力 (通常は画面への表示)

[figure]

  cout << 式; あるいは cout << 式_{1} << 式_{2} << · · · << 式_{n};

(なお行末を表わす名前 endl (END of Line?) を覚えておこう。 )

[figure]

標準入力からの入力 (通常はキーボードからの入力)

[figure]

  cin >> 変数名;

[figure]

複素数データの処理

[figure]

実験に用いる C++ 処理系では、 complex<double> というクラス^{a}が用意されている。

    • プログラムの先頭部分で\underline{#include <complex.h>} とする。

• \underline{complex<double> 変数名;} で複素数を表現する変数が定義できる。

• complex<double> クラスのデータに対して、通常の演算子 +, -, *, / で四則演算が

      出来る。比較演算子 ==, != もある。また real(), imag(), abs(), sqrt() などの関

数も用意されている (意味は多分名前で分かってもらえると思う)。

• (もちろん) <<, >> 演算子を用いての入出力が可能である。その際、 x + \sqrt{−1}y (x, y ∈ R) は (x,y) のように丸括弧 (,) を用いて 2 次元ベクトルのように表現する。例えば虚数単位\sqrt{−1} は (0,1) となる。

aクラスというのはオブジェクト指向言語における概念であるが、ここでは説明しない。型を拡張したも

のぐらいに思っておけばサンプル・プログラムは理解可能である。 [figure]

2これが何であるかは説明しない。

B.3.4 g++ — 実験に用いる C++ 処理系

6701 号室のワークステーションでは、 GCC というプログラミング言語処理系があり、 C, C++, Objective C, Fortran 77 などのプログラミング言語で書かれたプログラムを扱うことができる。

• GCC で C++ プログラムをコンパイルするには、通常 g++ というコマンドを使う。オプション等の使い方は C コンパイラーである gcc に準じる。例えば、サンプル・プロ

グラム quadratic-equation.C のコンパイル＆実行には、

[figure]

   oyabun% \underline{g++ quadratic-equation.C}
   oyabun% \underline{./a.out}

[figure]

または

[figure]

   oyabun%  \underline{g++ -o quadratic-equation quadratic-equation.C}
   oyabun% \underline{./quadratic-equation}

[figure]

のようにすれば良い。

• GCC では、 C++ で書かれたプログラムのファイル名末尾は .C (大文字の C であることに注意) とする習慣である。

B.4 グラフィックス・ライブラリィ GLSC について

サンプル・プログラム DK3.c では、収束の様子を図示するために、 C, C++ や Fortran プログラムから利用可能なグラフィックス・ライブラリィ GLSC^{3} を利用している。

GLSC について詳しいことが知りたければ、

     http://www.math.meiji.ac.jp/~mk/labo/howto/index.html#GLSC

を参照すること。

3龍谷大学で開発された。 ftp://ftp.st.ryukoku.ac.jp/pub/ryukoku/software/math/ で公開されている。

サンプル・プログラムの中で GLSC を使っているところ

[figure]

// GLSC のヘッダーファイルをインクルード

extern "C" {
#define G_DOUBLE
#include <glsc.h>
}
.......

// グラフィックス・ライブラリィ GLSC の初期化

g_init("DKGRAPH", 220.0, 220.0);
g_device(G_BOTH);

// 座標系の指定

g_def_scale(0,
            real(g) - 1.1 * R, real(g) + 1.1 * R,
            imag(g) - 1.1 * R, imag(g) + 1.1 * R,
            10.0, 10.0, 200.0, 200.0);
g_sel_scale(0);
.......

// 線種、マーカーの指定

g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
g_sel_line(0);
g_def_marker(0, G_RED, 2, G_MARKER_CIRC);
g_sel_marker(0);
.......

// 初期値と初期値を置く円を描く

g_circle(real(g), imag(g), r0, G_YES, G_NO);
for (i = 0; i < n; i++)
  g_marker(real(x[i]), imag(x[i]));
.......

// 図示する

    g_move(real(x[i]), imag(x[i]));
    g_plot(real(newx[i]), imag(newx[i]));
    g_marker(real(newx[i]), imag(newx[i]));
.......
// GLSC 終了の処理

cout << "終了にはグラフィックスのウィンドウをクリックして下さい。 " <<endl;

g_sleep(-1.0);
g_term();

[figure]

上のプログラム断片を理解するための注

• g は根の重心の座標

• R は {z ∈ C; |z − g| ≤ R} がすべての根を含むことが保証できる正数 (講義で解説済み)

• x[], newx[] はそれぞれ {x_{i}^{(k)}}_{i=1}^{n}, {x_{i}^{(k+1)}}_{i=1}^{n} を記憶する変数

B.5 素朴な Newton 法の実行例

 complex-newton.C

[figure]

 1  /*
 2   * complex-newton.C -- Newton 法で方程式 f(x)=0 を解く
 3   * コンパイル: g++ -o complex-newton complex-newton.C
 4   * 実行: ./complex-newton
 5   */
 6
 7  #include <iostream.h>
 8  #include <complex.h>
 9
10  int main()
11  {
12      int i, maxitr = 100;
13      complex<double> f(complex<double>), dfdz(complex<double>), x, dx;
14      double eps;
15
16      cout << " 初期値 x0, 許容精度ε=";
17      cin >> x >> eps;
18
19      cout.precision(16);
20      for (i = 0; i < maxitr; i++) {
21          dx = - f(x) / dfdz(x);
22          x += dx;
23          cout << "f(" << x << ")=" << f(x) << endl;
24          if (abs(dx) <= eps) break;
25      }
26      return 0;
27  }
28
29  /*

30 * この関数の例は杉原・室田『数値計算法の数理』岩波書店, p.67 による

31   *   f(z) = z^3-2z+2
32   *   1 実根 (≒-1.769292354238631),2 虚根 (≒ 0.8846461771193157 ± 0.5897428050222054)

33 * を持つが、原点の近くに初期値を取ると、 0 と 1 の間を振動する。

34   */
35
36  complex<double> f(complex<double> z)
37  {
38    /* return z * z * z - 2 * z + 2; */
39    return z * (z * z - 2) + 2;
40  }
41

42 /* 関数 f の導関数 (df/dx のつもりで名前をつけた) */

43  complex<double> dfdz(complex<double> z)
44  {
45    return 3 * z * z - 2;
46  }

[figure]

 complex-newton 実行結果

[figure]

oyabun% g++ -o complex-newton complex-newton.C
oyabun% ./complex-newton
 初期値 x0, 許容精度ε=(0.07,0.07) 1e-15
f((1.000479890500462,0.01402105439035678))=(0.9998905285720316,0.01405867910352982)
f((0.008690071144917044,0.08327942929892811))=(1.982439704951537,-0.1671175729050391)
f((0.9899874424723547,0.002680506223014176))=(0.9902658531083521,0.002520260799472426)
f((-0.06300187281112724,0.01783194417259117))=(2.125813775994601,-0.03545722093219378)
(中略)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
f((0.9999999999999999,0))=(1,0)
f((-9.992007221626409e-16,0))=(2.000000000000002,0)
oyabun%

[figure]

B.6 DK 法のサンプル・プログラム

B.6.1 ソースプログラム DK3.C

1  // DK3.C -- DKA 法で代数方程式の根を求める
2  //

3 // 数学科 6701 号室のワークステーション環境でのコンパイル 4 // (グラフィックス・ライブラリィ GLSC を使っているので 5 // 情報科学センターではコンパイルできない。悪しからず。 )

6  //     g++ -o DK3 DK3.C -I/usr/local/include -lglscd -lX11 -lsocket

7 // あるいは

8  //     ccmg+ DK3.C
9  //

10 // Win32 環境下の glscwin でも (若干の修正の下に) 動作させられることが 11 // 分かったので、画面に表示するメッセージを変更しました (2000/10/18)。

12
13  #include <iostream.h>
14  #include <iomanip.h>   // setprecision() のため
15  #include <math.h>
16  #include <complex.h>   // complex<double> のため
17

18 // GLSC のヘッダーファイルをインクルード

19  extern "C" {
20  #define G_DOUBLE
21  #include <glsc.h>
22  }

24 // プロトタイプ宣言

25  complex<double> polynomial(int, complex<double> *, complex<double>);
26  complex<double> bunbo(int, complex<double> *, int);
27

28 // 解きたい代数方程式の次数の最高値

29  #define MAXN 100
30
31  int main()
32  {
33    int i, k, n;
34    complex<double> a[MAXN+1], x[MAXN], newx[MAXN], dx[MAXN];
35    complex<double> g, I(0,1);
36    double r0, R, max, pi;
37

38 // 数学定数 (円周率) の準備

39    pi = 4 * atan(1.0);
40

41 // 表\示の桁数を 16 桁にする

42    cout << setprecision(16);
43

44 // 方程式の入力

45    cout << "次数 n を入力してください (1 ≦ n ≦" << MAXN << "): ";
46    cin >> n;
47    if (n > MAXN || n <= 0) {
48      cerr << "次数は" << MAXN << "以下の自然数として下さい。 " << endl;
49      exit(0);
50    }
51    for (i = 0; i <= n; i++) {
52      cout << (n-i) << "次の係数を入力してください: ";
53      cin >> a[i];
54    }
55

56 // 多項式を最高次の係数で割って monic にする

57    cout << "monic にします。 " << endl;
58    for (i = 1; i <= n; i++)
59      a[i] /= a[0];
60    a[0] = 1;
61    cout << "修正した係数" << endl;
62    for (i = 0; i <= n; i++)
63      cout << "a[" << i << "]=" << a[i] << endl;
64

65 // Aberth の初期値を配置する円の決定

66    g = - a[1] / n;
67    cout << "根の重心" << g << endl;
68    max = 0;
69    for (i = 1; i <= n; i++)
70      if (abs(a[i]) > max)
71        max = abs(a[i]);
72    cout << "max|a_i|=" << max << endl;
73    r0 = abs(g) + 1 + max;
74    cout << "根は重心中心で、半径 r0=" << r0 << "の円盤内にある" << endl;
75

76 cout << "円の半径 (分からなければ上の値を指定してください): ";

77    R = r0;
78    cin >> r0;
79    if (r0 > R)
80      R = r0;

81 cout << "図は根の重心を中心として半径 " << R

 82  #ifdef __CYGWIN__

83 << "の円が表\示できるようにします。 " << endl;

 84  #else

85 << "の円が表示できるようにします。 " << endl;

 86  #endif
 87
 88    // Aberth の初期値
 89    cout << "初期値" << endl;
 90    for (i = 0; i < n; i++) {
 91      double theta;
 92      theta = 2 * i * pi / n + pi / (2 * n);
 93      x[i] = g + r0 * exp(I * theta);
 94      cout << x[i] << endl;
 95    }
 96

97 // グラフィックス・ライブラリィ GLSC の初期化

 98    g_init("DKGRAPH", 120.0, 120.0);
 99    g_device(G_BOTH);

100 // 座標系の指定

101    g_def_scale(0,
102        real(g) - 1.1 * R, real(g) + 1.1 * R,
103        imag(g) - 1.1 * R, imag(g) + 1.1 * R,
104        10.0, 10.0, 100.0, 100.0);
105    g_sel_scale(0);

106 // 線種、マーカーの指定

107    g_def_line(0, G_BLACK, 0, G_LINE_SOLID);
108    g_sel_line(0);
109    g_def_marker(0, G_RED, 2, G_MARKER_CIRC);
110    g_sel_marker(0);
111

112 // 初期値と初期値を置く円を描く

113    g_circle(real(g), imag(g), r0, G_YES, G_NO);
114    for (i = 0; i < n; i++)
115      g_marker(real(x[i]), imag(x[i]));
116
117    // 反復
118    for (k = 1; k <= 1000; k++) {
119      double error;
120      cout << "第" << k << "反復" << endl;
121      for (i = 0; i < n; i++) {
122        dx[i] = polynomial(n, a, x[i]) / bunbo(n, x, i);
123        newx[i] = x[i] - dx[i];
124        // 図示する
125        g_move(real(x[i]), imag(x[i]));
126        g_plot(real(newx[i]), imag(newx[i]));
127        g_marker(real(newx[i]), imag(newx[i]));
128      }
129      // 更新
130      for (i = 0; i < n; i++) {
131        x[i] = newx[i];
132        cout << x[i] << endl;
133      }

134 // 変化量を計算する (ここは非常に素朴)

135      error = 0.0;
136      for (i = 0; i < n; i++)
137        error += abs(dx[i]);
138      cout << "変化量=" << error << endl;

139 // 変化量が小さければ収束したと判断して反復を終了する。

140      if (error < 1e-12)
141        break;
142    }
143
144    // GLSC 終了の処理
145    cout << "終了には" << endl;

146 cout << " グラフィックスのウィンドウをクリック (X の場合)。 " << endl;

147    cout << " グラフィックスのウィンドウをクローズ (Windows の場合)。 " << endl;
148    g_sleep(-1.0);
149    g_term();
150
151    return 0;
152  }
153
154  // 多項式の値の計算 (Horner 法)
155  complex<double> polynomial(int n,
156     complex<double> *a,
157     complex<double> z)
158  {
159    int i;
160    complex<double> w;
161    w = a[0];
162    for (i = 1; i <= n; i++)
163      w = (w * z + a[i]);
164    return w;
165  }
166
167  // Π (ｚ i-ｚ j) の計算
168  // j ≠ i
169  complex<double> bunbo(int n,
170        complex<double> *z,
171        int i)
172  {
173    int j;
174    complex<double> w(1,0);
175    for (j = 0; j < n; j++)
176      if (j != i)
177        w *= (z[i] - z[j]);
178    return w;
179  }

B.6.2 コンパイル＆実行例

z^{3} − 1 = 0 を解いてみた例を示す。

[figure]

oyabun% ccmg+ DK3.C
oyabun% ./DK3

次数 n を入力してください (1 ≦ n ≦ 100): 3 3 次の係数を入力してください: 1 2 次の係数を入力してください: 0 1 次の係数を入力してください: 0 0 次の係数を入力してください: -1

monic にします。

修正した係数

a[0]=(1,0)
a[1]=(0,0)
a[2]=(0,0)
a[3]=(-1,0)
根の重心 (-0,-0)
max|a_i|=1

根は重心中心で、半径 r0=2 の円盤内にある円の半径 (分からなければ上の値を指定してください): 2 図は根の重心を中心として半径 2 の円が表示できるようにします。初期値

(1.732050807568877,0.9999999999999999)
(-1.732050807568877,1.000000000000001)
(-3.673940397442059e-16,-2)

第 1 反復

(1.196367205045918,0.5944978830179634)
(-1.113033871712585,0.7388354503153702)
(-0.08333333333333359,-1.333333333333333)
変化量=2.015564437074637

第 2 反復

(0.9103764382644928,0.2474704563982415)
(-0.6695039210592528,0.6646738943447257)
(-0.2408725172052399,-0.912144350742967)
変化量=1.34906218596283

(途中省略)

第 6 反復

(0.9999999987110025,4.085881563506899e-10)
(-0.4999999997093489,0.86602540246384)
(-0.4999999990016535,-0.8660254028724281)
変化量=0.0001103168208393921

第 7 反復

(1,-1.053338202674773e-18)
(-0.5,0.8660254037844386)
(-0.5,-0.8660254037844386)
変化量=4.056615492443548e-09

第 8 反復

(1,0)
(-0.5,0.8660254037844387)
(-0.5,-0.8660254037844387)
変化量=1.665556150858263e-16

終了にはグラフィックスのウィンドウをクリック (X の場合)。グラフィックスのウィンドウをクローズ (Windows の場合)。

oyabun%

[figure]

付録C 情報処理IIから

C.1 レポート課題

課題 1

色々な方程式を二分法、 Newton 法で解いてみよ。初期値の選び方を変えて、収束するかしないか、試しなさい。収束の判定条件には注意を払うこと (特に理由なく低い精度の答を出しただけの場合は減点)。最終的に得られる精度や、必要な反復の回数はどの程度になるか。 Newton 法の収束のための十分条件を本などで探すことができたら、それも説明すること。

ここでは解説しないが、初等関数なども計算機の中では、四則演算などの簡単な演算の組合せで計算されていることが多い。

課題 2

平方根\sqrt{a} や立法根 a^{1/3} を方程式の解の形に定式化して、 Newton 法で解いて見よ。その結果を C 言語のライブラリィ関数 sqrt() や pow() ^{1} で計算した値と比較してみよ。

特に逆関数の計算にも使える (y が与えられているとして方程式 f(x) = y を x について解けば、 x = f^{−1}(y) が得られるはず) ことに注意しよう。

課題 3

C 言語のライブラリィ関数 asin(), acos(), atan() は使わずに、 arcsin, arccos, arctan を計算するプログラムを作り、実験せよ。

課題 4

連立方程式

\begin{aligned}x^{2} − y^{2} + x + 1 & = & 0 \\ 2xy + y & = & 0 \end{aligned}

を Newton 法を用いて解くプログラムを作れ。ヒント:

x⃗ = ( \begin{matrix}x \\ y \end{matrix} ) ,

1pow(a,b) で a の b 乗が計算できる。例えば pow(2.0, 1.0/3.0) で 2 の立法根が計算できる。

f(x⃗) = ( \begin{matrix}x^{2} − y^{2} + x + 1 \\ 2xy + y \end{matrix} )

とおくと、方程式は f(x⃗) = 0 と書ける。 f の x⃗ における Jacobi 行列を f′(x⃗) とすると、 Newton 法の式は

x⃗_{n+1} = x⃗_{n} − [f′(x⃗_{n})]^{−1}f(x⃗_{n})

となる。初期値 x⃗_{0} を ( \begin{matrix}1 \\ 1 \end{matrix} ), ( \begin{matrix}1 \\ −1 \end{matrix} ) として実験せよ。

C.2 例題を解くプログラム例

二分法のプログラム bisection.c^{2}, newton.c^{3} は情報処理 II の WWW ページ http://www.

math.meiji.ac.jp/~mk/syori2/ から入手できる。 WWW ブラウザー「ファイル・メニュー」

の「名前をつけて保存」を用いてセーブする。 2003 年度情報科学センターの教育用情報処理室の Windows 環境からは、 z: ドライブ (UNIX 環境のホームディレクトリィ) にセーブするのが便利であろう。

C.2.1 Newton 法の場合

1 を初期値として、許容精度 10^{−15} を指示して Newton 法で解かせたのが以下の結果 (入力

は\underline{1 1e-15})。

[figure]

  isc-xas06% gcc -o newton newton.c -lm
  isc-xas06% ./newton
   初期値 x0, 許容精度ε=1 1e-15
  f(   7.503638678402439e-01)=-1.89e-02
  f(   7.391128909113617e-01)=-4.65e-05
  f(   7.390851333852840e-01)=-2.85e-10
  f(   7.390851332151607e-01)= 0.00e+00
  f(   7.390851332151607e-01)= 0.00e+00
  isc-xas06%

[figure]

注意 C.2.1 Newton 法の繰り返しを停止させるための良いルールを独力で発見するのはかなり難しい。上の例題プログラムの採用したルールは、多くのプログラムで採用されているやり方ではあるが、いつでもうまく行く方法とは言えない。現時点で「これが良い方法」と見なされている方法は一応あるが、結構複雑なのでここでは説明しない。

C.2.2 二分法の場合

区間 (0, 1) 内に解があることがわかるから、二分法で許容精度 10^{−15} を指示して解かせたのが以下の結果 (入力は\underline{0 1 1e-15})。関数値が、区間の左端では正、右端では負になったまま区間が縮小して行くのを理解しよう。

2http://www.math.meiji.ac.jp/~mk/syori2/cprog/bisection.c

3http://www.math.meiji.ac.jp/~mk/syori2/cprog/newton.c

[figure]

isc-xas06% gcc -o bisection bisection.c -lm
isc-xas06% ./bisection

探す区間の左端α, 右端β, 許容精度ε=0 1 1e-15

f(   5.000000000000000e-01)= 3.78e-01, f(   1.000000000000000e+00)=-4.60e-01
f(   5.000000000000000e-01)= 3.78e-01, f(   7.500000000000000e-01)=-1.83e-02
f(   6.250000000000000e-01)= 1.86e-01, f(   7.500000000000000e-01)=-1.83e-02
f(   6.875000000000000e-01)= 8.53e-02, f(   7.500000000000000e-01)=-1.83e-02
f(   7.187500000000000e-01)= 3.39e-02, f(   7.500000000000000e-01)=-1.83e-02
f(   7.343750000000000e-01)= 7.87e-03, f(   7.500000000000000e-01)=-1.83e-02
f(   7.343750000000000e-01)= 7.87e-03, f(   7.421875000000000e-01)=-5.20e-03
f(   7.382812500000000e-01)= 1.35e-03, f(   7.421875000000000e-01)=-5.20e-03

中略

f(   7.390851332151556e-01)= 8.55e-15, f(   7.390851332151627e-01)=-3.44e-15
f(   7.390851332151591e-01)= 2.55e-15, f(   7.390851332151627e-01)=-3.44e-15
f(   7.390851332151591e-01)= 2.55e-15, f(   7.390851332151609e-01)=-4.44e-16
f(   7.390851332151600e-01)= 1.11e-15, f(   7.390851332151609e-01)=-4.44e-16
f(   7.390851332151600e-01)= 1.11e-15
isc-xas06%

[figure]

付録D 多次元の Newton 法の例

非線形の微分方程式を離散化して解く際に Newton 法を使う例。

D.1 ターゲット問題

次の非線形 2 点境界値問題の解を差分法で求める。

(D.1)  −u′′(x) = u(x)^{2} (x ∈ (0, 1)),
(D.2)  u(0) = u(1) = 0.

D.2 復習

対応する線形問題は 1 次元 Poisson 方程式の Dirichlet 境界値問題である。

(D.3)  −u′′(x) = f(x) (x ∈ (0, 1))
(D.4)  u(0) = u(1) = 0

h = 1/N, x_{i} = ih, u_{i} = u(x_{i}) (i = 0, 1, 2, ⋯, N) とおき、 u_{i} の近似値 U_{i} を差分法で求めることにする。

A = \frac{1}{h^{2}} ( \begin{matrix}2 & −1 \\ −1 & 2 & −1 \\ & ... & ... & ... \\ & & −1 & 2 & −1 \\ & & & −1 & 2 \end{matrix} ) , U⃗ = ( U_{1} ; U_{2} ; ⋮ ; U_{N−1} ) , f⃗ = ( f_{1} ; f_{2} ; ⋮ ; f_{N−1} )

とおくと、差分方程式は

AU⃗ = f⃗

になる。

D.3 ターゲット問題の差分近似

これは次の方程式で U⃗ を定めるのが自然だろう。

(D.5)  AU⃗ = ( (U_{1})^{2} ; (U_{2})^{2} ; ⋮ ; (U_{N−1})^{2} ) .

(D.5) は非線形方程式である

D.4 Newton 法

F⃗ (U⃗) := AU⃗ − ( (U_{1})^{2} ; (U_{2})^{2} ; ⋮ ; (U_{N−1})^{2} )

とおくと、 (D.5) は

(D.6)  F⃗ (U⃗) = 0

となる。 Newton 法とは、初期値 U⃗^{(0)} を適当に選んで、後は漸化式

(D.7)  U⃗^{(k+1)} = U⃗^{(k)} − (F⃗′(U^{(k)}))^{−1} F⃗ (U⃗^{(k)})

でベクトル列 {U⃗^{(k)}}_{k∈N} を定めるというものである。

F⃗′(U⃗) = A − 2 ( \begin{matrix}U_{1} \\ & U_{2} \\ & & . & . & . \\ & & & & & U_{N−1} \end{matrix} )

であるから、ヤコビ行列 F⃗′(U⃗) は三重対角行列である。 (D.7) において逆行列が現れるが、逆行列を計算せずに、連立 1 次方程式を解く形で計算を

遂行すべきことに注意しよう。

D.5 サンプル・プログラム

/*
 * Newton.c

* 非線形 2 点境界値問題

 *    -u''=u^2 in (0,1), u(0)=u(1)=0

* を差分法で離散化して得られる非線形方程式を Newton 法で解く。

*/
#include <stdio.h>
#include <math.h>
#include <matrix.h>
#include "fplot.h"
#include "trid-lu.h"
void mul_mv(int n,
          vector ab,
          vector al, vector ad, vector au,
          vector b)
{
int i, nm1 = n - 1;
ab[0]   =                  ad[0] * b[0] + au[0] * b[1];
for (i = 1; i < nm1; i++)
  ab[i] = al[i] * b[i-1] + ad[i] * b[i] + au[i] * b[i+1];
ab[nm1]   = al[nm1] * b[nm1-1] + ad[nm1] * b[nm1];
}
double norm(int n, vector x)
{
return sqrt(dotprod(n, x, x));
}
int main()
{
int N, i, k;
vector al,ad,au,akl,akd,aku;
vector U, x;
double h, h2, du, H;
N = 100;
h = 1.0 / N;
h2 = h * h;
al = new_vector(N+1); ad = new_vector(N+1); au = new_vector(N+1);
akl = new_vector(N+1); akd = new_vector(N+1); aku = new_vector(N+1);
U = new_vector(N+1);
x = new_vector(N+1);
/* 初期値 */
printf("H (10 位で OK)="); scanf("%lg", &H);
for (i = 0; i <= N; i++)
  U[i] = H;
U[0] = U[N] = 0.0;
/* A */
for (i = 1; i < N; i++) {
  al[i] = - 1.0 / h2; ad[i] = 2.0 / h2; au[i] = - 1.0 / h2;
}
openpl(); fspace2(-0.2, -2.0, 1.2, 20.0);
fmove(0.3, 15.0);
label("-u''=u^2 in (0,1), u(0)=u(1)=0");
linemod("dotted");
fline(-0.2, 0.0, 1.2, 0.0); fline(0.0, -2.0, 0.0, 20.0);
linemod("solid");
for (k = 1; k < 100; k++) {

/* A U^k の計算 */
mul_mv(N - 1, x+1, al+1, ad+1, au+1, U+1);
/* F(U^k) の計算 */
for (i = 1; i < N; i++)
  x[i] -= U[i] * U[i];
/* F'(U^k) の計算 */
for (i = 1; i < N; i++) {
  akl[i] = al[i]; akd[i] = ad[i] - 2 * U[i]; aku[i] = au[i];
}
/* F'(U^k)^{-1} U(U^k) の計算 */
trid(N-1, akl+1, akd+1, aku+1, x+1);
/* */
du = norm(N-1, x+1);
printf("du=%g\n", du);
/* U^{k+1} の計算 */
for (i = 1; i < N; i++)
  U[i] -= x[i];
/* */
#ifdef NONE
for (i = 0; i <= N; i++)
  printf("U[%d]=%g\n", i, U[i]);
#endif
fmove(0.0, U[0]);
for (i = 1; i <= N; i++)
  fcont(i * h, U[i]);
if (du < 1.0e-12)
  break;
}
{
double min, max;
max = U[0]; min = U[0];
for (i = 1; i <= N; i++) {
  if (U[i] > max)
    max = U[i];
  else if (U[i] < min)
    min = U[i];
}
printf("min=%g, max=%g\n", min, max);
}
mkplot("Newton.plot");
closepl();
return 0;
}

61

関連図書

[1] 伊理^{いり}正夫^{まさお}, 数値計算, 朝倉書店 (1981).

[2] 伊理正夫, 藤野和建, 数値計算の常識, 共立出版 (1985).

[3] 齋藤友克^{ともかつ}, 竹島卓^{たく}, 平野照比古, グレブナー基底の計算実践篇 Risa/Asir で解く, 東京大学出版会 (2003).

[4] 篠原能材^{よしたね}, 数値解析の基礎, 日新出版 (初版 1978, 5 版 1997).

[5] 杉原正顯^{まさあき}, 室田^{むろた}一雄, 数値計算法の数理, 岩波書店 (1994).

[6] 高木貞治, 代数学講義改訂新版, 共立出版 (1965).

[7] 名取亮^{まこと}編, 長谷川秀彦^{ひでひこ}, 櫻井^{さくらい}鉄也, 桧山^{ひやま}澄子^{すみこ}, 周^{しゅう}偉東^{いとう}, 花田孝郎^{たかお}, 北川高嗣^{たかし}, 数値計算法第 5 章「非線形方程式」 (櫻井鉄也執筆), オーム社 (1998).

[8] 一松^{ひとつまつ}信^{しん}, 数値解析、朝倉書店 (1982).

[9] 一松信, 数学とコンピュータ, 共立出版 (1995).

[10] 増田久弥, 非線型数学, 朝倉書店 (1985).

[11] 室田一雄, 平野の変形 Newton 法の大域的収束性, 情報処理学会論文誌, Vol.21 (1980), pp.469–474.

[12] 森正武^{まさたけ}, 数値解析, 共立出版 (1973).

[13] 森口繁一, 数値計算術, 共立出版 (19??).

[14] 森口繁一, 数値計算工学, 岩波書店 (1989).

[15] 山本哲朗^{てつろう}, 数値解析入門 [新訂版], サイエンス社 (2003).

(1976 年初版発行の定番本の待望の改訂版)

[16] 山本哲朗, Newton 法とその周辺, 数学, Vol.37 (1985), pp.1–15.

[17] 山本哲朗, 北川高嗣^{たかし}, 数値解析演習, サイエンス社 (1991).

[18] 渡部^{わたなべ}力^{つとむ}, 名取亮^{まこと}, 小国^{おぐに}力^{つとむ}監修, Fortran77 による数値計算ソフトウェア第 11 章「代数方程式」 (小野寺嘉孝執筆), 丸善 (1989).

[19] 山崎圭次郎, 基礎代数, 岩波講座応用数学, 岩波書店 (1994).

[20] M.A.Jenkins and J.F.Traub, Commun. ACM, 15 (1972), pp.97–99.

[21] A.Ralston and P.Rabinowitz, A First Course in Numerical Analysis, McGraw-Hill, second edition (1978).

[22] L.B.Rall, Computational Solution of Nonlinear Operator Equations, John Wiley, New York (1969).

[23] L.Schwartz, 解析学 I, 東京図書 (19??).

[24] E.Zeidler, Nonlinear Functional Analysis and Its Applications I, Springer–Verlag Berlin, 1986.