コンテンツにスキップ
Rikai

定常反復法：Jacobi 法・Gauss–Seidel 法・SOR 法の収束と最適な緩和係数

数学数値線形代数 第9章修士無料公開約 38 分2026-10-06

対象読者: 本テーマの行列ノルム・条件数の章と、ゲルシュゴリンの定理・ペロン・フロベニウスの定理などの固有値評価を学んだ人。反復法の知識は仮定しません。

前提となる記事:行列のノルムと線形方程式の条件数：誘導ノルム・後退誤差・特異行列までの距離、行列の固有値評価と不等式：ゲルシュゴリンからペロン・フロベニウス、カントロヴィッチまで

生 Markdown

0. この記事の要点

Section titled “0. この記事の要点”

大きな連立方程式を、一度に解くのではなく「いまの答えの見積もりを少しずつ直す」操作の繰り返しで解く方法を扱います。たとえば 2 つの未知数 x1x_1, x2x_2 があるとき、x2x_2 の今の値から x1x_1 を計算し直し、その x1x_1 から x2x_2 を計算し直す、という往復です。この記事を読むと、そうした繰り返しがいつ正しい答えに近づくのか、近づく速さは何で決まるのか、更新の幅をどう選ぶと速くなるのかが分かります。

  • 連立一次方程式 Ax=bA\boldsymbol{x} = \boldsymbol{b} の係数行列を A=M−NA = M - N と分け、Mxk+1=Nxk+bM\boldsymbol{x}_{k+1} = N\boldsymbol{x}_k + \boldsymbol{b} を繰り返すのが定常反復法です。任意の初期値で収束することと、反復行列のスペクトル半径（固有値の絶対値の最大）が ρ(M−1N)<1\rho(M^{-1}N) < 1 であることは同値で、1 回あたり誤差はおよそ ρ\rho 倍になります。
  • 狭義対角優位（各行で対角成分が他の成分の絶対値の和より大きいこと）なら Jacobi 法と Gauss–Seidel 法は収束します。正定値 Hermite 行列なら、SOR 法は 0<ω<20 < \omega < 2 で収束します（Ostrowski–Reich）。逆に、どんな行列でも ρ(Lω)≥∣ω−1∣\rho(\mathcal{L}_\omega) \ge |\omega - 1| なので、0<ω<20 < \omega < 2 は収束に必要です（Kahan）。
  • 矛盾なく順序づけられた行列（未知数の番号づけが、隣り合う未知数の間で規則正しく上下する行列。§4 で定義します）、たとえば三重対角行列や、格子の自然な順序・赤黒順序（格子点を市松模様に塗り分け、赤の点を先に並べる順序）の離散 Laplace 行列では、Jacobi 法の固有値 μ\mu と SOR 法の固有値 λ\lambda が (λ+ω−1)2=λω2μ2(\lambda + \omega - 1)^2 = \lambda\omega^2\mu^2 で結ばれます（Young の定理）。ここから ρGS=ρJ2\rho_{\mathrm{GS}} = \rho_{\mathrm{J}}^2 と最適な緩和係数 ωopt=2/(1+1−ρJ2)\omega_{\mathrm{opt}} = 2/(1 + \sqrt{1 - \rho_{\mathrm{J}}^2}) が出ます。
  • 離散 Laplace 方程式（格子幅 hh）では ρJ=cos⁡πh\rho_{\mathrm{J}} = \cos\pi h で、10 桁縮めるのに必要な反復回数は Jacobi 法・Gauss–Seidel 法で h−2h^{-2}、最適 SOR 法で h−1h^{-1} に比例します。
  • この遅さが、Krylov 部分空間法（共役勾配法）と前処理、多重格子法へ進む理由です。

1. 動機

Section titled “1. 動機”

偏微分方程式を差分法で離散化すると、未知数が 10610^6 を超える疎な連立一次方程式が現れます。Gauss 消去は fill-in（もとは 00 だった成分が消去の途中で非零になること）のためにメモリも演算量も膨らみますが、行列とベクトルの積や「各未知数を隣接する値から更新する」操作なら、非零成分の個数に比例する手間で済みます。

Laplace 方程式 −Δu=f-\Delta u = f の 5 点差分では、各格子点の値は「上下左右の値の平均 + h2f/4+\ h^2f/4」です。この関係を、いまの近似値を右辺に入れて左辺を更新する、という形で繰り返すのが Jacobi 法です。更新した値をすぐ使えば Gauss–Seidel 法、更新幅を ω\omega 倍に伸ばせば SOR 法（逐次過緩和法）になります。

本章の問いは 3 つです。いつ収束するか。どれくらいの速さか。ω\omega をどう選べばよいか。答えはすべて反復行列のスペクトル半径で書けます。スペクトル半径と行列の冪の関係は 行列のノルムと線形方程式の条件数 の結果を使い、固有値の計算は 冪乗法・逆反復の章 と同じく固有ベクトルでの展開が基本になります。

2. 分割と定常反復

Section titled “2. 分割と定常反復”

定義 2.1（行列の分割と定常反復）↗

A∈Cn×nA \in \mathbb{C}^{n\times n} を正則とする。A=M−NA = M - N（MM は正則）を AA の分割という。x0∈Cn\boldsymbol{x}_0 \in \mathbb{C}^n から

Mxk+1=Nxk+b,すなわちxk+1=Gxk+M−1b,G=M−1NM\boldsymbol{x}_{k+1} = N\boldsymbol{x}_k + \boldsymbol{b},\qquad\text{すなわち}\qquad \boldsymbol{x}_{k+1} = G\boldsymbol{x}_k + M^{-1}\boldsymbol{b},\quad G = M^{-1}N

で xk\boldsymbol{x}_k を作る反復を定常反復、GG を反復行列という。

「定常」とは、GG と M−1bM^{-1}\boldsymbol{b} が kk によらないことを指します。MM は「My=cM\boldsymbol{y} = \boldsymbol{c} が安く解ける」ように選びます。AA の対角部分を DD、狭義下三角部分の符号を変えたものを LL、狭義上三角部分の符号を変えたものを UU として A=D−L−UA = D - L - U と書くと、代表的な 3 つの方法は次の分割です（DD は正則、つまり対角成分はすべて 00 でないとする）。

方法MM反復行列
Jacobi 法DDGJ=D−1(L+U)G_{\mathrm{J}} = D^{-1}(L + U)
Gauss–Seidel 法D−LD - LGGS=(D−L)−1UG_{\mathrm{GS}} = (D - L)^{-1}U
SOR 法（緩和係数 ω≠0\omega \ne 0）ω−1D−L\omega^{-1}D - LLω=(D−ωL)−1((1−ω)D+ωU)\mathcal{L}_\omega = (D - \omega L)^{-1}\bigl((1 - \omega)D + \omega U\bigr)

Gauss–Seidel 法は ω=1\omega = 1 の SOR 法です。成分で書くと、Gauss–Seidel 法は i=1,2,…,ni = 1, 2, \ldots, n の順に

xi(k+1)=1aii(bi−∑j<iaijxj(k+1)−∑j>iaijxj(k))x_i^{(k+1)} = \frac{1}{a_{ii}}\Bigl(b_i - \sum_{j<i}a_{ij}x_j^{(k+1)} - \sum_{j>i}a_{ij}x_j^{(k)}\Bigr)

と、更新済みの成分をすぐに使います。SOR 法はこの右辺を x~i\tilde{x}_i として xi(k+1)=(1−ω)xi(k)+ωx~ix_i^{(k+1)} = (1 - \omega)x_i^{(k)} + \omega\tilde{x}_i とします。D−LD - L と D−ωLD - \omega L は対角成分が 00 でない下三角行列で、行列式は対角成分の積（命題 3.3［行列式とその性質］ の 2）なので正則です。

小さな例で動きを見ます。A=(2−1−12)A = \begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix}, b=(1,1)T\boldsymbol{b} = (1, 1)^{\mathsf{T}} の解は x∗=(1,1)T\boldsymbol{x}^{*} = (1, 1)^{\mathsf{T}} です。x0=0\boldsymbol{x}_0 = \boldsymbol{0} から Jacobi 法を 3 回回すと、x1=x2x_1 = x_2 が 0.5, 0.75, 0.8750.5,\ 0.75,\ 0.875 と進み、誤差は毎回ちょうど半分になります。Gauss–Seidel 法では、更新したばかりの x1x_1 をすぐ x2x_2 の計算に使うので、(x1,x2)(x_1, x_2) は (0.5,0.75)(0.5, 0.75), (0.875,0.9375)(0.875, 0.9375), (0.96875,0.984375)(0.96875, 0.984375) と進み、誤差は毎回 1/41/4 になります。同じ 3 回でも、誤差は Jacobi 法で 0.1250.125、Gauss–Seidel 法で 0.031250.03125 以下です。この「毎回何倍になるか」が、次の定理が扱う反復行列のスペクトル半径で、この例では Jacobi 法が 1/21/2、Gauss–Seidel 法が 1/41/4 です（演習で確かめます）。

次の定理は、反復が収束するかどうかが反復行列 GG のスペクトル半径だけで決まることを述べます。3 の ε\varepsilon は、スペクトル半径がほぼ「1 回あたりの縮み率」になることを正確に言うための小さな余裕です。

定理 2.2（定常反復の収束の同値条件）✓

定義 2.1 の定常反復について、x∗=A−1b\boldsymbol{x}^{*} = A^{-1}\boldsymbol{b} とする。

  1. 誤差 ek=xk−x∗\boldsymbol{e}_k = \boldsymbol{x}_k - \boldsymbol{x}^{*} は ek=Gke0\boldsymbol{e}_k = G^k\boldsymbol{e}_0 を満たす。
  2. すべての初期値 x0∈Cn\boldsymbol{x}_0 \in \mathbb{C}^n について xk→x∗\boldsymbol{x}_k \to \boldsymbol{x}^{*} となるための必要十分条件は ρ(G)<1\rho(G) < 1 である。AA, MM, b\boldsymbol{b} が実のとき、「すべての実の初期値」としても同じである。
  3. ρ(G)<1\rho(G) < 1 のとき、任意の ε>0\varepsilon > 0 に対し、e0\boldsymbol{e}_0 によらない定数 cεc_\varepsilon があって ∥ek∥2≤cε(ρ(G)+ε)k∥e0∥2\|\boldsymbol{e}_k\|_2 \le c_\varepsilon(\rho(G) + \varepsilon)^k\|\boldsymbol{e}_0\|_2。一方、GG の固有値 λ\lambda で ∣λ∣=ρ(G)|\lambda| = \rho(G) となるものの固有ベクトルを e0\boldsymbol{e}_0 にとると ∥ek∥2=ρ(G)k∥e0∥2\|\boldsymbol{e}_k\|_2 = \rho(G)^k\|\boldsymbol{e}_0\|_2 である。
証明（定理 2.2）

  1. Mx∗=Nx∗+bM\boldsymbol{x}^{*} = N\boldsymbol{x}^{*} + \boldsymbol{b}（Ax∗=bA\boldsymbol{x}^{*} = \boldsymbol{b} の書き直し）を反復の式から引くと Mek+1=NekM\boldsymbol{e}_{k+1} = N\boldsymbol{e}_k、すなわち ek+1=Gek\boldsymbol{e}_{k+1} = G\boldsymbol{e}_k です。

  2. ρ(G)<1\rho(G) < 1 なら Gk→OG^k \to O（系 5.3［行列のノルムと線形方程式の条件数］）なので、1 から任意の e0\boldsymbol{e}_0 で ek→0\boldsymbol{e}_k \to \boldsymbol{0} です。逆に ρ(G)≥1\rho(G) \ge 1 なら、∣λ∣≥1|\lambda| \ge 1 の固有値 λ\lambda と固有ベクトル z\boldsymbol{z} をとり x0=x∗+z\boldsymbol{x}_0 = \boldsymbol{x}^{*} + \boldsymbol{z} とすると ek=λkz\boldsymbol{e}_k = \lambda^k\boldsymbol{z} で、∥ek∥=∣λ∣k∥z∥↛0\|\boldsymbol{e}_k\| = |\lambda|^k\|\boldsymbol{z}\| \not\to 0 です。実の場合、すべての実の初期値で収束するなら、z=zR+izI\boldsymbol{z} = \boldsymbol{z}_{\mathrm{R}} + i\boldsymbol{z}_{\mathrm{I}}（実部と虚部）について GkzR→0G^k\boldsymbol{z}_{\mathrm{R}} \to \boldsymbol{0}、GkzI→0G^k\boldsymbol{z}_{\mathrm{I}} \to \boldsymbol{0}（GG は実行列）なので Gkz→0G^k\boldsymbol{z} \to \boldsymbol{0} となり、上と矛盾します。

  3. スペクトル半径と誘導ノルムの定理（定理 5.2［行列のノルムと線形方程式の条件数］ の 2）により、あるノルム ∥⋅∥ε\|\cdot\|_\varepsilon とその誘導ノルムで ∥G∥ε≤ρ(G)+ε\|G\|_\varepsilon \le \rho(G) + \varepsilon です。よって ∥ek∥ε≤(ρ(G)+ε)k∥e0∥ε\|\boldsymbol{e}_k\|_\varepsilon \le (\rho(G) + \varepsilon)^k\|\boldsymbol{e}_0\|_\varepsilon で、有限次元のノルムの同値性（補題 3.2［行列のノルムと線形方程式の条件数］）から a∥⋅∥2≤∥⋅∥ε≤b∥⋅∥2a\|\cdot\|_2 \le \|\cdot\|_\varepsilon \le b\|\cdot\|_2 となる正の定数 aa, bb があるので、cε=b/ac_\varepsilon = b/a とすればよい。後半は 2 の後半と同じ計算です。

3 の上下の評価から、誤差を 10−d10^{-d} 倍にするのに必要な反復回数は、最悪の初期値についておよそ dln⁡10/(−ln⁡ρ(G))d\ln 10/(-\ln\rho(G)) です。この分母 −ln⁡ρ(G)-\ln\rho(G) を漸近収束率といいます。ρ(G)=1−η\rho(G) = 1 - \eta（η\eta は小さい）なら −ln⁡ρ(G)≈η-\ln\rho(G) \approx \eta で、回数は 1/η1/\eta に比例します。

3. 対角優位な行列と正定値な行列

Section titled “3. 対角優位な行列と正定値な行列”

定理 3.1（狭義対角優位なら Jacobi 法と Gauss–Seidel 法は収束する）✓

A=(aij)A = (a_{ij}) が行について狭義対角優位、すなわちすべての ii で ∣aii∣>∑j≠i∣aij∣|a_{ii}| > \sum_{j\ne i}|a_{ij}| とする。このとき ρ(GGS)≤∥GJ∥∞<1\rho(G_{\mathrm{GS}}) \le \|G_{\mathrm{J}}\|_\infty < 1 かつ ρ(GJ)≤∥GJ∥∞<1\rho(G_{\mathrm{J}}) \le \|G_{\mathrm{J}}\|_\infty < 1 であり、Jacobi 法と Gauss–Seidel 法は任意の初期値で収束する。

証明（定理 3.1）

AA は狭義対角優位なので正則で（系 2.3［行列の固有値評価と不等式］ の 1）、対角成分は 00 でありません。si−=∑j<i∣aij∣/∣aii∣s_i^- = \sum_{j<i}|a_{ij}|/|a_{ii}|、si+=∑j>i∣aij∣/∣aii∣s_i^+ = \sum_{j>i}|a_{ij}|/|a_{ii}| とおくと、仮定は si−+si+<1s_i^- + s_i^+ < 1 です。

Jacobi 法. GJG_{\mathrm{J}} の第 ii 行は −aij/aii-a_{ij}/a_{ii}（j≠ij \ne i）と対角の 00 なので、行和の最大の公式（命題 4.2［行列のノルムと線形方程式の条件数］ の 2）から ∥GJ∥∞=max⁡i(si−+si+)<1\|G_{\mathrm{J}}\|_\infty = \max_i(s_i^- + s_i^+) < 1 です。スペクトル半径は誘導ノルム以下（定理 5.2［行列のノルムと線形方程式の条件数］ の 1）なので ρ(GJ)<1\rho(G_{\mathrm{J}}) < 1 です。

Gauss–Seidel 法. λ\lambda を GGSG_{\mathrm{GS}} の固有値、z\boldsymbol{z} を固有ベクトルとし、∣zi∣=max⁡j∣zj∣>0|z_i| = \max_j|z_j| > 0 となる ii をとります。Uz=λ(D−L)zU\boldsymbol{z} = \lambda(D - L)\boldsymbol{z} の第 ii 成分は

λaiizi=−λ∑j<iaijzj−∑j>iaijzj\lambda a_{ii}z_i = -\lambda\sum_{j<i}a_{ij}z_j - \sum_{j>i}a_{ij}z_j

です。両辺の絶対値をとって ∣aii∣∣zi∣|a_{ii}||z_i| で割ると、∣zj∣≤∣zi∣|z_j| \le |z_i| から ∣λ∣≤∣λ∣si−+si+|\lambda| \le |\lambda|s_i^- + s_i^+、すなわち ∣λ∣(1−si−)≤si+|\lambda|(1 - s_i^-) \le s_i^+ です。1−si−>si+≥01 - s_i^- > s_i^+ \ge 0 なので ∣λ∣≤si+/(1−si−)<1|\lambda| \le s_i^+/(1 - s_i^-) < 1 で、さらに si+/(1−si−)≤si−+si+s_i^+/(1 - s_i^-) \le s_i^- + s_i^+（両辺に 1−si−1 - s_i^- を掛けると si+≤si−+si+−si−(si−+si+)s_i^+ \le s_i^- + s_i^+ - s_i^-(s_i^- + s_i^+)、つまり si−(si−+si+)≤si−s_i^-(s_i^- + s_i^+) \le s_i^- で、これは si−+si+<1s_i^- + s_i^+ < 1 から成り立つ）なので ∣λ∣≤∥GJ∥∞|\lambda| \le \|G_{\mathrm{J}}\|_\infty です。収束は 定理 2.2 から従います。

正定値な行列では、ノルムとしてエネルギーノルム ∥e∥A=(e∗Ae)1/2\|\boldsymbol{e}\|_A = (\boldsymbol{e}^{*}A\boldsymbol{e})^{1/2} を使うと見通しがよくなります。AA が正定値なので e∗Ae\boldsymbol{e}^{*}A\boldsymbol{e} は e≠0\boldsymbol{e} \ne \boldsymbol{0} で正であり、これを長さの 2 乗とみなしたものです。次の補題は、誤差がこの長さで見て毎回確実に縮むための分割の条件を与えます。証明では ∥e∥A2−∥Ge∥A2\|\boldsymbol{e}\|_A^2 - \|G\boldsymbol{e}\|_A^2 を計算し、それがちょうど M+M∗−AM + M^{*} - A の形の二次形式になることから、この条件が出てきます。

補題 3.2（エネルギーノルムが減る分割）✓

AA を正定値 Hermite 行列、A=M−NA = M - N を分割とし、M+M∗−AM + M^{*} - A が正定値 Hermite 行列であるとする。このとき G=M−1N=I−M−1AG = M^{-1}N = I - M^{-1}A について、任意の e≠0\boldsymbol{e} \ne \boldsymbol{0} で ∥Ge∥A<∥e∥A\|G\boldsymbol{e}\|_A < \|\boldsymbol{e}\|_A であり、ρ(G)<1\rho(G) < 1 である。

証明（補題 3.2）

y=M−1Ae\boldsymbol{y} = M^{-1}A\boldsymbol{e} とおくと Ge=e−yG\boldsymbol{e} = \boldsymbol{e} - \boldsymbol{y}、Ae=MyA\boldsymbol{e} = M\boldsymbol{y} です。

∥Ge∥A2=e∗Ae−y∗Ae−e∗Ay+y∗Ay\|G\boldsymbol{e}\|_A^2 = \boldsymbol{e}^{*}A\boldsymbol{e} - \boldsymbol{y}^{*}A\boldsymbol{e} - \boldsymbol{e}^{*}A\boldsymbol{y} + \boldsymbol{y}^{*}A\boldsymbol{y}

で、y∗Ae=y∗My\boldsymbol{y}^{*}A\boldsymbol{e} = \boldsymbol{y}^{*}M\boldsymbol{y}、e∗Ay=(Ae)∗y=y∗M∗y\boldsymbol{e}^{*}A\boldsymbol{y} = (A\boldsymbol{e})^{*}\boldsymbol{y} = \boldsymbol{y}^{*}M^{*}\boldsymbol{y} なので

∥e∥A2−∥Ge∥A2=y∗(M+M∗−A)y\|\boldsymbol{e}\|_A^2 - \|G\boldsymbol{e}\|_A^2 = \boldsymbol{y}^{*}(M + M^{*} - A)\boldsymbol{y}

です。AA, MM は正則なので e≠0\boldsymbol{e} \ne \boldsymbol{0} なら y≠0\boldsymbol{y} \ne \boldsymbol{0} で、右辺は正です。固有値については、Ge=λeG\boldsymbol{e} = \lambda\boldsymbol{e}（e≠0\boldsymbol{e} \ne \boldsymbol{0}）なら ∥Ge∥A=∣λ∣∥e∥A\|G\boldsymbol{e}\|_A = |\lambda|\|\boldsymbol{e}\|_A で、∥e∥A>0\|\boldsymbol{e}\|_A > 0 なので ∣λ∣<1|\lambda| < 1 です。

定理 3.3（Ostrowski–Reich の定理）✓

AA を正定値 Hermite 行列とする。SOR 法は 0<ω<20 < \omega < 2 のとき任意の初期値で収束し、反復ごとに誤差のエネルギーノルムは真に減少する。

証明（定理 3.3）

AA は Hermite なので U=L∗U = L^{*} で、A=D−L−L∗A = D - L - L^{*} です。DD の対角成分 aii=ei∗Aeia_{ii} = \boldsymbol{e}_i^{*}A\boldsymbol{e}_i は正定値性から正です。SOR 法の M=ω−1D−LM = \omega^{-1}D - L について

M+M∗−A=2ωD−L−L∗−(D−L−L∗)=(2ω−1)DM + M^{*} - A = \frac{2}{\omega}D - L - L^{*} - (D - L - L^{*}) = \Bigl(\frac{2}{\omega} - 1\Bigr)D

で、0<ω<20 < \omega < 2 なら 2/ω−1>02/\omega - 1 > 0 なので正定値です。補題 3.2 と 定理 2.2 から主張が従います。

系 3.4（正定値行列の Gauss–Seidel 法）✓

正定値 Hermite 行列に対する Gauss–Seidel 法は任意の初期値で収束する。

証明（系 3.4）

Gauss–Seidel 法は ω=1\omega = 1 の SOR 法なので、定理 3.3 の特別な場合です。

Jacobi 法には同じ主張がありません。M=DM = D なので条件は「2D−A2D - A が正定値」で、これは正定値な AA でも成り立つとは限りません（演習を参照）。

定理 3.5（Kahan の必要条件）✓

AA の対角成分がすべて 00 でないとき、任意の ω≠0\omega \ne 0 について ρ(Lω)≥∣ω−1∣\rho(\mathcal{L}_\omega) \ge |\omega - 1| である。とくに SOR 法が任意の初期値で収束するなら 0<ω<20 < \omega < 2（ω\omega が実数のとき）である。

証明（定理 3.5）

三角行列の行列式は対角成分の積（命題 3.3［行列式とその性質］ の 2）なので det⁡(D−ωL)=det⁡D\det(D - \omega L) = \det D、det⁡((1−ω)D+ωU)=(1−ω)ndet⁡D\det((1 - \omega)D + \omega U) = (1 - \omega)^n\det D です。積の定理（定理 6.1［行列式とその性質］）から

det⁡Lω=det⁡(D−ωL)−1det⁡((1−ω)D+ωU)=(1−ω)n.\det\mathcal{L}_\omega = \det(D - \omega L)^{-1}\det\bigl((1 - \omega)D + \omega U\bigr) = (1 - \omega)^n .

行列式は重複を込めた固有値の積です（系 4.5［固有値と固有ベクトル］。複素数上で特性多項式は一次式の積に分解する）。nn 個の固有値の絶対値の積が ∣1−ω∣n|1 - \omega|^n なので、そのうち最大のものは ∣1−ω∣|1 - \omega| 以上です。後半は 定理 2.2 の 2 から、∣ω−1∣<1|\omega - 1| < 1 が必要だからです。

例 3.6（Jacobi 法は発散し Gauss–Seidel 法は収束する例）

A=(2−11222−1−12)A = \begin{pmatrix} 2 & -1 & 1 \\ 2 & 2 & 2 \\ -1 & -1 & 2 \end{pmatrix}

は対角優位でも正定値（対称）でもありません。Jacobi 法の反復行列は

GJ=(01/2−1/2−10−11/21/20),det⁡(tI−GJ)=t3+54tG_{\mathrm{J}} = \begin{pmatrix} 0 & 1/2 & -1/2 \\ -1 & 0 & -1 \\ 1/2 & 1/2 & 0 \end{pmatrix},\qquad \det(tI - G_{\mathrm{J}}) = t^3 + \tfrac{5}{4}t

で、固有値は 0,±i5/20, \pm i\sqrt{5}/2、ρ(GJ)=5/2≈1.118>1\rho(G_{\mathrm{J}}) = \sqrt{5}/2 \approx 1.118 > 1 なので発散します。一方 Gauss–Seidel 法では

GGS=(D−L)−1U=(01/2−1/20−1/2−1/200−1/2)G_{\mathrm{GS}} = (D - L)^{-1}U = \begin{pmatrix} 0 & 1/2 & -1/2 \\ 0 & -1/2 & -1/2 \\ 0 & 0 & -1/2 \end{pmatrix}

が上三角で、固有値は 0,−1/2,−1/20, -1/2, -1/2、ρ(GGS)=1/2\rho(G_{\mathrm{GS}}) = 1/2 なので収束します（GGSG_{\mathrm{GS}} は (D−L)GGS=U(D - L)G_{\mathrm{GS}} = U を第 1 行から順に解けば求まります）。どちらが速いかは行列しだいで、逆向きの例（Jacobi 法だけが収束する）もあります。

4. 矛盾なく順序づけられた行列と Young の定理

Section titled “4. 矛盾なく順序づけられた行列と Young の定理”

SOR 法の ω\omega をどう選ぶかに答えるには、Lω\mathcal{L}_\omega の固有値を Jacobi 法の固有値で表す必要があります。それが可能になる行列の構造が次の性質です。以下 A=D−L−UA = D - L - U、DD は正則とし、α≠0\alpha \ne 0 に対して

B(α)=αD−1L+α−1D−1UB(\alpha) = \alpha D^{-1}L + \alpha^{-1}D^{-1}U

とおきます。B(1)=GJB(1) = G_{\mathrm{J}} です。

SOR 法の反復行列には LL と UU が別々の重みで入っています。これを Jacobi 法の反復行列と結びつけるには、LL の部分を α\alpha 倍、UU の部分を α−1\alpha^{-1} 倍しても固有値が変わらない、という性質があると都合がよくなります。次の定義はその性質を取り出したもので、α\alpha はこのためだけに導入する補助のパラメータです。具体的には、三重対角行列や格子の行列がこの性質をもつことを、次の補題で確かめます。

定義 4.1（矛盾なく順序づけられた行列）↗

すべての α≠0\alpha \ne 0 について B(α)B(\alpha) の特性多項式 det⁡(tI−B(α))\det(tI - B(\alpha)) が α\alpha によらないとき、AA は矛盾なく順序づけられている（consistently ordered）という。

定義は抽象的ですが、次の補題で判定できます。

補題 4.2（番号づけによる判定）✓

整数 γ(1),…,γ(n)\gamma(1), \ldots, \gamma(n) があって、i≠ji \ne j、aij≠0a_{ij} \ne 0 のとき、j<ij < i なら γ(i)−γ(j)=1\gamma(i) - \gamma(j) = 1、j>ij > i なら γ(i)−γ(j)=−1\gamma(i) - \gamma(j) = -1 が成り立つとする。このとき AA は矛盾なく順序づけられている。とくに次の行列はそうである。

  1. 三重対角行列（γ(i)=i\gamma(i) = i）。
  2. m×mm\times m の格子の点 (p,q)(p, q)（1≤p,q≤m1 \le p, q \le m）を辞書式（自然な順序）に並べ、上下左右の隣接点の間にだけ非対角成分をもつ行列（γ(p,q)=p+q\gamma(p, q) = p + q）。
  3. 2 と同じ格子で、p+qp + q が偶数の点（赤）をすべて先に、奇数の点（黒）を後に並べた行列（赤黒順序。γ=0\gamma = 0（赤）、11（黒））。
証明（補題 4.2）

S=diag⁡(αγ(1),…,αγ(n))S = \operatorname{diag}(\alpha^{\gamma(1)}, \ldots, \alpha^{\gamma(n)}) とおくと、SB(1)S−1S B(1) S^{-1} の (i,j)(i, j) 成分は αγ(i)−γ(j)\alpha^{\gamma(i) - \gamma(j)} と B(1)B(1) の (i,j)(i, j) 成分の積です。B(1)B(1) の (i,j)(i, j) 成分は j<ij < i なら D−1LD^{-1}L の、j>ij > i なら D−1UD^{-1}U の成分で、aij≠0a_{ij} \ne 0 のときだけ零でありません。仮定からその位置で αγ(i)−γ(j)\alpha^{\gamma(i) - \gamma(j)} は j<ij < i なら α\alpha、j>ij > i なら α−1\alpha^{-1} なので、SB(1)S−1=B(α)SB(1)S^{-1} = B(\alpha) です。相似な行列の特性多項式は等しい（積の定理 定理 6.1［行列式とその性質］ から det⁡(tI−SBS−1)=det⁡Sdet⁡(tI−B)det⁡S−1\det(tI - SB S^{-1}) = \det S\det(tI - B)\det S^{-1}）ので、主張が従います。

1 では非零の非対角成分は j=i±1j = i \pm 1 だけです。2 では隣接点 (p′,q′)(p', q') は (p±1,q)(p \pm 1, q) か (p,q±1)(p, q \pm 1) で、辞書式の順序で前にあるのは (p−1,q)(p - 1, q) と (p,q−1)(p, q - 1)、すなわち γ\gamma が 11 小さい点です。3 では赤の点の隣接点はすべて黒、黒の隣接点はすべて赤（p+qp + q の偶奇が変わる）で、黒は赤より後に並ぶので、j<ij < i の非零成分は（黒の行，赤の列）、j>ij > i は（赤の行，黒の列）に限られます。

赤黒
格子点を市松模様に塗り分けた赤黒順序。赤の点（p + q が偶数）の上下左右の隣は必ず黒で、黒の隣は必ず赤である。

定理 4.3（Young の定理）✓

AA は矛盾なく順序づけられているとし、ω≠0\omega \ne 0 とする。

  1. μ\mu が GJG_{\mathrm{J}} の固有値なら −μ-\mu もそうである。
  2. μ\mu が GJG_{\mathrm{J}} の固有値で、λ\lambda が
(λ+ω−1)2=λ ω2μ2(\lambda + \omega - 1)^2 = \lambda\,\omega^2\mu^2

を満たすなら、λ\lambda は Lω\mathcal{L}_\omega の固有値である。逆に λ≠0\lambda \ne 0 が Lω\mathcal{L}_\omega の固有値なら、上の関係を満たす GJG_{\mathrm{J}} の固有値 μ\mu がある。 3. とくに ρ(GGS)=ρ(GJ)2\rho(G_{\mathrm{GS}}) = \rho(G_{\mathrm{J}})^2 である。

証明の方針は次のとおりです。det⁡(λI−Lω)\det(\lambda I - \mathcal{L}_\omega) を、B(α)B(\alpha) の特性多項式で書き直します。そのとき LL の係数が λω\lambda\omega、UU の係数が ω\omega と食い違うので、λ\lambda の平方根 ss を使って LL に ss、UU に s−1s^{-1} を掛けた形に直し、共通の因子 ωs\omega s をくくり出します。こうすると B(s)B(s) が現れ、矛盾なく順序づけられているという仮定で ss によらない特性多項式に置き換えられます。

証明（定理 4.3）

p(t)=det⁡(tI−GJ)p(t) = \det(tI - G_{\mathrm{J}}) とおくと、定義から det⁡(tI−B(α))=p(t)\det(tI - B(\alpha)) = p(t) がすべての α≠0\alpha \ne 0 で成り立ちます。

  1. B(−1)=−GJB(-1) = -G_{\mathrm{J}} なので p(t)=det⁡(tI+GJ)=(−1)np(−t)p(t) = \det(tI + G_{\mathrm{J}}) = (-1)^n p(-t) で、p(μ)=0p(\mu) = 0 なら p(−μ)=0p(-\mu) = 0 です。

  2. 積の定理（定理 6.1［行列式とその性質］）と det⁡(D−ωL)=det⁡D\det(D - \omega L) = \det D（三角行列、命題 3.3［行列式とその性質］ の 2）から

det⁡(λI−Lω)=det⁡(λ(D−ωL)−(1−ω)D−ωU)det⁡D=det⁡((λ+ω−1)I−λωD−1L−ωD−1U)\det(\lambda I - \mathcal{L}_\omega) = \frac{\det\bigl(\lambda(D - \omega L) - (1 - \omega)D - \omega U\bigr)}{\det D} = \det\bigl((\lambda + \omega - 1)I - \lambda\omega D^{-1}L - \omega D^{-1}U\bigr)

です。λ≠0\lambda \ne 0 のとき、λ\lambda の平方根の 1 つを ss とすると λωD−1L+ωD−1U=ωs B(s)\lambda\omega D^{-1}L + \omega D^{-1}U = \omega s\,B(s) なので

det⁡(λI−Lω)=det⁡((λ+ω−1)I−ωsB(s))=(ωs)n p(λ+ω−1ωs)\det(\lambda I - \mathcal{L}_\omega) = \det\bigl((\lambda + \omega - 1)I - \omega sB(s)\bigr) = (\omega s)^n\,p\Bigl(\frac{\lambda + \omega - 1}{\omega s}\Bigr)

です。よって λ≠0\lambda \ne 0 が Lω\mathcal{L}_\omega の固有値であることと、ν=(λ+ω−1)/(ωs)\nu = (\lambda + \omega - 1)/(\omega s) が GJG_{\mathrm{J}} の固有値であることは同値です。ν\nu が固有値なら μ=ν\mu = \nu が関係式を満たし（両辺を 2 乗）、逆に関係式を満たす固有値 μ\mu があれば ν=±μ\nu = \pm\mu で、1 から ν\nu も固有値です。これで λ≠0\lambda \ne 0 の場合が示せました。λ=0\lambda = 0 が関係式を満たすのは ω=1\omega = 1 のときで、このとき L1=(D−L)−1U\mathcal{L}_1 = (D - L)^{-1}U は det⁡U=0\det U = 0（UU は対角成分が 00 の三角行列）から特異なので、00 は固有値です。

  1. ω=1\omega = 1 では関係式は λ2=λμ2\lambda^2 = \lambda\mu^2 で、零でない固有値は λ=μ2\lambda = \mu^2（μ\mu は GJG_{\mathrm{J}} の固有値）と一致します。L1\mathcal{L}_1 の固有値は 00 と {μ2}\{\mu^2\} の部分集合で、2 の前半から各 μ2\mu^2 は実際に固有値なので、ρ(GGS)=max⁡∣μ∣2=ρ(GJ)2\rho(G_{\mathrm{GS}}) = \max|\mu|^2 = \rho(G_{\mathrm{J}})^2 です。

3 の意味は、Gauss–Seidel 法は Jacobi 法のちょうど 2 倍の漸近収束率をもつ、ということです（−ln⁡ρ2=−2ln⁡ρ-\ln\rho^2 = -2\ln\rho）。さらに ω\omega を最適に選ぶと、収束率は桁違いに改善します。

系 4.4（最適な緩和係数）✓

AA は矛盾なく順序づけられており、GJG_{\mathrm{J}} の固有値はすべて実数で ρJ=ρ(GJ)<1\rho_{\mathrm{J}} = \rho(G_{\mathrm{J}}) < 1 とする。0<ω<20 < \omega < 2 について

ρ(Lω)={14(ωρJ+ω2ρJ2−4(ω−1))2(0<ω≤ωopt)ω−1(ωopt≤ω<2),ωopt=21+1−ρJ2\rho(\mathcal{L}_\omega) = \begin{cases} \dfrac{1}{4}\Bigl(\omega\rho_{\mathrm{J}} + \sqrt{\omega^2\rho_{\mathrm{J}}^2 - 4(\omega - 1)}\Bigr)^2 & (0 < \omega \le \omega_{\mathrm{opt}}) \\[2mm] \omega - 1 & (\omega_{\mathrm{opt}} \le \omega < 2) \end{cases},\qquad \omega_{\mathrm{opt}} = \frac{2}{1 + \sqrt{1 - \rho_{\mathrm{J}}^2}}

であり、ρ(Lω)\rho(\mathcal{L}_\omega) は ω=ωopt\omega = \omega_{\mathrm{opt}} でだけ最小値 ωopt−1\omega_{\mathrm{opt}} - 1 をとる。

証明の方針は次のとおりです。Young の定理により、Lω\mathcal{L}_\omega の固有値は Jacobi 法の固有値 μ\mu ごとに、λ\lambda の 2 次方程式の根として決まります。まず 1 つの μ\mu について根の絶対値を調べ、それが ∣μ∣|\mu| とともに増えることから、最大の ρJ\rho_{\mathrm{J}} の場合だけを見ればよいと分かります。ω\omega が小さいうちは根が実数で、ω\omega を大きくすると下がります。ωopt\omega_{\mathrm{opt}} を超えると根が複素数になり、絶対値が ω−1\omega - 1 に固定されて今度は ω\omega とともに増えるので、境目の ωopt\omega_{\mathrm{opt}} が最小になります。

証明（系 4.4）

1 つの μ\mu に対する根. 実数 μ\mu を固定し、m=∣μ∣m = |\mu| とします。関係式 (λ+ω−1)2=λω2μ2(\lambda + \omega - 1)^2 = \lambda\omega^2\mu^2 は λ\lambda の 2 次方程式

λ2−(ω2m2−2(ω−1))λ+(ω−1)2=0\lambda^2 - \bigl(\omega^2m^2 - 2(\omega - 1)\bigr)\lambda + (\omega - 1)^2 = 0

で、判別式は ω2m2 d\omega^2m^2\,d、d=ω2m2−4(ω−1)d = \omega^2m^2 - 4(\omega - 1) です。

  • d<0d < 0 のとき（このとき ω>1\omega > 1）、2 根は共役な複素数で、絶対値の 2 乗は根の積 (ω−1)2(\omega - 1)^2 に等しいので、どちらも絶対値 ω−1\omega - 1 です。
  • d≥0d \ge 0 のとき、2 根は実数で、積 (ω−1)2≥0(\omega - 1)^2 \ge 0 と和 ω2m2−2(ω−1)≥0\omega^2m^2 - 2(\omega - 1) \ge 0（ω≥1\omega \ge 1 なら d≥0d \ge 0 から和 ≥2(ω−1)\ge 2(\omega - 1)、ω<1\omega < 1 なら各項が非負）から非負です。大きい方の根は、展開して確かめられるとおり ((ωm+d)/2)2\bigl((\omega m + \sqrt{d})/2\bigr)^2 です。

そこで g(ω,m)g(\omega, m) を、d≥0d \ge 0 なら ((ωm+d)/2)2\bigl((\omega m + \sqrt{d})/2\bigr)^2、d<0d < 0 なら ω−1\omega - 1 と定めると、gg は 2 根の絶対値の大きい方です。ω\omega を固定すると、d≥0d \ge 0 の範囲では ωm\omega m も dd も mm について増加し、境界 d=0d = 0 では g=ω2m2/4=ω−1g = \omega^2m^2/4 = \omega - 1 なので、g(ω,⋅)g(\omega, \cdot) は m≥0m \ge 0 について非減少です。

ρ(Lω)=g(ω,ρJ)\rho(\mathcal{L}_\omega) = g(\omega, \rho_{\mathrm{J}}). 定理 4.3 の 2 の後半から、Lω\mathcal{L}_\omega の零でない固有値は、ある固有値 μ\mu（実数で ∣μ∣≤ρJ|\mu| \le \rho_{\mathrm{J}}）についての根なので、絶対値は g(ω,∣μ∣)≤g(ω,ρJ)g(\omega, |\mu|) \le g(\omega, \rho_{\mathrm{J}}) 以下です。一方 ρJ\rho_{\mathrm{J}} または −ρJ-\rho_{\mathrm{J}} は固有値なので、2 の前半から絶対値 g(ω,ρJ)g(\omega, \rho_{\mathrm{J}}) の根は Lω\mathcal{L}_\omega の固有値です。

ω\omega についての形. ρ=ρJ\rho = \rho_{\mathrm{J}} とし、d(ω)=ρ2ω2−4ω+4d(\omega) = \rho^2\omega^2 - 4\omega + 4 とおきます。dd の零点は ω=2/(1±1−ρ2)\omega = 2/(1 \pm \sqrt{1 - \rho^2}) で、小さい方が ωopt∈[1,2)\omega_{\mathrm{opt}} \in [1, 2)、大きい方は 22 以上です（ρ=0\rho = 0 なら d=4−4ωd = 4 - 4\omega で零点は 1=ωopt1 = \omega_{\mathrm{opt}} だけ）。よって 0<ω≤ωopt0 < \omega \le \omega_{\mathrm{opt}} で d≥0d \ge 0、ωopt≤ω<2\omega_{\mathrm{opt}} \le \omega < 2 で d≤0d \le 0 となり、主張の式を得ます（ω=ωopt\omega = \omega_{\mathrm{opt}} では両方の式が ω−1\omega - 1 に一致します）。

単調性. ωopt≤ω<2\omega_{\mathrm{opt}} \le \omega < 2 では ω−1\omega - 1 は増加です。0<ω<ωopt0 < \omega < \omega_{\mathrm{opt}} では d>0d > 0 で、s(ω)=(ωρ+d)/2>0s(\omega) = (\omega\rho + \sqrt{d})/2 > 0 の導関数は

s′(ω)=12(ρ−2−ρ2ωd)s'(\omega) = \frac{1}{2}\Bigl(\rho - \frac{2 - \rho^2\omega}{\sqrt{d}}\Bigr)

です。(2−ρ2ω)2−ρ2d=4(1−ρ2)>0(2 - \rho^2\omega)^2 - \rho^2d = 4(1 - \rho^2) > 0 と 2−ρ2ω>02 - \rho^2\omega > 0 から 2−ρ2ω>ρd2 - \rho^2\omega > \rho\sqrt{d} なので s′<0s' < 0 で、s2s^2 は減少です。s2s^2 は ωopt\omega_{\mathrm{opt}} で連続に ωopt−1\omega_{\mathrm{opt}} - 1 につながるので、最小は ωopt\omega_{\mathrm{opt}} でだけとられます。

00.5100.511.52ωω_optGauss–Seidel（0.854）最小値 ω_opt − 1 ≈ 0.446
ρ_J = cos(π/8) ≈ 0.924（7×7 格子の離散 Laplace 行列）のときの ρ(𝓛_ω)。ω_opt ≈ 1.446 の左では急に下がり（平方根の特異性）、右では傾き 1 の直線 ω − 1 になる。ω を最適値より少し大きめにとる方が安全である。

5. モデル問題での収束率

Section titled “5. モデル問題での収束率”

正方形 (0,1)2(0,1)^2 上の −Δu=f-\Delta u = f、境界で u=0u = 0 を、格子幅 h=1/(m+1)h = 1/(m+1) の 5 点差分で離散化します。内部の格子点 (p,q)(p, q)（1≤p,q≤m1 \le p, q \le m）を自然な順序に並べると、係数行列は（h2h^2 を掛けて）

A=Im⊗T+T⊗Im,T=(2−1−12⋱⋱⋱−1−12)∈Rm×mA = I_m \otimes T + T \otimes I_m,\qquad T = \begin{pmatrix} 2 & -1 & & \\ -1 & 2 & \ddots & \\ & \ddots & \ddots & -1 \\ & & -1 & 2 \end{pmatrix} \in \mathbb{R}^{m\times m}

です（⊗\otimes は Kronecker 積）。1 次元のモデル問題は A=TA = T です。

命題 5.1（離散 Laplace 方程式の収束率）✓

1 次元（A=TA = T）と 2 次元（A=I⊗T+T⊗IA = I \otimes T + T \otimes I）のどちらでも、h=1/(m+1)h = 1/(m+1) として

ρJ=cos⁡πh,ρGS=cos⁡2πh,ωopt=21+sin⁡πh,ρ(Lωopt)=1−sin⁡πh1+sin⁡πh\rho_{\mathrm{J}} = \cos\pi h,\qquad \rho_{\mathrm{GS}} = \cos^2\pi h,\qquad \omega_{\mathrm{opt}} = \frac{2}{1 + \sin\pi h},\qquad \rho(\mathcal{L}_{\omega_{\mathrm{opt}}}) = \frac{1 - \sin\pi h}{1 + \sin\pi h}

が成り立つ。2 次元では自然な順序でも赤黒順序でも同じである。h→0h \to 0 で −ln⁡ρJ≈π2h2/2-\ln\rho_{\mathrm{J}} \approx \pi^2h^2/2、−ln⁡ρGS≈π2h2-\ln\rho_{\mathrm{GS}} \approx \pi^2h^2、−ln⁡ρ(Lωopt)≈2πh-\ln\rho(\mathcal{L}_{\omega_{\mathrm{opt}}}) \approx 2\pi h である。

証明（命題 5.1）

Jacobi 法の固有値. θk=kπh\theta_k = k\pi h（k=1,…,mk = 1, \ldots, m）とおきます。三重対角 Toeplitz 行列の固有値の補題（補題 8.3［行列のノルムと線形方程式の条件数］）を α=2\alpha = 2, β=−1\beta = -1 で使うと、TT の固有値は相異なる τk=2−2cos⁡θk\tau_k = 2 - 2\cos\theta_k で、固有ベクトルは vk=(sin⁡jθk)j\boldsymbol{v}_k = (\sin j\theta_k)_j です。

1 次元では D=2ID = 2I なので GJ=I−T/2G_{\mathrm{J}} = I - T/2 の固有値は 1−τk/2=cos⁡θk1 - \tau_k/2 = \cos\theta_k です。

2 次元では Kronecker 積の性質 (P⊗Q)(x⊗y)=Px⊗Qy(P \otimes Q)(\boldsymbol{x} \otimes \boldsymbol{y}) = P\boldsymbol{x} \otimes Q\boldsymbol{y}（定義から成分を比べればよい）により

A(vj⊗vk)=vj⊗Tvk+Tvj⊗vk=(τj+τk) vj⊗vkA(\boldsymbol{v}_j \otimes \boldsymbol{v}_k) = \boldsymbol{v}_j \otimes T\boldsymbol{v}_k + T\boldsymbol{v}_j \otimes \boldsymbol{v}_k = (\tau_j + \tau_k)\,\boldsymbol{v}_j \otimes \boldsymbol{v}_k

です。TT は対称なので相異なる固有値の固有ベクトルは直交し（τjvjTvk=(Tvj)Tvk=vjTTvk=τkvjTvk\tau_j\boldsymbol{v}_j^{\mathsf{T}}\boldsymbol{v}_k = (T\boldsymbol{v}_j)^{\mathsf{T}}\boldsymbol{v}_k = \boldsymbol{v}_j^{\mathsf{T}}T\boldsymbol{v}_k = \tau_k\boldsymbol{v}_j^{\mathsf{T}}\boldsymbol{v}_k）、(x⊗y)T(x′⊗y′)=(xTx′)(yTy′)(\boldsymbol{x} \otimes \boldsymbol{y})^{\mathsf{T}}(\boldsymbol{x}' \otimes \boldsymbol{y}') = (\boldsymbol{x}^{\mathsf{T}}\boldsymbol{x}')(\boldsymbol{y}^{\mathsf{T}}\boldsymbol{y}') なので、m2m^2 個の vj⊗vk\boldsymbol{v}_j \otimes \boldsymbol{v}_k は零でなく互いに直交します。直交する零でないベクトルは一次独立なので、これは Rm2\mathbb{R}^{m^2} の基底で、AA の固有値は τj+τk\tau_j + \tau_k で尽きます。D=4ID = 4I なので GJ=I−A/4G_{\mathrm{J}} = I - A/4 の固有値は (cos⁡θj+cos⁡θk)/2(\cos\theta_j + \cos\theta_k)/2 です。

どちらの場合も固有値は実数で、絶対値の最大は j=k=1j = k = 1（または j=k=mj = k = m）の cos⁡πh\cos\pi h です。2 次元の赤黒順序は自然な順序の行列を置換行列 PP で PTAPP^{\mathsf{T}}AP と並べ替えたもので、対角成分はすべて 44 のままなので GJG_{\mathrm{J}} も PTGJPP^{\mathsf{T}}G_{\mathrm{J}}P になり、固有値は変わりません。

SOR 法. 1 次元の TT は三重対角、2 次元は自然な順序でも赤黒順序でも 補題 4.2 の条件を満たすので、矛盾なく順序づけられています。定理 4.3 の 3 から ρGS=cos⁡2πh\rho_{\mathrm{GS}} = \cos^2\pi h、系 4.4 から ωopt=2/(1+1−cos⁡2πh)=2/(1+sin⁡πh)\omega_{\mathrm{opt}} = 2/(1 + \sqrt{1 - \cos^2\pi h}) = 2/(1 + \sin\pi h) と ρ(Lωopt)=ωopt−1=(1−sin⁡πh)/(1+sin⁡πh)\rho(\mathcal{L}_{\omega_{\mathrm{opt}}}) = \omega_{\mathrm{opt}} - 1 = (1 - \sin\pi h)/(1 + \sin\pi h) です。

漸近形. ln⁡cos⁡x=−x2/2+O(x4)\ln\cos x = -x^2/2 + O(x^4)、ln⁡((1−s)/(1+s))=−2s+O(s3)\ln\bigl((1 - s)/(1 + s)\bigr) = -2s + O(s^3) と sin⁡πh=πh+O(h3)\sin\pi h = \pi h + O(h^3) から従います。

誤差を 10−d10^{-d} 倍にする反復回数は、Jacobi 法で約 2dln⁡10/(π2h2)2d\ln 10/(\pi^2h^2)、Gauss–Seidel 法でその半分、最適 SOR 法で約 dln⁡10/(2πh)d\ln 10/(2\pi h) です。格子を 2 倍細かくすると、前の 2 つは 4 倍、最適 SOR 法は 2 倍に増えます。1 回の反復の手間は未知数の個数 m2≈h−2m^2 \approx h^{-2} に比例するので、全体の手間は h−4h^{-4} と h−3h^{-3} です。

次のコードで確かめます。Gauss–Seidel 法と SOR 法は、配列演算で書けるよう赤黒順序で実装しました（命題 5.1 により収束率は自然な順序と同じです）。

# Python 3.12.10 / NumPy 2.5.2（x86-64）
import numpy as np


def nbr(u):                 # 上下左右の和（境界は 0）
    return u[:-2, 1:-1] + u[2:, 1:-1] + u[1:-1, :-2] + u[1:-1, 2:]


def count(m, method, w=1.0, tol=1e-8):
    b = np.full((m, m), 1.0 / (m + 1) ** 2)     # f = 1 の右辺
    u = np.zeros((m + 2, m + 2))
    red = np.indices((m, m)).sum(axis=0) % 2 == 0
    for k in range(1, 10**6):
        if method == "J":
            u[1:-1, 1:-1] = (b + nbr(u)) / 4
        else:                # 赤黒順序の SOR（w = 1 で Gauss-Seidel）
            for c in (red, ~red):
                v = u[1:-1, 1:-1]
                v[c] = (1 - w) * v[c] + w * ((b + nbr(u)) / 4)[c]
        r = b - 4 * u[1:-1, 1:-1] + nbr(u)
        if np.linalg.norm(r) <= tol * np.linalg.norm(b):
            return k


print(" 1/h    J(予測)       GS(予測)     SOR(予測)")
for m in (15, 31, 63):
    h = 1 / (m + 1)
    rj = np.cos(np.pi * h)
    w = 2 / (1 + np.sin(np.pi * h))
    pr = [np.log(1e-8) / np.log(x) for x in (rj, rj**2, w - 1)]
    k = [count(m, "J"), count(m, "S"), count(m, "S", w)]
    print(f"{m+1:4d}" + "".join(
        f"  {a:5d}({p:5.0f})" for a, p in zip(k, pr)))
 1/h    J(予測)       GS(予測)     SOR(予測)
  16    942(  949)    480(  475)     61(   47)
  32   3779( 3816)   1926( 1908)    125(   94)
  64  15122(15283)   7705( 7642)    257(  188)

検算. (a) 予測は閉じた形の ρ\rho による ln⁡10−8/ln⁡ρ\ln 10^{-8}/\ln\rho です。(b) 反復行列のスペクトル半径そのものは次のコードで固有値を直接計算して照合します。(c) 格子幅を半分にするごとの回数の比は、Jacobi 法で 4.01,4.004.01, 4.00、Gauss–Seidel 法で 4.01,4.004.01, 4.00（h−2h^{-2}）、最適 SOR 法で 2.05,2.062.05, 2.06（h−1h^{-1}）で、命題 5.1 の漸近形と合います。Jacobi 法・Gauss–Seidel 法の回数は予測と 1〜2% で一致します。最適 SOR 法の回数は予測の 1.3〜1.4 倍ですが、これは ωopt\omega_{\mathrm{opt}} では Lω\mathcal{L}_\omega の固有値が重根になり（系 4.4 の証明の d=0d = 0）、誤差が ρk\rho^k ではなく kρkk\rho^k のように減る成分をもつためです（注意 5.2）。回数が h−1h^{-1} に比例することは変わりません。

# Python 3.12.10 / NumPy 2.5.2（x86-64, OpenBLAS）
import numpy as np


m = 7                                  # 内部の格子点 7×7、h = 1/8
h = 1 / (m + 1)
T = 2 * np.eye(m) - np.eye(m, k=1) - np.eye(m, k=-1)
A = np.kron(np.eye(m), T) + np.kron(T, np.eye(m))   # 自然な順序
D = np.diag(np.diag(A))
L = -np.tril(A, -1)
U = -np.triu(A, 1)


def rho(G):
    return max(abs(np.linalg.eigvals(G)))


def sor(w):
    return np.linalg.solve(D - w * L, (1 - w) * D + w * U)


rj = np.cos(np.pi * h)
wopt = 2 / (1 + np.sin(np.pi * h))
rows = [("Jacobi", np.linalg.solve(D, L + U), rj),
        ("GS", sor(1.0), rj**2),
        ("SOR w=1.5", sor(1.5), 0.5),
        ("SOR w_opt", sor(wopt), wopt - 1)]
for name, G, th in rows:
    t = f"{th:.6f}"
    print(f"{name:10s} {rho(G):.6f}  {t}")
Jacobi     0.923880  0.923880
GS         0.853553  0.853553
SOR w=1.5  0.500000  0.500000
SOR w_opt  0.446463  0.446463

検算. (a) 右の列は 命題 5.1 の閉じた形と、ω=1.5>ωopt≈1.446\omega = 1.5 > \omega_{\mathrm{opt}} \approx 1.446 での 系 4.4 の値 ω−1\omega - 1 です。(b) 左の列は自然な順序の反復行列を密行列として作り、numpy.linalg.eigvals で求めた値で、6 桁一致します。(c) 4 行は 系 4.4 の曲線（上の図）の 3 点 ω=1,ωopt,1.5\omega = 1, \omega_{\mathrm{opt}}, 1.5 にあたります。ωopt\omega_{\mathrm{opt}} の行列は重根をもつため固有値の計算は丸めに敏感ですが、表示した 6 桁は既定のスレッド数と 1 スレッドで一致しました。

注意 5.2（最適 SOR 法の過渡的な振る舞い）

ω=ωopt\omega = \omega_{\mathrm{opt}} では、Jacobi 法の固有値 ±ρJ\pm\rho_{\mathrm{J}} に対応する Lω\mathcal{L}_\omega の固有値が重根になり、Lωopt\mathcal{L}_{\omega_{\mathrm{opt}}} は一般に対角化できません（2 次の Jordan 細胞をもつ）。このため ∥Lωk∥\|\mathcal{L}_\omega^k\| は ρk\rho^k ではなく kρk−1k\rho^{k-1} の程度で減り、実際の反復回数は ln⁡ε/ln⁡ρ\ln\varepsilon/\ln\rho の予測より多くなります。実務では ω\omega を ωopt\omega_{\mathrm{opt}} よりわずかに大きくとることが多く、ρJ\rho_{\mathrm{J}} が正確に分からないときも、系 4.4 の曲線の形（左側は急、右側は傾き 1）から大きめの側に外す方が損が小さくなります。Varga『Matrix Iterative Analysis』第 4 章、Young『Iterative Solution of Large Linear Systems』を参照してください。

注意 5.3（Chebyshev 加速）

対称な反復行列 GG の固有値が [−ρ,ρ][-\rho, \rho] にあると分かっているとき、x0,…,xk\boldsymbol{x}_0, \ldots, \boldsymbol{x}_k の 1 次結合で誤差を pk(G)e0p_k(G)\boldsymbol{e}_0（pk(1)=1p_k(1) = 1 の kk 次多項式）の形にできます。pkp_k を区間 [−ρ,ρ][-\rho, \rho] 上で最小最大にする Chebyshev 多項式に選ぶと、漸近収束率は最適 SOR 法と同じ h−1h^{-1} の程度まで改善します。AA が正定値 Hermite 行列なら、固有値の区間を知らずに同じことを自動で行うのが共役勾配法で、Krylov 部分空間法 I の章 の主題です。Golub–Van Loan『Matrix Computations』§11.2、Saad『Iterative Methods for Sparse Linear Systems』第 12 章を参照してください。

注意 5.4（多重格子法での平滑化）

命題 5.1 の証明の固有ベクトル vj⊗vk\boldsymbol{v}_j \otimes \boldsymbol{v}_k を見ると、Jacobi 法で遅く減るのは θj\theta_j, θk\theta_k が小さい（格子に対してなめらかな）成分だけです。振動の激しい成分は、減衰つき Jacobi 法（ω=4/5\omega = 4/5 など）や Gauss–Seidel 法の数回の反復で格子幅によらない割合で減ります。なめらかになった誤差を粗い格子で解き直すのが多重格子法で、反復回数が hh によらない解法になります。Briggs–Henson–McCormick『A Multigrid Tutorial』2nd ed., SIAM, 2000 を参照してください。

6. 演習

Section titled “6. 演習”

演習 6.1標準

1/2<a<11/2 < a < 1 とし、A=(1aaa1aaa1)A = \begin{pmatrix} 1 & a & a \\ a & 1 & a \\ a & a & 1 \end{pmatrix} とする。AA は正定値であるが Jacobi 法は発散し、Gauss–Seidel 法は収束することを示せ。補題 3.2 の条件 2D−A2D - A の正定値性はどうなっているか。

解答

A=(1−a)I+a11TA = (1 - a)I + a\boldsymbol{1}\boldsymbol{1}^{\mathsf{T}} で、1\boldsymbol{1} の方向の固有値は 1+2a1 + 2a、それに直交する 2 次元の空間では 1−a1 - a です。どちらも正なので AA は正定値で、系 3.4 から Gauss–Seidel 法は収束します。D=ID = I なので GJ=I−AG_{\mathrm{J}} = I - A の固有値は −2a-2a と aa（2 重）で、ρ(GJ)=2a>1\rho(G_{\mathrm{J}}) = 2a > 1 から Jacobi 法は発散します（定理 2.2）。2D−A=2I−A2D - A = 2I - A の固有値は 1−2a<01 - 2a < 0 と 1+a1 + a なので正定値でなく、補題の仮定が崩れています。

演習 6.2易

2 次元のモデル問題で h=1/100h = 1/100 とする。誤差を 10−610^{-6} 倍にするのに必要な反復回数を、Jacobi 法・Gauss–Seidel 法・最適 SOR 法について 命題 5.1 の ρ\rho から見積もれ。

解答

πh≈0.0314159\pi h \approx 0.0314159 です。ρJ=cos⁡πh≈0.9995066\rho_{\mathrm{J}} = \cos\pi h \approx 0.9995066 で −ln⁡ρJ≈4.936×10−4-\ln\rho_{\mathrm{J}} \approx 4.936\times10^{-4}、回数は 6ln⁡10/(4.936×10−4)≈13.82/(4.936×10−4)≈2.8×1046\ln 10/(4.936\times10^{-4}) \approx 13.82/(4.936\times10^{-4}) \approx 2.8\times10^4 です。Gauss–Seidel 法はその半分の約 1.4×1041.4\times10^4 回です。最適 SOR 法は sin⁡πh≈0.0314108\sin\pi h \approx 0.0314108、ρ=(1−0.0314108)/(1+0.0314108)≈0.939092\rho = (1 - 0.0314108)/(1 + 0.0314108) \approx 0.939092、−ln⁡ρ≈0.06284-\ln\rho \approx 0.06284 なので約 220220 回です。

演習 6.3標準

A=(2−1−12)A = \begin{pmatrix} 2 & -1 \\ -1 & 2 \end{pmatrix} について、ρJ\rho_{\mathrm{J}}、ρGS\rho_{\mathrm{GS}}、ωopt\omega_{\mathrm{opt}}、ρ(Lωopt)\rho(\mathcal{L}_{\omega_{\mathrm{opt}}}) を求めよ。さらに Lω\mathcal{L}_\omega の特性多項式を直接計算し、ω=ωopt\omega = \omega_{\mathrm{opt}} で重根になることを確かめよ。

解答

GJ=(01/21/20)G_{\mathrm{J}} = \begin{pmatrix} 0 & 1/2 \\ 1/2 & 0 \end{pmatrix} の固有値は ±1/2\pm 1/2 なので ρJ=1/2\rho_{\mathrm{J}} = 1/2 です。AA は三重対角なので矛盾なく順序づけられており（補題 4.2）、定理 4.3 から ρGS=1/4\rho_{\mathrm{GS}} = 1/4、系 4.4 から ωopt=2/(1+3/2)=4(2−3)≈1.0718\omega_{\mathrm{opt}} = 2/(1 + \sqrt{3}/2) = 4(2 - \sqrt{3}) \approx 1.0718、ρ(Lωopt)=7−43≈0.0718\rho(\mathcal{L}_{\omega_{\mathrm{opt}}}) = 7 - 4\sqrt{3} \approx 0.0718 です。直接計算すると D−ωL=(20−ω2)D - \omega L = \begin{pmatrix} 2 & 0 \\ -\omega & 2 \end{pmatrix}、(1−ω)D+ωU=(2(1−ω)ω02(1−ω))(1 - \omega)D + \omega U = \begin{pmatrix} 2(1 - \omega) & \omega \\ 0 & 2(1 - \omega) \end{pmatrix} で、

det⁡(λ(D−ωL)−(1−ω)D−ωU)=4(λ+ω−1)2−λω2\det\bigl(\lambda(D - \omega L) - (1 - \omega)D - \omega U\bigr) = 4(\lambda + \omega - 1)^2 - \lambda\omega^2

です。det⁡(D−ωL)=4\det(D - \omega L) = 4 で割ると特性多項式は (λ+ω−1)2−λω2/4(\lambda + \omega - 1)^2 - \lambda\omega^2/4 で、これは Young の関係式で μ2=1/4\mu^2 = 1/4 としたものです。λ\lambda の 2 次式として判別式は (ω2/4)(ω2/4−4(ω−1))(\omega^2/4)\bigl(\omega^2/4 - 4(\omega - 1)\bigr) で、ω2−16ω+16=0\omega^2 - 16\omega + 16 = 0 の小さい方の根 ω=8−43=4(2−3)\omega = 8 - 4\sqrt{3} = 4(2 - \sqrt{3}) で 00 になり、重根 λ=ω−1=7−43\lambda = \omega - 1 = 7 - 4\sqrt{3} をもちます。

参考文献

Section titled “参考文献”
  • R. S. Varga『Matrix Iterative Analysis』2nd ed., Springer, 2000 — 第 3 章（定常反復の基本）、第 4 章（SOR 法と矛盾なく順序づけられた行列）。
  • D. M. Young『Iterative Solution of Large Linear Systems』Academic Press, 1971 — SOR 法の理論の原典的な教科書。
  • Y. Saad『Iterative Methods for Sparse Linear Systems』2nd ed., SIAM, 2003 — 第 4 章（基本的な反復法）、第 12 章（多項式前処理。Chebyshev 加速を含む）。
  • G. H. Golub, C. F. Van Loan『Matrix Computations』4th ed., Johns Hopkins University Press, 2013 — §11.2（古典的な反復法）。
  • W. L. Briggs, V. E. Henson, S. F. McCormick『A Multigrid Tutorial』2nd ed., SIAM, 2000 — 第 2 章（平滑化と誤差の周波数）。

このテーマの 100 本ノック（100 問）

関連記事

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

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

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