連続最適化および関連分野に関する夏季学校
連続最適化への応用に向けた常微分方程式の数値解析入門
演習問題 解答例
佐藤 峻

2023 年 8 月 9 日
問題 1. Runge–Kutta 法の安定性関数が，R(z) =det(I − zA + z1bT )det(I − zA)を満たすことを示せ．
解答 1. まず，安定性関数 R を素直に計算する．Dahlquist のテスト方程式に Runge–Kutta 法を適用すると，

Xi = 1 + h
x1 = 1 + hs
X
j=1s
X

aij λXj
bi λXii=1となり，x1 を z = hλ の関数で表したものが安定性関数 R である．内部段 Xi を束ねたベクトルを X (つまhり，X = X1···iTXs ) と書くことにすると，第一式は，

X = 1 + λhAXと書ける．すなわち，X = (I − zA)−1 1 である．これを第二式に代入すると，

R(z) = 1 + zbT (I − zA)−1 1(1)を得る．
ブロック行列 M をM=1−bTz1 I − zA

と定め，M の行列式を 2 通りの方法で評価すると，det M = det(1) × det(I − zA + z1bT ),−1det M = det(I − zA) × det 1 + zbT (I − zA) 1を得るため，これを整理すると，題意を得る．
問題 2. 任意の 1 次以上の精度をもつ陽的 Runge–Kutta 法について，その安定領域は有界であることを示せ．
解答 2. 1 次以上の精度をもつ陽的 Runge–Kutta 法の安定性関数 R は，R(z) = 1 + z + α2 z 2 + · · · + αs z sとかける (|αs | ̸= 0 を満たすように s を選ぶ)．極座標表示 z = r exp(iθ) を用いると，|R(z)| ≥ |αs |rs − 1 + r exp(iθ) + · · · + αs−1 rs−1 exp((s − 1)iθ)≥ |αs |rs − 1 + r + |α2 |r2 + · · · + |αs−1 |rs−11が成立する．最右辺は，実係数の多項式であり，最高次の係数が正であるため，ある R が存在して，任意のr > R について，

|αs |rs − 1 + r + |α2 |r2 + · · · + |αs−1 |rs−1 > 1が成立する．これより，安定領域 S は原点を中心とし，半径 R の円に含まれるため，有界である (したがって，A 安定ではない)．
問題 3. µ-強凸 L-平滑関数に対する Nesterov の加速勾配法に相当する線形 2 段階法
x(k+2) − (1 + γ)x(k+1) + γx(k) = h (1 − γ 2 )g x(k+1) − γ(1 − γ)g x(k)の絶対安定領域 S ⊆ C について，S ∩ R を L, µ を用いて表せ (L ≥ µ > 0 とする)．ここで，γ :=pµ1−pLµ ∈ [0, 1)1+Lである．
解答 3. 安定領域 S の定義は，

S = { z ∈ C | π(ζ; z) が根条件を満たす．}であったため，

S ∩ R = { z ∈ R | π(ζ; z) が根条件を満たす．}
を考えればよい．与えられた線形 2 段階法について，

π(ζ; z) = ζ 2 − (1 + γ)(1 + z(1 − γ))ζ + γ(1 + z(1 − γ))である．この (ζ に関する) 2 次関数の根を調べればよい．記述を簡単にするために，ẑ = 1 + z(1 − γ) を導入し，π̂(ζ; ẑ) = ζ 2 − (1 + γ)ẑζ + γ ẑを考えることにする．
まずは，根条件の 2 つ目に当たる，重解の場合を処理しておく．重解になるのは，(1 + γ)2 ẑ 2 = 4γ ẑ の場合4γであり，これは，(A) ẑ = 0, (B) ẑ = (1+γ)
2 の 2 通りのみである．(A) の場合については，根も 0 であるため，2γ絶対値は 1 未満であり，考えなくてよい．(B) の場合については，根は 1+γであり，任意の γ ∈ [0, 1) に対して 1 未満であるため，考えなくてよい．
続いて，重解でない場合に，根の絶対値が 1 以下になる条件を考える．2 次関数 π̂ の根 ζ± は，

ζ± =(1 + γ)ẑ ±4γと書ける (複号同順)．ここで， (i) ẑ ∈ 0, (1+γ)
24γ(iii) ẑ > (1+γ)
2p(1 + γ)2 ẑ 2 − 4γ ẑ2のとき根 ζ± は (共役) 複素数であり，(ii) ẑ < 0 あるいはのときは実根である．それぞれの場合について検討する．

4γ(i) ẑ ∈ 0, (1+γ)
2の場合．

根の絶対値の 2 乗を考えると，

2|ζ± | =(1 + γ)ẑ2

2+4γ ẑ − (1 + γ)2 ẑ 2= γ ẑ44γ4γ1であるため，|ζ± | ≤ 1 は，ẑ ≤ 1/γ と同値である． (1+γ)
0, (1+γ)
2 < 1 < γ であるため，ẑ ∈2

す全ての ẑ において，根の絶対値は 1 以下である．

2を満た

(ii) ẑ < 0 の場合．
符号に注意すると，−(1 + γ)ẑ +max{|ζ+ |, |ζ− |} = |ζ− | =p(1 + γ)2 ẑ 2 − 4γ ẑ2h11である．これが 1 以下になる必要十分条件は，整理すると ẑ ≥ − 1+2γである．つまり，ẑ ∈ − 1+2γ,0を満たす全ての ẑ において，根の絶対値は 1 以下である．
4γ(iii) ẑ > (1+γ)
2 の場合．

符号に注意すると，

max{|ζ+ |, |ζ− |} = |ζ+ | =(1 + γ)ẑ +p(1 + γ)2 ẑ 2 − 4γ ẑ2である．これが 1 以下になる必要十分条件は，整理すると ẑ ≤ 1 である．つまり，ẑ ∈i4γ(1+γ)2 , 1を満たす全ての ẑ において，根の絶対値は 1 以下である．

hi

1以上の考察を総合すると，ẑ ∈ − 1+2γ
, 1 を満たす全ての ẑ において，根条件は満たされる．よって，ẑ の定義より，


1S∩R= z ∈R −≤ 1 + z(1 − γ) ≤ 1
1 + 2γ2(1 + γ)= z∈R −≤z≤0
(1 + 2γ)(1 − γ)s()pµ
2 1+ LLpµ= z∈R −≤z≤0 .
µ3− L問題 4. 写像 g : Rd → Rd は L-Lipschitz 連続とする．任意の x0 , y0 ∈ Rd に対して，x, y : [0, T ) → Rd が，ẋ(t) = g(x(t)) (t ∈ (0, T )),x(0) = x0 ,ẏ(t) = g(y(t)) (t ∈ (0, T )),y(0) = y0を満たすとき，∥x(t) − y(t)∥ ≤ exp(Lt)∥x0 − y0 ∥(t ∈ (0, T ))が成立することを示せ．
2解答 4. 関数 e(t) := ∥x(t) − y(t)∥ について，de(t) = 2⟨x(t) − y(t), ẋ(t) − ẏ(t)⟩dt= 2⟨x(t) − y(t), g(x(t)) − g(y(t))⟩
2≤ 2L∥x(t) − y(t)∥= 2Le(t)が成立する．よって，e(t) ≤ exp(2Lt)e(0) が成立する．関数 e の定義より，これは題意の成立を意味する．
問題 5. 2 段 2 次の陽的 Runge–Kutta 法を全て求めよ．また，具体的な計算により，どの手法についても，安2定性関数は R(z) = 1 + z + z2 であることを確認せよ．
解答 5. 2 段陽的 Runge–Kutta 法の非零パラメータは，a21 , b1 , b2 の 3 つである．このときの，2 次精度の条件は，b1 + b2 = 1,b2 a21 =12の 2 つである．よって，α をパラメータとして，

a21 = α,b1 = 1 −31,2αb2 =12αのようにパラメータを選べば，この Runge–Kutta 法は 2 次精度である．
安定性関数 R について，式 (1) を用いると，R(z) = 1 + z=1+z

11 − 2α12α11 − 2α12α


1−zα1zα−1  0111 0 11 1


11 + zα
1 + zα1=1+z 1−+2α2α1= 1 + z + z2.
2=1+z1
1 − 2α12α問題 6. 陽的な線形 3 段階法で，5 次精度を達成するものを求めよ．また，求めた線形 3 段階法が零安定でないことを示せ．(この問題については，次数条件を満たすパラメータを求める際に数値実験を用いてよい．)
解答 6. 陽的 3 段法は，α3 = 1 とおくと (一般性を失わない)，x(k+3) + α2 x(k+2) + α1 x(k+1) + α0 x(k) = β2 g x(k+2) + β1 g x(k+1) + β0 g x(k)とかける．5 次精度の次数条件を書き下すと，

1
0

0

0

00111111124816
32010000  
α0−10
  1 α1   −3 α2   −9 4  =   12  β0   −27 32  β1   −81 −24380
β2012345となる．これの解は，

α0 = −10,

α1 = −9,α2 = 18,β0 = −3,β1 = −18,β2 = −9√となる*1 ．このとき，第一特性多項式は ρ(z) = z 3 + 18z 2 − 9z − 10 である．この ρ の根は 1 と −19± 219 −402であるため，この線形 3 段法は零安定ではない．
問題 7. 凸関数に対する加速勾配法に対応する ODE

3ẍ + ẋ + ∇f (x) = 0,tx(0) = x0 , ẋ(0) = 0を数値的に解いて挙動を観察せよ．ただし，目的関数 f や，数値解法は自由に定めてよい．(汎用の ODE ソルバーを用いても構わない．ただし，その場合は，実装に用いられている手法を確認すること)．
解答 7. まずは，1 階の ODE に書き換える．つまり，v = ẋ を導入して，ẋ = v,3v̇ = − v − ∇f (x)t2を初期値 (x(0), v(0)) = (x0 , 0) の下で考える．目的関数は，f (x1 , x2 ) = 21 x21 + 10
2 x2 とする．このとき，具体的な ODE は， 0x1
d x 2 =  0dt  v1  −1v2001000− 3t−10 0 0x1
x2 1  0   v1 v2− 3t*1 手で解くのは大変だが，適当な線形方程式のソルバーに解いてもらって，その解を有理数 (今回は整数) に丸めて確認すればよい．

41101RK4 (h = 1/2)

RK4 (h = 1)RK4 (h = 1/5)
RK4(h = 1/10)

RK4(h = 1/10)

10−2
f (x(t))x20.5

RK4 (h = 1/5)

0−0.510−510−8−1−0.200.2

図10.4
x10.60.810−1110246RK4 の軌道

図2RK4 の軌道上の目的関数値10

陰的 Euler (h = 1)

陰的 Euler (h = 1)

陰的 Euler (h = 1/5)

陰的 Euler (h = 1/5)

陰的 Euler(h = 1/10)

陰的 Euler(h = 1/10)

10−2f (x(t))

x210110.58t0−0.5

10−510−8−1−0.20

図30.20.4
x10.60.810−1110246810t

陰的 Euler の軌道

図4

陰的 Euler の軌道上の目的関数値となる．ここに現れる行列の固有値は，−3 ±√9 − 4t2,2t−3 ±√9 − 40t22tの 4 つであり，t → ∞ でこれらの固有値は虚軸に漸近する*2 ．このため，陽的 Euler 法との相性が致命的に悪く，陽的 Euler 法を使うとステップ幅をかなり小さくしてもまともに動かない．
そこで，虚軸上を少し含む陽的解法である，4 段 4 次 Runge–Kutta 法 (RK4) を適用してみると，h = 1 では発散してしまうが，h をもう少し小さくすると，発散はせず，有限の解を出力する (図 1)*3 ．目的関数値も，単調減少ではないが，加速勾配法でよく見られるような，振動しながら減少していく様が観察される (図 2)．
続いて，陰的 Euler 法の数値実験も行う．陰的 Euler 法は，非常に安定な数値解法であるため，h = 1 でも計算は破綻しない．軌道 (図 3) を見ると，あまり振動せずに最小解に向かっているように見えるが，目的関数値(図 4) を見ると，RK4 と同様に振動しつつ減少していく様が観察された．

*2 t → 0 で負の実軸方向に発散するが，これは有限の大きさのステップ幅をとっている限り，あまり大きな問題にはならない (計算は破綻しない)．ただし，h → 0 としたときに厳密解に収束するかは怪しい．
*3 t = 0 で特異性があるが，これは，t = 0 で v = 0 と定めることで処理した．
t5