§E20.29線形多段法

最終更新

一段法は、次の近似値yn+1y_{n+1}を直前の近似値yny_nだけから定める。これに対して、連続するk+1k+1個の時刻における近似値と右辺の値の間に一つの一次関係を課し、それによって新しい近似値を定める方法を線形多段法という。θ 法はk=1k=1の場合としてこれに含まれる。たとえば 2 段 Adams–Bashforth 法

yn+2=yn+1+h2(3f(tn+1,yn+1)−f(tn,yn))y_{n+2}=y_{n+1}+\frac h2\bigl(3f(t_{n+1},y_{n+1})-f(t_n,y_n)\bigr)

では、n≥1n\ge1のときf(tn,yn)f(t_n,y_n)は前の歩で評価済みであり、一歩ごとに新たに評価する右辺の値はf(tn+1,yn+1)f(t_{n+1},y_{n+1})だけである。この点は前進 Euler 法と同じであるが、この方法の次数は22である。

一方で、線形多段法の近似値の誤差はkk項にわたる漸化式に従って伝わる。このため、厳密解を代入したときの一歩の食い違いが刻みの高い冪で小さくなることだけからは、刻みを細かくしたときに近似値が解へ近づくとは限らない。線形多段法を扱うには、一歩の精度と誤差の伝わり方を分けて調べる必要がある。

本記事では、線形多段法の次数と収束を係数から調べ、代表的な方法の安定性を一段法と比べる。

1 線形多段法と次数

定義 1.1.k∈N≥1k\in\NNとし、実数α0,…,αk,β0,…,βk\alpha_0,\dots,\alpha_k,\beta_0,\dots,\beta_kがαk=1\alpha_k=1を満たすとする。

  1. 係数の組(α0,…,αk;β0,…,βk)(\alpha_0,\dots,\alpha_k;\beta_0,\dots,\beta_k)をkk段の 線形多段法 (linear multistep method) といい、 ρ(ζ):=∑j=0kαjζj,σ(ζ):=∑j=0kβjζj\rho(\zeta):=\sum_{j=0}^k\alpha_j\zeta^j,\qquad\sigma(\zeta):=\sum_{j=0}^k\beta_j\zeta^j をそれぞれ線形多段法の 第一特性多項式 (first characteristic polynomial)、第二特性多項式 (second characteristic polynomial) という。βk=0\beta_k=0のとき線形多段法は 陽的 (explicit) であるといい、βk≠0\beta_k\ne0のとき 陰的 (implicit) であるという。
  2. d∈N≥1d\in\NNとし、I⊆RI\subseteq\Rを開区間、f ⁣:I×Rd→Rdf\colon I\times\R^d\to\R^dを連続写像とする。t0∈It_0\in Iとh>0h>0を取り、tn:=t0+nht_n:=t_0+nhと置き、tk−1∈It_{k-1}\in Iとする。y0,…,yk−1∈Rdy_0,\dots,y_{k-1}\in\R^dを取る。n∈N≥0n\in\Nについて、yn,…,yn+k−1y_n,\dots,y_{n+k-1}が定まり、tn+k∈It_{n+k}\in Iであり、v∈Rdv\in\R^dについての方程式 v+∑j=0k−1αjyn+j=hβkf(tn+k,v)+h∑j=0k−1βjf(tn+j,yn+j)v+\sum_{j=0}^{k-1}\alpha_jy_{n+j}=h\beta_kf(t_{n+k},v)+h\sum_{j=0}^{k-1}\beta_jf(t_{n+j},y_{n+j}) がただ一つの解をもつとき、その解をyn+ky_{n+k}とする。y0,…,yk−1y_0,\dots,y_{k-1}を 開始値 (starting values) といい、こうして定まるyny_nを線形多段法の近似値という。yn+ky_{n+k}が定まるとき、fm:=f(tm,ym)f_m:=f(t_m,y_m)と書くと∑j=0kαjyn+j=h∑j=0kβjfn+j\sum_{j=0}^k\alpha_jy_{n+j}=h\sum_{j=0}^k\beta_jf_{n+j}である。

定義 1.2.kk段の線形多段法(α0,…,αk;β0,…,βk)(\alpha_0,\dots,\alpha_k;\beta_0,\dots,\beta_k)を取る。

  1. d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。t0<Tt_0<Tとし、J⊇[t0,T]J\supseteq[t_0,T]を開区間、y ⁣:J→Rdy\colon J\to\R^dをC1C^1級写像とする。0<h≤(T−t0)/k0<h\le(T-t_0)/kとt∈[t0,T−kh]t\in[t_0,T-kh]に対して dh(t):=∑j=0kαjy(t+jh)−h∑j=0kβjy′(t+jh)d_h(t):=\sum_{j=0}^k\alpha_jy(t+jh)-h\sum_{j=0}^k\beta_jy'(t+jh) を、yyの点tt、刻みhhにおける 局所欠陥 (local defect) という。
  2. 00:=10^0:=1として C0:=∑j=0kαj,Cq:=∑j=0kαjjqq!−∑j=0kβjjq−1(q−1)!(q∈N≥1)C_0:=\sum_{j=0}^k\alpha_j,\qquad C_q:=\sum_{j=0}^k\frac{\alpha_jj^q}{q!}-\sum_{j=0}^k\frac{\beta_jj^{q-1}}{(q-1)!}\quad(q\in\NN) と置く。
  3. 任意のd∈N≥1d\in\NN、Rd\R^dの任意のノルム、任意のt0<Tt_0<T、[t0,T][t_0,T]を含む任意の開区間JJと任意のy∈C1(J;Rd)y\in C^1(J;\R^d)について lim⁡h→+0sup⁡{∥dh(t)∥h ∣ t∈[t0,T−kh]}=0\lim_{h\to+0}\sup\Bigl\{\frac{\|d_h(t)\|}{h}\ \Big|\ t\in[t_0,T-kh]\Bigr\}=0 が成り立つとき、線形多段法は 整合的 (consistent) であるという。
  4. p∈N≥1p\in\NNとする。任意のd∈N≥1d\in\NN、Rd\R^dの任意のノルム、任意のt0<Tt_0<T、[t0,T][t_0,T]を含む任意の開区間JJと任意のy∈Cp+1(J;Rd)y\in C^{p+1}(J;\R^d)に対して、C≥0C\ge0が存在し、0<h≤(T−t0)/k0<h\le(T-t_0)/kとt∈[t0,T−kh]t\in[t_0,T-kh]を満たす任意のh,th,tについて∥dh(t)∥≤Chp+1\|d_h(t)\|\le Ch^{p+1}が成り立つとき、線形多段法の 次数 (order) はpp以上であるという。次数がpp以上でありp+1p+1以上でないとき、次数はppであるという。
  5. 次数がppでありσ(1)≠0\sigma(1)\ne0であるとき、Cp+1/σ(1)C_{p+1}/\sigma(1)を線形多段法の 誤差定数 (error constant) という。

補題 1.3.kk段の線形多段法を取り、C0,C1,…C_0,C_1,\dotsを定義 1.2 (2)の係数とする。

  1. q∈N≥1q\in\NNとし、 Kq:=∑j=0k(∣αj∣jq+1(q+1)!+∣βj∣jqq!)K_q:=\sum_{j=0}^k\Bigl(\frac{|\alpha_j|j^{q+1}}{(q+1)!}+\frac{|\beta_j|j^q}{q!}\Bigr) と置く。d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。t0<Tt_0<Tとし、J⊇[t0,T]J\supseteq[t_0,T]を開区間、y∈Cq+1(J;Rd)y\in C^{q+1}(J;\R^d)、M:=max⁡τ∈[t0,T]∥y(q+1)(τ)∥M:=\max_{\tau\in[t_0,T]}\|y^{(q+1)}(\tau)\|とする。このとき、0<h≤(T−t0)/k0<h\le(T-t_0)/kとt∈[t0,T−kh]t\in[t_0,T-kh]を満たす任意のh,th,tについて ∥dh(t)−∑i=0qCihiy(i)(t)∥≤KqMhq+1\Bigl\|d_h(t)-\sum_{i=0}^qC_ih^iy^{(i)}(t)\Bigr\|\le K_qMh^{q+1} が成り立つ。
  2. i∈N≥0i\in\N、s∈Rs\in\Rとし、y(τ):=(τ−s)iy(\tau):=(\tau-s)^i(τ∈R\tau\in\R)とする。任意のh>0h>0について∑j=0kαjy(s+jh)−h∑j=0kβjy′(s+jh)=i! Cihi\sum_{j=0}^k\alpha_jy(s+jh)-h\sum_{j=0}^k\beta_jy'(s+jh)=i!\,C_ih^iである。

証明.(1)を示す。0≤j≤k0\le j\le kとすると[t,t+jh]⊆[t0,T][t,t+jh]\subseteq[t_0,T]である。§E20.27 補題 1.5 (1)を、yyとr=qr=q、展開点tt、増分jhjhに適用し、y′y'とr=q−1r=q-1、同じ展開点と増分に適用すると

∥y(t+jh)−∑i=0q(jh)ii!y(i)(t)∥≤(jh)q+1(q+1)!M,∥y′(t+jh)−∑i=1q(jh)i−1(i−1)!y(i)(t)∥≤(jh)qq!M\Bigl\|y(t+jh)-\sum_{i=0}^q\frac{(jh)^i}{i!}y^{(i)}(t)\Bigr\|\le\frac{(jh)^{q+1}}{(q+1)!}M,\qquad \Bigl\|y'(t+jh)-\sum_{i=1}^q\frac{(jh)^{i-1}}{(i-1)!}y^{(i)}(t)\Bigr\|\le\frac{(jh)^q}{q!}M

である。CiC_iの定義により

∑j=0kαj∑i=0q(jh)ii!y(i)(t)−h∑j=0kβj∑i=1q(jh)i−1(i−1)!y(i)(t)=∑i=0qCihiy(i)(t)\sum_{j=0}^k\alpha_j\sum_{i=0}^q\frac{(jh)^i}{i!}y^{(i)}(t)-h\sum_{j=0}^k\beta_j\sum_{i=1}^q\frac{(jh)^{i-1}}{(i-1)!}y^{(i)}(t)=\sum_{i=0}^qC_ih^iy^{(i)}(t)

であるから、第一の評価に∣αj∣|\alpha_j|を、第二の評価にh∣βj∣h|\beta_j|を掛けてjjについて加えると主張を得る。

(2)を示す。y(s+jh)=jihiy(s+jh)=j^ih^iであり、i≥1i\ge1ならばy′(s+jh)=iji−1hi−1y'(s+jh)=ij^{i-1}h^{i-1}、i=0i=0ならばy′=0y'=0である。これを代入すると左辺はi=0i=0でC0C_0、i≥1i\ge1でhi(∑jαjji−i∑jβjji−1)=i! Cihih^i\bigl(\sum_j\alpha_jj^i-i\sum_j\beta_jj^{i-1}\bigr)=i!\,C_ih^iである。▨

命題 1.4.kk段の線形多段法を取り、ρ,σ\rho,\sigmaをその特性多項式、C0,C1,…C_0,C_1,\dotsを定義 1.2 (2)の係数とする。C0=ρ(1)C_0=\rho(1)、C1=ρ′(1)−σ(1)C_1=\rho'(1)-\sigma(1)である。

  1. 線形多段法が整合的であることと、次数が11以上であることと、ρ(1)=0\rho(1)=0かつρ′(1)=σ(1)\rho'(1)=\sigma(1)であることは同値である。
  2. p∈N≥1p\in\NNとする。次数がpp以上であることとC0=⋯=Cp=0C_0=\dots=C_p=0は同値である。次数がppであることと、C0=⋯=Cp=0C_0=\dots=C_p=0かつCp+1≠0C_{p+1}\ne0であることは同値である。
  3. p∈N≥1p\in\NNとし、次数がpp以上であるとする。d∈N≥1d\in\NN、Rd\R^dのノルム∥⋅∥\|\cdot\|、t0<Tt_0<T、開区間J⊇[t0,T]J\supseteq[t_0,T]を取り、y∈Cp+2(J;Rd)y\in C^{p+2}(J;\R^d)、M:=max⁡τ∈[t0,T]∥y(p+2)(τ)∥M:=\max_{\tau\in[t_0,T]}\|y^{(p+2)}(\tau)\|とする。0<h≤(T−t0)/k0<h\le(T-t_0)/kとt∈[t0,T−kh]t\in[t_0,T-kh]を満たす任意のh,th,tについて ∥dh(t)−Cp+1hp+1y(p+1)(t)∥≤Kp+1Mhp+2\bigl\|d_h(t)-C_{p+1}h^{p+1}y^{(p+1)}(t)\bigr\|\le K_{p+1}Mh^{p+2} が成り立つ。ここでKp+1K_{p+1}は補題 1.3 (1)の定数である。

証明.C0=ρ(1)C_0=\rho(1)は定義であり、C1=∑jαjj−∑jβj=ρ′(1)−σ(1)C_1=\sum_j\alpha_jj-\sum_j\beta_j=\rho'(1)-\sigma(1)である。i∈N≥0i\in\Nとし、d=1d=1、t0=0t_0=0、T=kT=k、J=RJ=\R、y(τ):=τiy(\tau):=\tau^iとすると、補題 1.3 (2)により0<h≤10<h\le1についてdh(0)=i! Cihid_h(0)=i!\,C_ih^iである。

(2)を示す。C0=⋯=Cp=0C_0=\dots=C_p=0ならば、補題 1.3 (1)をq=pq=pとして適用すると∥dh(t)∥≤KpMhp+1\|d_h(t)\|\le K_pMh^{p+1}であり、次数はpp以上である。逆に次数がpp以上であるとし、0≤i≤p0\le i\le pとする。上のy(τ)=τiy(\tau)=\tau^iに次数の定義のCCを取ると、0<h≤10<h\le1についてi! ∣Ci∣hi≤Chp+1i!\,|C_i|h^i\le Ch^{p+1}であり、∣Ci∣≤Chp+1−i/i!|C_i|\le Ch^{p+1-i}/i!の右辺はh→+0h\to+0で00に収束するからCi=0C_i=0である。次数がppであることについての主張は、この同値をppとp+1p+1に用いると従う。

(1)を示す。(2)により、次数が11以上であることとC0=C1=0C_0=C_1=0は同値であり、これはρ(1)=0\rho(1)=0かつρ′(1)=σ(1)\rho'(1)=\sigma(1)と同値である。

C0=C1=0C_0=C_1=0とし、d∈N≥1d\in\NN、Rd\R^dのノルム∥⋅∥\|\cdot\|、t0<Tt_0<T、開区間J⊇[t0,T]J\supseteq[t_0,T]、y∈C1(J;Rd)y\in C^1(J;\R^d)を取る。A:=∑j=0kj∣αj∣+∑j=0k∣βj∣A:=\sum_{j=0}^kj|\alpha_j|+\sum_{j=0}^k|\beta_j|と置く。[t0,T][t_0,T]は§E2.9 定理 4.3によりコンパクトであるから、§E2.9 定理 5.1によりy′y'は[t0,T][t_0,T]上で一様連続である。ε>0\varepsilon>0を取り、s,s′∈[t0,T]s,s'\in[t_0,T]が∣s−s′∣<δ|s-s'|<\deltaを満たすならば∥y′(s)−y′(s′)∥<ε\|y'(s)-y'(s')\|<\varepsilonであるδ>0\delta>0を取る。0<h≤(T−t0)/k0<h\le(T-t_0)/k、kh<δkh<\delta、t∈[t0,T−kh]t\in[t_0,T-kh]とし、0≤j≤k0\le j\le kとすると[t,t+jh]⊆[t0,T][t,t+jh]\subseteq[t_0,T]である。x(τ):=y(τ)−τy′(t)x(\tau):=y(\tau)-\tau y'(t)(τ∈J\tau\in J)はC1C^1級でありx′(τ)=y′(τ)−y′(t)x'(\tau)=y'(\tau)-y'(t)であるから、§E20.27 補題 1.5 (1)をxxとr=0r=0、展開点tt、増分jhjhに適用すると∥y(t+jh)−y(t)−jhy′(t)∥≤jhε\|y(t+jh)-y(t)-jhy'(t)\|\le jh\varepsilonである。また∥y′(t+jh)−y′(t)∥<ε\|y'(t+jh)-y'(t)\|<\varepsilonである。C0=C1=0C_0=C_1=0により

dh(t)=∑j=0kαj(y(t+jh)−y(t)−jhy′(t))−h∑j=0kβj(y′(t+jh)−y′(t))d_h(t)=\sum_{j=0}^k\alpha_j\bigl(y(t+jh)-y(t)-jhy'(t)\bigr)-h\sum_{j=0}^k\beta_j\bigl(y'(t+jh)-y'(t)\bigr)

であるから、∥dh(t)∥≤Aεh\|d_h(t)\|\le A\varepsilon hである。したがって0<h≤(T−t0)/k0<h\le(T-t_0)/kかつkh<δkh<\deltaならばsup⁡t∥dh(t)∥/h≤Aε\sup_t\|d_h(t)\|/h\le A\varepsilonであり、線形多段法は整合的である。

整合的であるとする。y(τ)=1y(\tau)=1についてはdh(0)=C0d_h(0)=C_0であるからsup⁡t∣dh(t)∣/h≥∣C0∣/h\sup_t|d_h(t)|/h\ge|C_0|/hであり、左辺はh→+0h\to+0で00に収束するのでC0=0C_0=0である。y(τ)=τy(\tau)=\tauについてはdh(0)=C1hd_h(0)=C_1hであるからsup⁡t∣dh(t)∣/h≥∣C1∣\sup_t|d_h(t)|/h\ge|C_1|であり、C1=0C_1=0である。

(3)は、(2)によりC0=⋯=Cp=0C_0=\dots=C_p=0であるから、補題 1.3 (1)をq=p+1q=p+1として適用すると従う。▨

注意 1.5.θ∈[0,1]\theta\in[0,1]とし、k=1k=1、ρ(ζ)=ζ−1\rho(\zeta)=\zeta-1、σ(ζ)=θζ+1−θ\sigma(\zeta)=\theta\zeta+1-\thetaの線形多段法を考える。近似値を定める方程式は θ 法の方程式(§E20.28 定義 2.1)と同じであり、その解がただ一つである限りyn+1=Ψh(tn,yn)y_{n+1}=\Psi_h(t_n,y_n)である。係数はC0=C1=0C_0=C_1=0、C2=12−θC_2=\frac12-\thetaであるから、命題 1.4 (2)により次数はθ≠12\theta\ne\frac12ならば11、θ=12\theta=\frac12ならば22以上である。局所欠陥dh(t)d_h(t)は§E20.28 補題 2.2 (2)のη\etaに等しく、同じ補題の仮定の下で一段法の局所打切り誤差δy(t,h)=y(t+h)−Ψh(t,y(t))\delta_y(t,h)=y(t+h)-\Psi_h(t,y(t))との差は∥δy(t,h)−dh(t)∥≤θhL∥dh(t)∥/(1−θhL)\|\delta_y(t,h)-d_h(t)\|\le\theta hL\|d_h(t)\|/(1-\theta hL)を満たす。

注意 1.6.kk段の線形多段法(ρ,σ)(\rho,\sigma)とa∈R∖{1}a\in\R\setminus\{1\}に対して、ρ^(ζ):=(ζ−a)ρ(ζ)\hat\rho(\zeta):=(\zeta-a)\rho(\zeta)、σ^(ζ):=(ζ−a)σ(ζ)\hat\sigma(\zeta):=(\zeta-a)\sigma(\zeta)を特性多項式とするk+1k+1段の線形多段法を考える。ρ^\hat\rhoの最高次の係数は11である。α−1=αk+1=β−1=βk+1:=0\alpha_{-1}=\alpha_{k+1}=\beta_{-1}=\beta_{k+1}:=0として係数をα^j=αj−1−aαj\hat\alpha_j=\alpha_{j-1}-a\alpha_j、β^j=βj−1−aβj\hat\beta_j=\beta_{j-1}-a\beta_jと書いて比べると、任意のC1C^1級写像yyについて局所欠陥はd^h(t)=dh(t+h)−adh(t)\hat d_h(t)=d_h(t+h)-ad_h(t)を満たす。

(ρ,σ)(\rho,\sigma)の次数がppでありσ(1)≠0\sigma(1)\ne0であるとする。i∈N≥0i\in\N、i≤p+1i\le p+1、s∈Rs\in\Rとし、y(τ):=(τ−s)iy(\tau):=(\tau-s)^iとする。q:=max⁡{i,1}q:=\max\{i,1\}とするとy(q+1)=0y^{(q+1)}=0であるから、補題 1.3 (1)によりdh(t)=∑l=0qClhly(l)(t)d_h(t)=\sum_{l=0}^qC_lh^ly^{(l)}(t)である。命題 1.4 (2)によりC0=⋯=Cp=0C_0=\dots=C_p=0であるから、i≤pi\le pならばdh≡0d_h\equiv0であり、i=p+1i=p+1ならばdh(t)=(p+1)! Cp+1hp+1d_h(t)=(p+1)!\,C_{p+1}h^{p+1}である。補題 1.3 (2)を(ρ^,σ^)(\hat\rho,\hat\sigma)に適用するとi! C^ihi=d^h(s)=dh(s+h)−adh(s)i!\,\hat C_ih^i=\hat d_h(s)=d_h(s+h)-ad_h(s)であるから、i≤pi\le pでC^i=0\hat C_i=0、i=p+1i=p+1でC^p+1=(1−a)Cp+1\hat C_{p+1}=(1-a)C_{p+1}である。

(ρ,σ)(\rho,\sigma)の次数がppであるから命題 1.4 (2)によりCp+1≠0C_{p+1}\ne0であり、a≠1a\ne1によりC^p+1≠0\hat C_{p+1}\ne0である。したがって(ρ^,σ^)(\hat\rho,\hat\sigma)の次数もppである。σ^(1)=(1−a)σ(1)≠0\hat\sigma(1)=(1-a)\sigma(1)\ne0であり、C^p+1/σ^(1)=Cp+1/σ(1)\hat C_{p+1}/\hat\sigma(1)=C_{p+1}/\sigma(1)であるから、二つの線形多段法の誤差定数は等しい。

2 ゼロ安定性と収束

定義 2.1.ppを次数11以上の複素係数多項式とする。ppの複素数の根がすべて∣ζ∣≤1|\zeta|\le1を満たし、∣ζ∣=1|\zeta|=1を満たす根がすべて単根であるとき、ppは 根条件 (root condition) を満たすという。線形多段法は、第一特性多項式ρ\rhoが根条件を満たすとき ゼロ安定 (zero-stable) であるという。

補題 2.2.k∈N≥1k\in\NNとし、a0,…,ak−1∈Ca_0,\dots,a_{k-1}\in\C、ak:=1a_k:=1、p(ζ):=∑j=0kajζjp(\zeta):=\sum_{j=0}^ka_j\zeta^jとする。任意のn∈N≥0n\in\Nについて∑j=0kajun+j=0\sum_{j=0}^ka_ju_{n+j}=0を満たす複素数列(un)n∈N≥0(u_n)_{n\in\N}を差分方程式の解という。0≤j<k0\le j<kについて、0≤i<k0\le i<kでui=1u_i=1(i=ji=j)、ui=0u_i=0(i≠ji\ne j)を満たす解をγ(j)=(γn(j))n∈N≥0\gamma^{(j)}=(\gamma^{(j)}_n)_{n\in\N}とする。次の二条件は同値である。

  1. 差分方程式の任意の解は有界である。
  2. ppは根条件を満たす。

さらに次が成り立つ。

  1. 上の二条件が成り立つとき、Γ:=sup⁡n∈N≥0∑j=0k−1∣γn(j)∣\Gamma:=\sup_{n\in\N}\sum_{j=0}^{k-1}|\gamma^{(j)}_n|は有限かつ11以上であり、差分方程式の任意の解(un)(u_n)と任意のn∈N≥0n\in\Nについて∣un∣≤Γmax⁡0≤j<k∣uj∣|u_n|\le\Gamma\max_{0\le j<k}|u_j|である。a0,…,ak−1a_0,\dots,a_{k-1}が実数ならば、任意のd∈N≥1d\in\NN、Rd\R^dの任意のノルム∥⋅∥\|\cdot\|と、任意のn∈N≥0n\in\Nについて∑j=0kajun+j=0\sum_{j=0}^ka_ju_{n+j}=0を満たすRd\R^dの任意の列(un)(u_n)について、∥un∥≤Γmax⁡0≤j<k∥uj∥\|u_n\|\le\Gamma\max_{0\le j<k}\|u_j\|である。
  2. 差分方程式の任意の解が00に収束することと、ppの根がすべて∣ζ∣<1|\zeta|<1を満たすことは同値である。

証明.A∈Ck×kA\in\C^{k\times k}を、U=(U0,…,Uk−1)TU=(U_0,\dots,U_{k-1})^{\mathsf T}に対して(AU)i=Ui+1(AU)_i=U_{i+1}(0≤i≤k−20\le i\le k-2)、(AU)k−1=−∑j=0k−1ajUj(AU)_{k-1}=-\sum_{j=0}^{k-1}a_jU_jとなる行列とする。複素数列(un)(u_n)が解であることと、Un:=(un,…,un+k−1)TU_n:=(u_n,\dots,u_{n+k-1})^{\mathsf T}が任意のnnについてUn+1=AUnU_{n+1}=AU_nを満たすことは同値であり、このときUn=AnU0U_n=A^nU_0である。解はU0U_0によってただ一つに定まり、∑j<kujγ(j)\sum_{j<k}u_j\gamma^{(j)}はU0U_0から定まる解であるから

un=∑j=0k−1ujγn(j)(n∈N≥0)u_n=\sum_{j=0}^{k-1}u_j\gamma^{(j)}_n\qquad(n\in\N)

である。

λ\lambdaをAAの固有値とし、v≠0v\ne0をAv=λvAv=\lambda vを満たすベクトルとする。0≤i≤k−20\le i\le k-2の成分からvi+1=λviv_{i+1}=\lambda v_iであり、vi=λiv0v_i=\lambda^iv_0、v0≠0v_0\ne0である。第k−1k-1成分から−∑j<kajλjv0=λkv0-\sum_{j<k}a_j\lambda^jv_0=\lambda^kv_0であり、p(λ)=0p(\lambda)=0である。さらにv0=1v_0=1とし、w∈Ckw\in\C^kが(A−λI)w=v(A-\lambda I)w=vを満たすとする。0≤i≤k−20\le i\le k-2の成分からwi+1=λwi+λiw_{i+1}=\lambda w_i+\lambda^iであり、iiについての帰納法によりwi=λiw0+iλi−1w_i=\lambda^iw_0+i\lambda^{i-1}(i=0i=0では第二項を00と読む)である。これを第k−1k-1成分の等式−∑j<kajwj−λwk−1=λk−1-\sum_{j<k}a_jw_j-\lambda w_{k-1}=\lambda^{k-1}に代入するとw0p(λ)+p′(λ)=0w_0p(\lambda)+p'(\lambda)=0を得る。p(λ)=0p(\lambda)=0であるからp′(λ)=0p'(\lambda)=0であり、λ\lambdaはppの重根である。

(2)⇒\Rightarrow(1)を示す。複素係数の多項式は一次式の積に分解するから、Ck\C^kの線形変換U↦AUU\mapsto AUの特性多項式は一次式の積に分解する。§E3.30 系 1.2により、Ck\C^kの基底B\mathcal Bが存在して[A]B=Jm1(λ1)⊕⋯⊕Jmr(λr)[A]_{\mathcal B}=J_{m_1}(\lambda_1)\oplus\dots\oplus J_{m_r}(\lambda_r)、[An]B=Jm1(λ1)n⊕⋯⊕Jmr(λr)n[A^n]_{\mathcal B}=J_{m_1}(\lambda_1)^n\oplus\dots\oplus J_{m_r}(\lambda_r)^nである。各λl\lambda_lはAAの固有値であるからppの根であり、根条件により∣λl∣≤1|\lambda_l|\le1である。

ml≥2m_l\ge2とし、Jml(λl)J_{m_l}(\lambda_l)に対応するB\mathcal Bの元が張る部分空間をWWとする。WWはAAで不変であり、A−λlIA-\lambda_lIのWWへの制限の表現行列Jml(λl)−λlIJ_{m_l}(\lambda_l)-\lambda_lIは零でない冪零行列である。(A−λlI)sx≠0(A-\lambda_lI)^sx\ne0かつ(A−λlI)s+1x=0(A-\lambda_lI)^{s+1}x=0を満たすx∈Wx\in Wとs≥1s\ge1を取り、b:=(A−λlI)sxb:=(A-\lambda_lI)^sx、b′:=(A−λlI)s−1xb':=(A-\lambda_lI)^{s-1}xと置くと、b≠0b\ne0、Ab=λlbAb=\lambda_lb、(A−λlI)b′=b(A-\lambda_lI)b'=bである。AAの固有ベクトルの第00成分は00でないからb0≠0b_0\ne0であり、v:=b/b0v:=b/b_0、w:=b′/b0w:=b'/b_0はv0=1v_0=1、Av=λlvAv=\lambda_lv、(A−λlI)w=v(A-\lambda_lI)w=vを満たすので、λl\lambda_lはppの重根である。したがって根条件により∣λl∣=1|\lambda_l|=1ならばml=1m_l=1である。

§E3.30 命題 1.1によりJm(λ)nJ_{m}(\lambda)^nの成分は00または(ni)λn−i\binom ni\lambda^{n-i}(0≤i≤min⁡{n,m−1}0\le i\le\min\{n,m-1\})である。ml=1m_l=1、∣λl∣=1|\lambda_l|=1ならば∣λln∣=1|\lambda_l^n|=1である。∣λl∣<1|\lambda_l|<1ならば、λl=0\lambda_l=0のときn>in>iで(ni)λln−i=0\binom ni\lambda_l^{n-i}=0であり、λl≠0\lambda_l\ne0のとき∣(ni)λln−i∣≤ni∣λl∣n−i|\binom ni\lambda_l^{n-i}|\le n^i|\lambda_l|^{n-i}であって、0<r<10<r<1についてnirn→0n^ir^n\to0であるから、いずれの場合も(ni)λln−i→0\binom ni\lambda_l^{n-i}\to0(n→∞n\to\infty)である。したがって[An]B[A^n]_{\mathcal B}の成分はnnについて有界であり、Un=AnU0U_n=A^nU_0のB\mathcal Bに関する座標は有界であるから、(un)(u_n)は有界である。

(1)⇒\Rightarrow(2)を示す。ppが根条件を満たさないとすると、∣ζ∣>1|\zeta|>1を満たす根ζ\zetaが存在するか、∣ζ∣=1|\zeta|=1かつp(ζ)=p′(ζ)=0p(\zeta)=p'(\zeta)=0を満たすζ\zetaが存在する。前者の場合un:=ζnu_n:=\zeta^nと置くと∑jajun+j=ζnp(ζ)=0\sum_ja_ju_{n+j}=\zeta^np(\zeta)=0であり、∣un∣=∣ζ∣n|u_n|=|\zeta|^nは有界でない。後者の場合un:=nζnu_n:=n\zeta^nと置くと

∑j=0kajun+j=nζnp(ζ)+ζn+1p′(ζ)=0\sum_{j=0}^ka_ju_{n+j}=n\zeta^np(\zeta)+\zeta^{n+1}p'(\zeta)=0

であり、∣un∣=n|u_n|=nは有界でない。

(1)を示す。(1)により各γ(j)\gamma^{(j)}は有界であるからΓ≤∑jsup⁡n∣γn(j)∣<∞\Gamma\le\sum_j\sup_n|\gamma^{(j)}_n|<\inftyであり、n=0n=0の項は∣γ0(0)∣=1|\gamma^{(0)}_0|=1であるからΓ≥1\Gamma\ge1である。上の表示により∣un∣≤∑j∣uj∣∣γn(j)∣≤Γmax⁡j∣uj∣|u_n|\le\sum_j|u_j||\gamma^{(j)}_n|\le\Gamma\max_j|u_j|である。a0,…,ak−1a_0,\dots,a_{k-1}が実数ならば、漸化式un+k=−∑j<kajun+ju_{n+k}=-\sum_{j<k}a_ju_{n+j}によりγ(j)\gamma^{(j)}は実数列である。Rd\R^dの列(un)(u_n)について、∑j<kγn(j)uj\sum_{j<k}\gamma^{(j)}_nu_jは同じ漸化式を満たし、n<kn<kでunu_nに一致するから、任意のnnでunu_nに等しい。したがって∥un∥≤∑j∣γn(j)∣∥uj∥≤Γmax⁡j∥uj∥\|u_n\|\le\sum_j|\gamma^{(j)}_n|\|u_j\|\le\Gamma\max_j\|u_j\|である。

(2)を示す。ppの根がすべて∣ζ∣<1|\zeta|<1を満たすならば、各λl\lambda_lはppの根であるから∣λl∣<1|\lambda_l|<1であり、上の評価により[An]B[A^n]_{\mathcal B}の成分はすべて00に収束する。したがってUn→0U_n\to0でありun→0u_n\to0である。∣ζ∣≥1|\zeta|\ge1を満たす根ζ\zetaが存在するならば、解un=ζnu_n=\zeta^nは∣un∣≥1|u_n|\ge1を満たし、00に収束しない。▨

補題 2.3.k∈N≥1k\in\NNとし、実数a0,…,ak−1a_0,\dots,a_{k-1}とak:=1a_k:=1についてp(ζ):=∑j=0kajζjp(\zeta):=\sum_{j=0}^ka_j\zeta^jが根条件を満たすとする。Γ\Gammaを補題 2.2 (1)の定数とする。d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。整数N≥kN\ge kとu0,…,uN∈Rdu_0,\dots,u_N\in\R^d、φ0,…,φN−k∈Rd\varphi_0,\dots,\varphi_{N-k}\in\R^dが、0≤m≤N−k0\le m\le N-kについて∑j=0kajum+j=φm\sum_{j=0}^ka_ju_{m+j}=\varphi_mを満たすとする。このとき0≤n≤N0\le n\le Nについて

∥un∥≤Γ(max⁡0≤j<k∥uj∥+∑m=0n−k∥φm∥)\|u_n\|\le\Gamma\Bigl(\max_{0\le j<k}\|u_j\|+\sum_{m=0}^{n-k}\|\varphi_m\|\Bigr)

が成り立つ。n<kn<kのとき和は00と読む。

証明.γ(k−1)\gamma^{(k-1)}を補題 2.2の解とし、i∈N≥0i\in\Nについてθi:=γi(k−1)\theta_i:=\gamma^{(k-1)}_i、負の整数iiについてθi:=0\theta_i:=0と置く。θ0=⋯=θk−2=0\theta_0=\dots=\theta_{k-2}=0、θk−1=1\theta_{k-1}=1である。整数i≥−ki\ge-kについて

∑j=0kajθi+j={1(i=−1),0(i≠−1)\sum_{j=0}^ka_j\theta_{i+j}=\begin{cases}1&(i=-1),\\0&(i\ne-1)\end{cases}

である。実際、i≥0i\ge0ではθ\thetaが差分方程式の解であることによる。−k≤i≤−1-k\le i\le-1ではi+j≤k−1i+j\le k-1であり、θi+j≠0\theta_{i+j}\ne0となりうるのはi+j=k−1i+j=k-1のときだけであるから、j≤kj\le kによりi=−1i=-1、j=kj=kに限られ、その項はakθk−1=1a_k\theta_{k-1}=1である。

0≤n≤N0\le n\le Nについてwn:=∑m=0n−kθn−m−1φmw_n:=\sum_{m=0}^{n-k}\theta_{n-m-1}\varphi_mと置く。m≥n−k+1m\ge n-k+1ならばn−m−1≤k−2n-m-1\le k-2でありθn−m−1=0\theta_{n-m-1}=0であるから、0≤l≤N−k0\le l\le N-kと0≤j≤k0\le j\le kについてwl+j=∑m=0lθl+j−m−1φmw_{l+j}=\sum_{m=0}^{l}\theta_{l+j-m-1}\varphi_mである。したがって

∑j=0kajwl+j=∑m=0lφm∑j=0kajθ(l−m−1)+j=φl\sum_{j=0}^ka_jw_{l+j}=\sum_{m=0}^l\varphi_m\sum_{j=0}^ka_j\theta_{(l-m-1)+j}=\varphi_l

であり、n<kn<kでwn=0w_n=0である。vn:=un−wnv_n:=u_n-w_nは0≤l≤N−k0\le l\le N-kで∑jajvl+j=0\sum_ja_jv_{l+j}=0を満たし、n<kn<kでvn=unv_n=u_nである。vN+1,vN+2,…v_{N+1},v_{N+2},\dotsを漸化式vn+k=−∑j<kajvn+jv_{n+k}=-\sum_{j<k}a_jv_{n+j}で定めると、補題 2.2 (1)により∥vn∥≤Γmax⁡j<k∥uj∥\|v_n\|\le\Gamma\max_{j<k}\|u_j\|である。Γ\Gammaの定義により∣θi∣≤Γ|\theta_i|\le\Gammaであるから∥wn∥≤Γ∑m=0n−k∥φm∥\|w_n\|\le\Gamma\sum_{m=0}^{n-k}\|\varphi_m\|であり、un=vn+wnu_n=v_n+w_nから主張を得る。▨

定義 2.4.kk段の線形多段法が次を満たすとき、線形多段法は 収束する (convergent) という。d∈N≥1d\in\NN、Rd\R^dのノルム∥⋅∥\|\cdot\|、開区間I⊆RI\subseteq\R、連続写像f ⁣:I×Rd→Rdf\colon I\times\R^d\to\R^dとL≥0L\ge0が、任意のs∈Is\in Iとu,v∈Rdu,v\in\R^dについて∥f(s,u)−f(s,v)∥≤L∥u−v∥\|f(s,u)-f(s,v)\|\le L\|u-v\|を満たすとし、t0<Tt_0<T、[t0,T]⊆J⊆I[t_0,T]\subseteq J\subseteq Iを満たす開区間JJ、任意のτ∈J\tau\in Jについてy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たすC1C^1級写像y ⁣:J→Rdy\colon J\to\R^dを取る。このようなすべての組について、整数N1≥kN_1\ge kが存在して次が成り立つ。N≥N1N\ge N_1についてhN:=(T−t0)/Nh_N:=(T-t_0)/N、tn:=t0+nhNt_n:=t_0+nh_Nと置くと、任意の開始値y0(N),…,yk−1(N)∈Rdy^{(N)}_0,\dots,y^{(N)}_{k-1}\in\R^dから、刻みhNh_Nによる定義 1.1 (2)の近似値y0(N),…,yN(N)y^{(N)}_0,\dots,y^{(N)}_Nがすべて定まる。さらに、各N≥N1N\ge N_1について開始値を取り

εN:=max⁡0≤j<k∥yj(N)−y(tj)∥,EN:=max⁡0≤n≤N∥yn(N)−y(tn)∥\varepsilon_N:=\max_{0\le j<k}\|y^{(N)}_j-y(t_j)\|,\qquad E_N:=\max_{0\le n\le N}\|y^{(N)}_n-y(t_n)\|

と置くと、εN→0\varepsilon_N\to0(N→∞N\to\infty)ならばEN→0E_N\to0である。

定理 2.5.kk段の線形多段法(α0,…,αk;β0,…,βk)(\alpha_0,\dots,\alpha_k;\beta_0,\dots,\beta_k)がゼロ安定であるとし、Γ\Gammaをaj:=αja_j:=\alpha_jとした補題 2.2 (1)の定数、βˉ:=∑j=0k∣βj∣\bar\beta:=\sum_{j=0}^k|\beta_j|とする。d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。I⊆RI\subseteq\Rを開区間、f ⁣:I×Rd→Rdf\colon I\times\R^d\to\R^dを連続写像とし、L≥0L\ge0が、任意のs∈Is\in Iとu,v∈Rdu,v\in\R^dについて∥f(s,u)−f(s,v)∥≤L∥u−v∥\|f(s,u)-f(s,v)\|\le L\|u-v\|を満たすとする。t0<Tt_0<Tとし、JJを[t0,T]⊆J⊆I[t_0,T]\subseteq J\subseteq Iを満たす開区間、y ⁣:J→Rdy\colon J\to\R^dを、任意のτ∈J\tau\in Jについてy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たすC1C^1級写像とする。整数N≥kN\ge kについてh:=(T−t0)/Nh:=(T-t_0)/N、tn:=t0+nht_n:=t_0+nhと置き、2ΓhL∣βk∣≤12\Gamma hL|\beta_k|\le1とする。このとき任意の開始値y0,…,yk−1∈Rdy_0,\dots,y_{k-1}\in\R^dから近似値y0,…,yNy_0,\dots,y_Nがすべて定まり、en:=yn−y(tn)e_n:=y_n-y(t_n)、dm:=dh(tm)d_m:=d_h(t_m)(0≤m≤N−k0\le m\le N-k)と置くと

max⁡0≤n≤N∥en∥≤2Γe2ΓLβˉ(T−t0)(max⁡0≤j<k∥ej∥+∑m=0N−k∥dm∥)\max_{0\le n\le N}\|e_n\|\le2\Gamma e^{2\Gamma L\bar\beta(T-t_0)}\Bigl(\max_{0\le j<k}\|e_j\|+\sum_{m=0}^{N-k}\|d_m\|\Bigr)

が成り立つ。

証明.0≤n≤N−k0\le n\le N-kとし、yn,…,yn+k−1y_n,\dots,y_{n+k-1}が定まっているとする。c:=−∑j<kαjyn+j+h∑j<kβjf(tn+j,yn+j)c:=-\sum_{j<k}\alpha_jy_{n+j}+h\sum_{j<k}\beta_jf(t_{n+j},y_{n+j})、Φ(v):=c+hβkf(tn+k,v)\Phi(v):=c+h\beta_kf(t_{n+k},v)と置くと、近似値の方程式の解はΦ\Phiの不動点と一致する。∥Φ(v)−Φ(v′)∥≤h∣βk∣L∥v−v′∥\|\Phi(v)-\Phi(v')\|\le h|\beta_k|L\|v-v'\|であり、Γ≥1\Gamma\ge1によりh∣βk∣L≤1/2h|\beta_k|L\le1/2である。Rd\R^dは空でない完備距離空間であるから、§E2.7 定理 2.2によりΦ\Phiはただ一つの不動点をもち、yn+ky_{n+k}が定まる。nnについての帰納法によりy0,…,yNy_0,\dots,y_Nはすべて定まる。

gi:=f(ti,yi)−f(ti,y(ti))g_i:=f(t_i,y_i)-f(t_i,y(t_i))と置くと∥gi∥≤L∥ei∥\|g_i\|\le L\|e_i\|である。0≤m≤N−k0\le m\le N-kについて、近似値の方程式と局所欠陥の定義∑jαjy(tm+j)=h∑jβjf(tm+j,y(tm+j))+dm\sum_j\alpha_jy(t_{m+j})=h\sum_j\beta_jf(t_{m+j},y(t_{m+j}))+d_mの差を取ると

∑j=0kαjem+j=h∑j=0kβjgm+j−dm\sum_{j=0}^k\alpha_je_{m+j}=h\sum_{j=0}^k\beta_jg_{m+j}-d_m

である。Es:=max⁡j<k∥ej∥E_{\mathrm s}:=\max_{j<k}\|e_j\|、D:=∑m=0N−k∥dm∥D:=\sum_{m=0}^{N-k}\|d_m\|と置き、k≤n≤Nk\le n\le Nとする。補題 2.3をun:=enu_n:=e_n、右辺をφm\varphi_mとして適用すると

∥en∥≤Γ(Es+D+hL∑m=0n−k∑j=0k∣βj∣∥em+j∥)\|e_n\|\le\Gamma\Bigl(E_{\mathrm s}+D+hL\sum_{m=0}^{n-k}\sum_{j=0}^k|\beta_j|\|e_{m+j}\|\Bigr)

である。二重和のうち∥en∥\|e_n\|を含む項はm=n−km=n-k、j=kj=kの項∣βk∣∥en∥|\beta_k|\|e_n\|だけであり、i<ni<nについて∥ei∥\|e_i\|に掛かる係数の和はβˉ\bar\beta以下である。ΓhL∣βk∣≤1/2\Gamma hL|\beta_k|\le1/2によりΓhL∣βk∣∥en∥\Gamma hL|\beta_k|\|e_n\|を左辺へ移すと

∥en∥≤2Γ(Es+D)+Λh∑i=0n−1∥ei∥,Λ:=2ΓLβˉ\|e_n\|\le2\Gamma(E_{\mathrm s}+D)+\Lambda h\sum_{i=0}^{n-1}\|e_i\|,\qquad\Lambda:=2\Gamma L\bar\beta

である。n<kn<kでは∥en∥≤Es≤2Γ(Es+D)\|e_n\|\le E_{\mathrm s}\le2\Gamma(E_{\mathrm s}+D)であるから、c′:=2Γ(Es+D)c':=2\Gamma(E_{\mathrm s}+D)、Sn:=∑i<n∥ei∥S_n:=\sum_{i<n}\|e_i\|と置くと、0≤n≤N0\le n\le Nについて∥en∥≤c′+ΛhSn\|e_n\|\le c'+\Lambda hS_nである。

0≤n<N0\le n<NについてSn+1=Sn+∥en∥≤(1+Λh)Sn+c′S_{n+1}=S_n+\|e_n\|\le(1+\Lambda h)S_n+c'である。§E20.28 補題 1.2 (1)を格子t0<⋯<tNt_0<\dots<t_N、En:=SnE_n:=S_n(E0=0E_0=0)、an:=c′a_n:=c'として適用し、§E20.28 補題 1.2 (3)をその主張のα\alphaをc′/hc'/hとして適用すると、Sn≤(c′/h)φΛ(tn−t0)S_n\le(c'/h)\varphi_\Lambda(t_n-t_0)である。1+ΛφΛ(τ)=eΛτ1+\Lambda\varphi_\Lambda(\tau)=e^{\Lambda\tau}であるから

∥en∥≤c′(1+ΛφΛ(tn−t0))=c′eΛ(tn−t0)≤2Γe2ΓLβˉ(T−t0)(Es+D)\|e_n\|\le c'\bigl(1+\Lambda\varphi_\Lambda(t_n-t_0)\bigr)=c'e^{\Lambda(t_n-t_0)}\le2\Gamma e^{2\Gamma L\bar\beta(T-t_0)}(E_{\mathrm s}+D)

を得る。▨

定理 2.6 (Dahlquist の等価定理). 整合的な線形多段法がゼロ安定であることは、収束するための必要十分条件である。

証明.kk段の線形多段法(α0,…,αk;β0,…,βk)(\alpha_0,\dots,\alpha_k;\beta_0,\dots,\beta_k)は整合的であるとする。

必要性を示す。線形多段法が収束し、ゼロ安定でないとする。補題 2.2をaj:=αja_j:=\alpha_jとして適用すると、∑jαjun+j=0\sum_j\alpha_ju_{n+j}=0(n∈N≥0n\in\N)を満たす有界でない複素数列(un)(u_n)が存在する。C\CをR2\R^2と同一視し、ノルムとして絶対値を用いる。d=2d=2、I=J=RI=J=\R、f=0f=0、L=0L=0、t0=0t_0=0、T=1T=1、y=0y=0とし、定義 2.4のN1N_1を取る。f=0f=0では近似値の方程式はv=−∑j<kαjyn+jv=-\sum_{j<k}\alpha_jy_{n+j}であり、ただ一つの解をもつ。mN:=max⁡0≤n≤N∣un∣m_N:=\max_{0\le n\le N}|u_n|と置き、mN>0m_N>0ならばcN:=mN−1/2c_N:=m_N^{-1/2}、mN=0m_N=0ならばcN:=0c_N:=0とし、開始値をyj(N):=cNujy^{(N)}_j:=c_Nu_j(0≤j<k0\le j<k)とする。αj\alpha_jは実数であるから(cNun)(c_Nu_n)は近似値の漸化式を満たし、0≤n≤N0\le n\le Nでyn(N)=cNuny^{(N)}_n=c_Nu_nである。(un)(u_n)は有界でないのでmN→∞m_N\to\inftyであり、mN>0m_N>0となるNNについて

εN=cNmax⁡j<k∣uj∣≤mN−1/2mk−1,EN=cNmN=mN1/2\varepsilon_N=c_N\max_{j<k}|u_j|\le m_N^{-1/2}m_{k-1},\qquad E_N=c_Nm_N=m_N^{1/2}

である。したがってεN→0\varepsilon_N\to0でありEN→∞E_N\to\inftyであるから、εN→0\varepsilon_N\to0ならばEN→0E_N\to0であるという収束の条件に反する。

十分性を示す。線形多段法がゼロ安定であるとする。定義 2.4のdd、ノルム、II、ff、LL、t0<Tt_0<T、JJ、yyを取り、Γ\Gamma、βˉ\bar\betaを定理 2.5のとおりとする。整数N1≥kN_1\ge kを2Γ(T−t0)L∣βk∣≤N12\Gamma(T-t_0)L|\beta_k|\le N_1を満たすように取ると、N≥N1N\ge N_1について2ΓhNL∣βk∣≤12\Gamma h_NL|\beta_k|\le1である。定理 2.5により近似値はすべて定まり、

EN≤2Γe2ΓLβˉ(T−t0)(εN+∑m=0N−k∥dhN(tm)∥)E_N\le2\Gamma e^{2\Gamma L\bar\beta(T-t_0)}\Bigl(\varepsilon_N+\sum_{m=0}^{N-k}\|d_{h_N}(t_m)\|\Bigr)

である。線形多段法は整合的でありyyはC1C^1級であるから、定義 1.2 (3)により

ηN:=sup⁡{∥dhN(t)∥hN ∣ t∈[t0,T−khN]}\eta_N:=\sup\Bigl\{\frac{\|d_{h_N}(t)\|}{h_N}\ \Big|\ t\in[t_0,T-kh_N]\Bigr\}

はN→∞N\to\inftyで00に収束する。0≤m≤N−k0\le m\le N-kについてtm∈[t0,T−khN]t_m\in[t_0,T-kh_N]であり、(N−k+1)hN≤T−t0(N-k+1)h_N\le T-t_0であるから、和は(T−t0)ηN(T-t_0)\eta_N以下である。したがってεN→0\varepsilon_N\to0ならばEN→0E_N\to0である。▨

系 2.7.定理 2.5の仮定の下で、線形多段法の次数がp∈N≥1p\in\NN以上であり、y∈Cp+1(J;Rd)y\in C^{p+1}(J;\R^d)であるとする。M:=max⁡τ∈[t0,T]∥y(p+1)(τ)∥M:=\max_{\tau\in[t_0,T]}\|y^{(p+1)}(\tau)\|とし、KpK_pを補題 1.3 (1)の定数とする。このとき

max⁡0≤n≤N∥en∥≤2Γe2ΓLβˉ(T−t0)(max⁡0≤j<k∥ej∥+(T−t0)KpMhp)\max_{0\le n\le N}\|e_n\|\le2\Gamma e^{2\Gamma L\bar\beta(T-t_0)}\Bigl(\max_{0\le j<k}\|e_j\|+(T-t_0)K_pMh^p\Bigr)

が成り立つ。特に、C′≥0C'\ge0について開始値がmax⁡j<k∥ej∥≤C′hp\max_{j<k}\|e_j\|\le C'h^pを満たすならば、max⁡n∥en∥≤2Γe2ΓLβˉ(T−t0)(C′+(T−t0)KpM)hp\max_n\|e_n\|\le2\Gamma e^{2\Gamma L\bar\beta(T-t_0)}(C'+(T-t_0)K_pM)h^pである。

証明.命題 1.4 (2)によりC0=⋯=Cp=0C_0=\dots=C_p=0であるから、補題 1.3 (1)をq=pq=pとして適用すると∥dm∥≤KpMhp+1\|d_m\|\le K_pMh^{p+1}である。(N−k+1)h≤T−t0(N-k+1)h\le T-t_0により∑m=0N−k∥dm∥≤(T−t0)KpMhp\sum_{m=0}^{N-k}\|d_m\|\le(T-t_0)K_pMh^pであり、定理 2.5に代入すると主張を得る。▨

例 2.8.k=2k=2、(α0,α1,α2)=(−5,4,1)(\alpha_0,\alpha_1,\alpha_2)=(-5,4,1)、(β0,β1,β2)=(2,4,0)(\beta_0,\beta_1,\beta_2)=(2,4,0)の陽的線形多段法

yn+2+4yn+1−5yn=h(4fn+1+2fn),ρ(ζ)=(ζ−1)(ζ+5),σ(ζ)=4ζ+2y_{n+2}+4y_{n+1}-5y_n=h(4f_{n+1}+2f_n),\qquad\rho(\zeta)=(\zeta-1)(\zeta+5),\quad\sigma(\zeta)=4\zeta+2

を考える。係数は

C0=0,C1=6−6=0,C2=4−4=0,C3=2−2=0,C4=2024−46=16C_0=0,\quad C_1=6-6=0,\quad C_2=4-4=0,\quad C_3=2-2=0,\quad C_4=\frac{20}{24}-\frac46=\frac16

であるから、命題 1.4 (2)によりこの方法の次数は33であり、命題 1.4 (1)により整合的である。ρ\rhoは根−5-5をもつのでゼロ安定でなく、定理 2.6により収束しない。

d=1d=1、f=0f=0、t0=0t_0=0、T=1T=1、h=1/Nh=1/N、解y=0y=0とし、開始値をy0=0y_0=0、y1=hy_1=hとする。近似値は漸化式yn+2=−4yn+1+5yny_{n+2}=-4y_{n+1}+5y_nで定まり、yn=h6(1−(−5)n)y_n=\frac h6\bigl(1-(-5)^n\bigr)である。開始値の誤差はhhであり、∣yN∣≥(5N−1)/(6N)|y_N|\ge(5^N-1)/(6N)はN→∞N\to\inftyで発散する。

d=1d=1、f(t,u)=uf(t,u)=u、解y(τ)=eτy(\tau)=e^\tau、t0=0t_0=0、T=1T=1、h=1/Nh=1/Nとし、開始値を厳密値y0=1y_0=1、y1=ehy_1=e^hとする。近似値は漸化式yn+2=−(4−4h)yn+1+(5+2h)yny_{n+2}=-(4-4h)y_{n+1}+(5+2h)y_nで定まる。値を有効数字7桁に丸めると次のとおりである。

NN yNy_N ∣e−yN∣\lvert e-y_N\rvert
55 2.7343272.734327 1.604520×10−21.604520\times10^{-2}
1010 −1.271974×10−1-1.271974\times10^{-1} 2.8454792.845479
2020 −1.622493×106-1.622493\times10^{6} 1.622496×1061.622496\times10^{6}
4040 −9.344203×1018-9.344203\times10^{18} 9.344203×10189.344203\times10^{18}

3 Adams 法と後退差分公式

定義 3.1.k∈N≥1k\in\NNとする。m∈N≥0m\in\Nについて、節点0,1,…,m0,1,\dots,mの Lagrange 基底をℓ0(m),…,ℓm(m)\ell^{(m)}_0,\dots,\ell^{(m)}_mとする。

  1. αk=1\alpha_k=1、αk−1=−1\alpha_{k-1}=-1、j<k−1j<k-1でαj=0\alpha_j=0とし、βj:=∫k−1kℓj(k−1)(s) ds\beta_j:=\int_{k-1}^k\ell^{(k-1)}_j(s)\,ds(0≤j<k0\le j<k)、βk:=0\beta_k:=0としたkk段の線形多段法を kk段 Adams–Bashforth 法 (kk-step Adams–Bashforth method) という。
  2. αk=1\alpha_k=1、αk−1=−1\alpha_{k-1}=-1、j<k−1j<k-1でαj=0\alpha_j=0とし、βj:=∫k−1kℓj(k)(s) ds\beta_j:=\int_{k-1}^k\ell^{(k)}_j(s)\,ds(0≤j≤k0\le j\le k)としたkk段の線形多段法を kk段 Adams–Moulton 法 (kk-step Adams–Moulton method) という。
  3. ℓk(k)(s)=∏i<k(s−i)/k!\ell^{(k)}_k(s)=\prod_{i<k}(s-i)/k!のs=ks=kにおける微分係数はck:=∑r=1k1rc_k:=\sum_{r=1}^k\frac1rである。αj:=ℓj(k)′(k)/ck\alpha_j:=\ell^{(k)\prime}_j(k)/c_k(0≤j≤k0\le j\le k)、βk:=1/ck\beta_k:=1/c_k、j<kj<kでβj:=0\beta_j:=0としたkk段の線形多段法を kk段後退差分公式 (kk-step backward differentiation formula) という。

命題 3.2.k∈N≥1k\in\NNとする。

  1. kk段 Adams–Bashforth 法の次数はkk以上である。
  2. kk段 Adams–Moulton 法の次数はk+1k+1以上である。
  3. kk段後退差分公式の次数はkk以上である。
  4. kk段 Adams–Bashforth 法とkk段 Adams–Moulton 法はゼロ安定である。

証明.i∈N≥0i\in\Nとし、y(τ):=τiy(\tau):=\tau^iと置く。補題 1.3 (2)をs=0s=0、h=1h=1として適用すると∑jαjy(j)−∑jβjy′(j)=i! Ci\sum_j\alpha_jy(j)-\sum_j\beta_jy'(j)=i!\,C_iである。命題 1.4 (2)により、次数がpp以上であることを示すには、i≤pi\le pについてこの左辺が00であることを示せばよい。

(1)を示す。i≤ki\le kとするとy′y'は次数k−1k-1以下の多項式である。∑j<ky′(j)ℓj(k−1)\sum_{j<k}y'(j)\ell^{(k-1)}_jも次数k−1k-1以下であり、相異なる節点0,…,k−10,\dots,k-1でy′y'と同じ値をとるから、§E20.12 定理 1.3 (2)によりy′=∑j<ky′(j)ℓj(k−1)y'=\sum_{j<k}y'(j)\ell^{(k-1)}_jである。[k−1,k][k-1,k]で積分するとy(k)−y(k−1)=∑j<kβjy′(j)y(k)-y(k-1)=\sum_{j<k}\beta_jy'(j)である。αk=1\alpha_k=1、αk−1=−1\alpha_{k-1}=-1、j<k−1j<k-1でαj=0\alpha_j=0、βk=0\beta_k=0であるから、∑jαjy(j)−∑jβjy′(j)=0\sum_j\alpha_jy(j)-\sum_j\beta_jy'(j)=0である。

(2)を示す。i≤k+1i\le k+1とするとy′y'は次数kk以下であり、節点0,…,k0,\dots,kについて同じ議論によりy(k)−y(k−1)=∑j≤kβjy′(j)y(k)-y(k-1)=\sum_{j\le k}\beta_jy'(j)である。

(3)を示す。i≤ki\le kとすると、§E20.12 定理 1.3 (2)によりy=∑j≤ky(j)ℓj(k)y=\sum_{j\le k}y(j)\ell^{(k)}_jであり、s=ks=kで微分するとy′(k)=∑jℓj(k)′(k)y(j)=ck∑jαjy(j)y'(k)=\sum_j\ell^{(k)\prime}_j(k)y(j)=c_k\sum_j\alpha_jy(j)である。したがって∑jαjy(j)−βky′(k)=0\sum_j\alpha_jy(j)-\beta_ky'(k)=0である。

(4)を示す。Adams 法の第一特性多項式はζk−1(ζ−1)\zeta^{k-1}(\zeta-1)であり、11はその単根である。k=1k=1ならば根は11だけであり、k≥2k\ge2ならば11以外の根は00だけである。したがって第一特性多項式は根条件を満たす。▨

命題 3.3.

  1. 2 段 Adams–Bashforth 法は yn+2−yn+1=h2(3fn+1−fn),ρ(ζ)=ζ2−ζ,σ(ζ)=3ζ−12y_{n+2}-y_{n+1}=\frac h2(3f_{n+1}-f_n),\qquad\rho(\zeta)=\zeta^2-\zeta,\quad\sigma(\zeta)=\frac{3\zeta-1}2 であり、C3=512C_3=\frac5{12}、σ(1)=1\sigma(1)=1である。
  2. 1 段 Adams–Moulton 法はρ(ζ)=ζ−1\rho(\zeta)=\zeta-1、σ(ζ)=ζ+12\sigma(\zeta)=\frac{\zeta+1}2であり、近似値の方程式は台形法の方程式(§E20.28 定義 2.1)に一致する。C3=−112C_3=-\frac1{12}、σ(1)=1\sigma(1)=1である。
  3. 2 段後退差分公式は yn+2−43yn+1+13yn=23hfn+2,ρ(ζ)=(ζ−1)(ζ−13),σ(ζ)=23ζ2y_{n+2}-\frac43y_{n+1}+\frac13y_n=\frac23hf_{n+2},\qquad\rho(\zeta)=(\zeta-1)\Bigl(\zeta-\frac13\Bigr),\quad\sigma(\zeta)=\frac23\zeta^2 であり、C3=−29C_3=-\frac29、σ(1)=23\sigma(1)=\frac23である。

三つの方法はいずれもゼロ安定であって次数は22であり、誤差定数はそれぞれ512\frac5{12}、−112-\frac1{12}、−13-\frac13である。

証明.(1)を示す。ℓ0(1)(s)=1−s\ell^{(1)}_0(s)=1-s、ℓ1(1)(s)=s\ell^{(1)}_1(s)=sであり、∫12(1−s) ds=−12\int_1^2(1-s)\,ds=-\frac12、∫12s ds=32\int_1^2s\,ds=\frac32である。C3=∑jαjj3/6−∑jβjj2/2=8−16−32⋅12=512C_3=\sum_j\alpha_jj^3/6-\sum_j\beta_jj^2/2=\frac{8-1}6-\frac32\cdot\frac12=\frac5{12}である。

(2)を示す。∫01(1−s) ds=∫01s ds=12\int_0^1(1-s)\,ds=\int_0^1s\,ds=\frac12であるから、近似値の方程式はyn+1=yn+h2(f(tn,yn)+f(tn+1,yn+1))y_{n+1}=y_n+\frac h2(f(t_n,y_n)+f(t_{n+1},y_{n+1}))である。C3=16−12⋅12=−112C_3=\frac16-\frac12\cdot\frac12=-\frac1{12}である。

(3)を示す。ℓ0(2)(s)=(s−1)(s−2)2\ell^{(2)}_0(s)=\frac{(s-1)(s-2)}2、ℓ1(2)(s)=−s(s−2)\ell^{(2)}_1(s)=-s(s-2)、ℓ2(2)(s)=s(s−1)2\ell^{(2)}_2(s)=\frac{s(s-1)}2のs=2s=2における微分係数は12\frac12、−2-2、32\frac32であり、c2=32c_2=\frac32であるから、(α0,α1,α2)=(13,−43,1)(\alpha_0,\alpha_1,\alpha_2)=(\frac13,-\frac43,1)、β2=23\beta_2=\frac23である。C3=8−4/36−23⋅42=109−43=−29C_3=\frac{8-4/3}6-\frac23\cdot\frac42=\frac{10}9-\frac43=-\frac29である。

命題 3.2により三つの方法の次数は22以上であり、C3≠0C_3\ne0であるから命題 1.4 (2)により次数は22である。Adams 法は命題 3.2 (4)によりゼロ安定であり、2 段後退差分公式の第一特性多項式の根11、13\frac13は根条件を満たす。誤差定数はC3/σ(1)C_3/\sigma(1)であり、2 段後退差分公式では−29/23=−13-\frac29\big/\frac23=-\frac13である。▨

4 絶対安定性

定義 4.1.kk段の線形多段法を取り、ρ,σ\rho,\sigmaをその特性多項式とする。z∈Cz\in\Cに対してπz(ζ):=ρ(ζ)−zσ(ζ)\pi_z(\zeta):=\rho(\zeta)-z\sigma(\zeta)を線形多段法の 安定多項式 (stability polynomial) という。πz\pi_zのζk\zeta^kの係数は1−zβk1-z\beta_kである。

S:={z∈C∣1−zβk≠0 かつ πz は根条件を満たす}S:=\{z\in\C\mid1-z\beta_k\ne0\ \text{かつ}\ \pi_z\ \text{は根条件を満たす}\}

を線形多段法の 絶対安定領域 (region of absolute stability) という。{z∈C∣Re⁡z≤0}⊆S\{z\in\C\mid\operatorname{Re}z\le0\}\subseteq Sであるとき、線形多段法は A 安定 (A-stable) であるという。

命題 4.2.kk段の線形多段法を取り、πz\pi_zをその安定多項式、SSを絶対安定領域とする。λ∈C\lambda\in\C、h>0h>0、z:=hλz:=h\lambdaとし、C\CをR2\R^2と同一視してf(t,u):=λuf(t,u):=\lambda u((t,u)∈R×C(t,u)\in\R\times\C)に線形多段法を刻みhhで適用する。

  1. 1−zβk≠01-z\beta_k\ne0ならば、任意の開始値y0,…,yk−1∈Cy_0,\dots,y_{k-1}\in\Cから近似値yny_nがすべてのn∈N≥0n\in\Nについて定まり、∑j=0k(αj−zβj)yn+j=0\sum_{j=0}^k(\alpha_j-z\beta_j)y_{n+j}=0を満たす。1−zβk=01-z\beta_k=0ならば、任意のyn,…,yn+k−1∈Cy_n,\dots,y_{n+k-1}\in\Cに対して近似値の方程式はただ一つの解をもつことがない。
  2. z∈Sz\in Sであることと、1−zβk≠01-z\beta_k\ne0であって任意の開始値から定まる近似値の列が有界であることは同値である。
  3. 1−zβk≠01-z\beta_k\ne0とする。πz\pi_zの根がすべて∣ζ∣<1|\zeta|<1を満たすことと、任意の開始値から定まる近似値の列が00に収束することは同値である。

証明. 近似値の方程式は(1−zβk)v=−∑j<k(αj−zβj)yn+j(1-z\beta_k)v=-\sum_{j<k}(\alpha_j-z\beta_j)y_{n+j}である。1−zβk≠01-z\beta_k\ne0ならばただ一つの解をもち、(1)の漸化式を得る。1−zβk=01-z\beta_k=0ならば左辺は00であり、右辺が00でなければ解をもたず、右辺が00ならば任意のvvが解である。1−zβk≠01-z\beta_k\ne0のとき、aj:=(αj−zβj)/(1−zβk)a_j:=(\alpha_j-z\beta_j)/(1-z\beta_k)として補題 2.2を適用する。∑jajζj=πz(ζ)/(1−zβk)\sum_ja_j\zeta^j=\pi_z(\zeta)/(1-z\beta_k)の根は重複度を込めてπz\pi_zの根に等しいから、(2)は二条件の同値から、(3)は補題 2.2 (2)から従う。▨

命題 4.3. 各方法の絶対安定領域をSSとする。

  1. 2 段 Adams–Bashforth 法について、z∈Cz\in\Cに対して安定多項式の根は ζ±(z)=1+32z±1+z+94z22\zeta_\pm(z)=\frac{1+\frac32z\pm\sqrt{1+z+\frac94z^2}}2 である。平方根の枝を替えるとζ+\zeta_+とζ−\zeta_-が入れ替わる。SSは∣ζ+(z)∣≤1|\zeta_+(z)|\le1かつ∣ζ−(z)∣≤1|\zeta_-(z)|\le1であり、1+z+94z2=01+z+\frac94z^2=0ならば∣ζ+(z)∣<1|\zeta_+(z)|<1であるzzの全体である。S∩R=[−1,0]S\cap\R=[-1,0]であり、2 段 Adams–Bashforth 法は A 安定でない。
  2. 1 段 Adams–Moulton 法について、z≠2z\ne2ならば安定多項式の根は台形法の安定関数の値R1/2(z)R_{1/2}(z)である。S={z∈C∣Re⁡z≤0}S=\{z\in\C\mid\operatorname{Re}z\le0\}であり、1 段 Adams–Moulton 法は A 安定である。
  3. 2 段後退差分公式について、z≠32z\ne\frac32ならば安定多項式の根は ζ±(z)=2±1+2z3−2z\zeta_\pm(z)=\frac{2\pm\sqrt{1+2z}}{3-2z} である。SSはz≠32z\ne\frac32、∣ζ+(z)∣≤1|\zeta_+(z)|\le1かつ∣ζ−(z)∣≤1|\zeta_-(z)|\le1であり、1+2z=01+2z=0ならば∣ζ+(z)∣<1|\zeta_+(z)|<1であるzzの全体である。2 段後退差分公式は A 安定である。

証明.(1)を示す。安定多項式はζ2−(1+32z)ζ+z2\zeta^2-(1+\frac32z)\zeta+\frac z2であり、最高次の係数は11である。判別式は(1+32z)2−2z=1+z+94z2(1+\frac32z)^2-2z=1+z+\frac94z^2であるから、根はζ±(z)\zeta_\pm(z)であり、判別式が00であることと重根をもつことは同値である。したがって根条件は述べた条件と同値である。z=x∈Rz=x\in\Rとする。1+x+94x21+x+\frac94x^2はxxの二次式として判別式1−9<01-9<0をもつので正であり、根r−<r+r_-<r_+は相異なる実数で、r−+r+=1+32xr_-+r_+=1+\frac32x、(1−r−)(1−r+)=−x(1-r_-)(1-r_+)=-x、(1+r−)(1+r+)=2+2x(1+r_-)(1+r_+)=2+2xである。x>0x>0ならば(1−r−)(1−r+)<0(1-r_-)(1-r_+)<0によりr−<1<r+r_-<1<r_+であり、x∉Sx\notin Sである。x<−1x<-1ならば(1+r−)(1+r+)<0(1+r_-)(1+r_+)<0によりr−<−1<r+r_-<-1<r_+であり、x∉Sx\notin Sである。−1≤x≤0-1\le x\le0とするとr−+r+∈[−12,1]r_-+r_+\in[-\frac12,1]である。(1−r−)(1−r+)≥0(1-r_-)(1-r_+)\ge0であるからr+≤1r_+\le1またはr−≥1r_-\ge1であり、r−≥1r_-\ge1ならばr−+r+>2r_-+r_+>2となるのでr+≤1r_+\le1である。(1+r−)(1+r+)≥0(1+r_-)(1+r_+)\ge0であるからr−≥−1r_-\ge-1またはr+≤−1r_+\le-1であり、r+≤−1r_+\le-1ならばr−+r+<−2r_-+r_+<-2となるのでr−≥−1r_-\ge-1である。二根は相異なるから根条件が成り立ち、x∈Sx\in Sである。したがってS∩R=[−1,0]S\cap\R=[-1,0]であり、Re⁡(−2)≤0\operatorname{Re}(-2)\le0かつ−2∉S-2\notin Sであるから A 安定でない。

(2)を示す。安定多項式は(1−z2)ζ−(1+z2)(1-\frac z2)\zeta-(1+\frac z2)であり、最高次の係数が00であることとz=2z=2は同値である。z≠2z\ne2ならば根は1+z/21−z/2=R1/2(z)\frac{1+z/2}{1-z/2}=R_{1/2}(z)であり(§E20.28 命題 5.5)、一次式の根は単根であるから、根条件は∣R1/2(z)∣≤1|R_{1/2}(z)|\le1と同値である。§E20.28 命題 5.5 (1)をθ=12\theta=\frac12として適用すると、これはRe⁡z≤0\operatorname{Re}z\le0と同値である。Re⁡2>0\operatorname{Re}2>0であるからS={z∣Re⁡z≤0}S=\{z\mid\operatorname{Re}z\le0\}である。

(3)を示す。安定多項式の33倍は(3−2z)ζ2−4ζ+1(3-2z)\zeta^2-4\zeta+1であり、最高次の係数が00であることとz=32z=\frac32は同値である。z≠32z\ne\frac32ならば判別式の14\frac14は4−(3−2z)=1+2z4-(3-2z)=1+2zであり、根はζ±(z)\zeta_\pm(z)である。SSの記述は(1)と同じ理由による。Re⁡z≤0\operatorname{Re}z\le0とするとz≠32z\ne\frac32である。ζ\zetaを根とすると、定数項が11であるからζ≠0\zeta\ne0であり、w:=1/ζw:=1/\zetaと置いて方程式をζ2\zeta^2で割ると3−2z−4w+w2=03-2z-4w+w^2=0、すなわちz=3−4w+w22z=\frac{3-4w+w^2}2である。u:=Re⁡wu:=\operatorname{Re}wとするとRe⁡(w2)=2u2−∣w∣2\operatorname{Re}(w^2)=2u^2-|w|^2であるから

Re⁡z=(1−u)2+1−∣w∣22\operatorname{Re}z=(1-u)^2+\frac{1-|w|^2}2

である。∣ζ∣>1|\zeta|>1ならば∣w∣<1|w|<1でありRe⁡z>0\operatorname{Re}z>0となるので、∣ζ∣≤1|\zeta|\le1である。∣ζ∣=1|\zeta|=1ならばRe⁡z=(1−u)2≤0\operatorname{Re}z=(1-u)^2\le0によりu=1u=1であり、∣w∣=1|w|=1と合わせてw=1w=1、ζ=1\zeta=1、z=0z=0である。z=0z=0の安定多項式はρ\rhoであり、根11は単根である。したがってRe⁡z≤0\operatorname{Re}z\le0ならばz∈Sz\in Sである。▨

注意 4.4. 2 段 Adams–Bashforth 法でn≥1n\ge1についてyn+2y_{n+2}を求めるとき、fnf_nは前の歩で評価済みであり、新たに評価するのはfn+1f_{n+1}だけである。この点で前進 Euler 法と同じである。前進 Euler 法の絶対安定領域と実軸の共通部分は[−2,0][-2,0]であり(§E20.28 系 5.6)、次数22の 2 段 Adams–Bashforth 法では[−1,0][-1,0]である。1 段 Adams–Moulton 法の絶対安定領域は、一段法としての台形法の絶対安定領域{z∈C∖{2}∣∣R1/2(z)∣≤1}\{z\in\C\setminus\{2\}\mid|R_{1/2}(z)|\le1\}に等しい。台形法と 2 段後退差分公式はともに A 安定であるが、実数x→−∞x\to-\inftyで台形法の根R1/2(x)R_{1/2}(x)は−1-1に収束し(§E20.28 命題 5.5 (3))、2 段後退差分公式の根の絶対値は00に収束する(問題 6.3)。

5 A 安定性と次数の上限

補題 5.1.P:={w∈C∣Re⁡w>0}P:=\{w\in\C\mid\operatorname{Re}w>0\}とする。a≥0a\ge0とし、GGをPP上の正則関数で、任意のw∈Pw\in PについてRe⁡G(w)≥0\operatorname{Re}G(w)\ge0を満たすものとする。K≥0K\ge0とR0>0R_0>0が存在して、∣w∣≥R0|w|\ge R_0を満たす任意のw∈Pw\in Pについて∣G(w)−aw∣≤K/∣w∣|G(w)-aw|\le K/|w|が成り立つとする。このとき任意のw∈Pw\in PについてRe⁡(G(w)−aw)≥0\operatorname{Re}(G(w)-aw)\ge0である。

証明.H(w):=G(w)−awH(w):=G(w)-awと置く。w0∈Pw_0\in Pを取り、0<ε<Re⁡w00<\varepsilon<\operatorname{Re}w_0、R>max⁡{R0,∣w0∣}R>\max\{R_0,|w_0|\}とする。D:={w∣Re⁡w>ε, ∣w∣<R}D:=\{w\mid\operatorname{Re}w>\varepsilon,\ |w|<R\}は凸な開集合でありw0w_0を含むので領域である。Dˉ:={w∣Re⁡w≥ε, ∣w∣≤R}\bar D:=\{w\mid\operatorname{Re}w\ge\varepsilon,\ |w|\le R\}はC=R2\C=\R^2の有界閉集合であるから§E2.9 定理 4.3によりコンパクトであり、D⊆Dˉ⊆PD\subseteq\bar D\subseteq Pである。F:=exp⁡∘(−H)F:=\exp\circ(-H)は§E5.3 命題 1.8によりPP上で正則であり、§E5.3 命題 1.7 (1)により∣F(w)∣=e−Re⁡H(w)|F(w)|=e^{-\operatorname{Re}H(w)}である。§E5.9 命題 1.1により∣F∣|F|はDˉ\bar D上で点w∗w^*において最大値をとる。

w∈Dˉ∖Dw\in\bar D\setminus DはRe⁡w=ε\operatorname{Re}w=\varepsilonまたは∣w∣=R|w|=Rを満たす。前者ではRe⁡H(w)=Re⁡G(w)−aε≥−aε\operatorname{Re}H(w)=\operatorname{Re}G(w)-a\varepsilon\ge-a\varepsilonであり、後者では∣w∣≥R0|w|\ge R_0により∣H(w)∣≤K/R|H(w)|\le K/RであるからRe⁡H(w)≥−K/R\operatorname{Re}H(w)\ge-K/Rである。μ:=max⁡{aε,K/R}\mu:=\max\{a\varepsilon,K/R\}と置くと、Dˉ∖D\bar D\setminus D上で∣F∣≤eμ|F|\le e^\muである。

w∗∉Dw^*\notin Dならば∣F(w0)∣≤∣F(w∗)∣≤eμ|F(w_0)|\le|F(w^*)|\le e^\muである。w∗∈Dw^*\in Dならば、DDは開集合であるからD(w∗,r)⊆DD(w^*,r)\subseteq Dを満たすr>0r>0が存在し、その上で∣F∣≤∣F(w∗)∣|F|\le|F(w^*)|であるので、§E5.9 定理 3.1によりFFはDD上で定数である。w1:=ε+iIm⁡w0w_1:=\varepsilon+i\operatorname{Im}w_0は∣w1∣≤∣w0∣<R|w_1|\le|w_0|<RによりDˉ∖D\bar D\setminus Dに属し、0<s<Re⁡w0−ε0<s<\operatorname{Re}w_0-\varepsilonについて∣w1+s∣≤∣w0∣|w_1+s|\le|w_0|であるからw1+s∈Dw_1+s\in Dである。したがって∣F(w1+s)∣=∣F(w0)∣|F(w_1+s)|=|F(w_0)|であり、s→+0s\to+0として∣F(w0)∣=∣F(w1)∣≤eμ|F(w_0)|=|F(w_1)|\le e^\muである。

いずれの場合もe−Re⁡H(w0)≤eμe^{-\operatorname{Re}H(w_0)}\le e^\mu、すなわちRe⁡H(w0)≥−max⁡{aε,K/R}\operatorname{Re}H(w_0)\ge-\max\{a\varepsilon,K/R\}である。R→∞R\to\inftyとし、次にε→+0\varepsilon\to+0とするとRe⁡H(w0)≥0\operatorname{Re}H(w_0)\ge0を得る。▨

定理 5.2 (Dahlquist の第二障壁). A 安定な整合的線形多段法について、次が成り立つ。

  1. 次数は33以上でない。
  2. 次数が22以上ならばσ(1)≠0\sigma(1)\ne0かつC3/σ(1)≤−112C_3/\sigma(1)\le-\frac1{12}である。特に、次数が22ならば誤差定数の絶対値は112\frac1{12}以上である。

1 段 Adams–Moulton 法は A 安定で次数が22の線形多段法であり、その誤差定数は−112-\frac1{12}である。

証明. 線形多段法をkk段とし、ρ,σ\rho,\sigmaをその特性多項式、C0,C1,…C_0,C_1,\dotsを係数とする。z=0z=0はRe⁡z≤0\operatorname{Re}z\le0を満たすから0∈S0\in Sであり、π0=ρ\pi_0=\rhoは根条件を満たす。命題 1.4 (1)によりρ(1)=0\rho(1)=0かつρ′(1)=σ(1)\rho'(1)=\sigma(1)である。11は絶対値11の根であるから単根であり、σ(1)=ρ′(1)≠0\sigma(1)=\rho'(1)\ne0である。

(2)を示す。次数が22以上であるとし、E:=C3/σ(1)E:=C_3/\sigma(1)、P:={w∈C∣Re⁡w>0}P:=\{w\in\C\mid\operatorname{Re}w>0\}と置く。実係数の多項式

ρ~(w):=∑j=0kαj(w+1)j(w−1)k−j,σ~(w):=∑j=0kβj(w+1)j(w−1)k−j\tilde\rho(w):=\sum_{j=0}^k\alpha_j(w+1)^j(w-1)^{k-j},\qquad\tilde\sigma(w):=\sum_{j=0}^k\beta_j(w+1)^j(w-1)^{k-j}

を考えると、ρ~(1)=2k\tilde\rho(1)=2^kである。w∈P∖{1}w\in P\setminus\{1\}についてζ:=w+1w−1\zeta:=\frac{w+1}{w-1}と置くと、∣w+1∣2−∣w−1∣2=4Re⁡w>0|w+1|^2-|w-1|^2=4\operatorname{Re}w>0により∣ζ∣>1|\zeta|>1であり、ρ~(w)=(w−1)kρ(ζ)\tilde\rho(w)=(w-1)^k\rho(\zeta)、σ~(w)=(w−1)kσ(ζ)\tilde\sigma(w)=(w-1)^k\sigma(\zeta)である。ρ\rhoは根条件を満たすのでρ(ζ)≠0\rho(\zeta)\ne0であり、ρ~\tilde\rhoはPP上で00にならない。したがってG:=σ~/ρ~G:=\tilde\sigma/\tilde\rhoはPP上で正則であり、w∈P∖{1}w\in P\setminus\{1\}についてG(w)=σ(ζ)/ρ(ζ)G(w)=\sigma(\zeta)/\rho(\zeta)である。

w∈P∖{1}w\in P\setminus\{1\}がRe⁡G(w)<0\operatorname{Re}G(w)<0を満たすとすると、z:=1/G(w)z:=1/G(w)はRe⁡z=Re⁡G(w)‾/∣G(w)∣2<0\operatorname{Re}z=\operatorname{Re}\overline{G(w)}/|G(w)|^2<0を満たし、πz(ζ)=ρ(ζ)(1−zG(w))=0\pi_z(\zeta)=\rho(\zeta)(1-zG(w))=0である。A 安定性によりz∈Sz\in Sであるからπz\pi_zの根は∣ζ∣≤1|\zeta|\le1を満たし、∣ζ∣>1|\zeta|>1と両立しない。したがってP∖{1}P\setminus\{1\}上でRe⁡G≥0\operatorname{Re}G\ge0であり、GGの連続性によりPP上でRe⁡G≥0\operatorname{Re}G\ge0である。

ρ(1+x)=∑q≥1rqxq\rho(1+x)=\sum_{q\ge1}r_qx^q、σ(1+x)=∑q≥0sqxq\sigma(1+x)=\sum_{q\ge0}s_qx^qと展開するとrq=∑jαj(jq)r_q=\sum_j\alpha_j\binom jq、sq=∑jβj(jq)s_q=\sum_j\beta_j\binom jqである。j22=(j2)+j2\frac{j^2}2=\binom j2+\frac j2、j36=(j3)+(j2)+j6\frac{j^3}6=\binom j3+\binom j2+\frac j6により

C1=r1−s0,C2=r2+r12−s1,C3=r3+r2+r16−s2−s12C_1=r_1-s_0,\qquad C_2=r_2+\frac{r_1}2-s_1,\qquad C_3=r_3+r_2+\frac{r_1}6-s_2-\frac{s_1}2

である。命題 1.4 (2)によりC1=C2=0C_1=C_2=0であるから、s0=r1=σ(1)s_0=r_1=\sigma(1)、s1=r2+r12s_1=r_2+\frac{r_1}2、C3=r3+r22−r112−s2C_3=r_3+\frac{r_2}2-\frac{r_1}{12}-s_2である。多項式Q(x):=ρ(1+x)/x=r1+r2x+r3x2+⋯Q(x):=\rho(1+x)/x=r_1+r_2x+r_3x^2+\cdotsはQ(0)=r1≠0Q(0)=r_1\ne0を満たす。c2:=−(112+E)c_2:=-(\frac1{12}+E)と置くと、σ(1+x)−(1+x2+c2x2)Q(x)\sigma(1+x)-(1+\frac x2+c_2x^2)Q(x)のx0,x1,x2x^0,x^1,x^2の係数は

s0−r1=0,s1−r2−r12=0,s2−r3−r22−c2r1=−C3+Er1=0s_0-r_1=0,\qquad s_1-r_2-\frac{r_1}2=0,\qquad s_2-r_3-\frac{r_2}2-c_2r_1=-C_3+Er_1=0

であるから、多項式VVが存在してσ(1+x)−(1+x2+c2x2)Q(x)=x3V(x)\sigma(1+x)-(1+\frac x2+c_2x^2)Q(x)=x^3V(x)である。

V(x)=∑ivixiV(x)=\sum_iv_ix^iと書き、MV:=∑i∣vi∣M_V:=\sum_i|v_i|と置く。δ∈(0,1]\delta\in(0,1]を、∣x∣≤δ|x|\le\deltaで∣Q(x)∣≥∣r1∣/2|Q(x)|\ge|r_1|/2となるように取る。0<∣x∣≤δ0<|x|\le\deltaならば∣V(x)∣≤MV|V(x)|\le M_V、ρ(1+x)=xQ(x)≠0\rho(1+x)=xQ(x)\ne0であり、

σ(1+x)ρ(1+x)=1x+12+c2x+x2V(x)Q(x),∣x2V(x)Q(x)∣≤2MV∣r1∣∣x∣2\frac{\sigma(1+x)}{\rho(1+x)}=\frac1x+\frac12+c_2x+\frac{x^2V(x)}{Q(x)},\qquad\Bigl|\frac{x^2V(x)}{Q(x)}\Bigr|\le\frac{2M_V}{|r_1|}|x|^2

である。R0:=max⁡{2,4/δ}R_0:=\max\{2,4/\delta\}とし、w∈Pw\in P、∣w∣≥R0|w|\ge R_0とする。x:=2w−1x:=\frac2{w-1}と置くと1+x=ζ1+x=\zetaであり、∣w−1∣≥∣w∣/2|w-1|\ge|w|/2により0<∣x∣≤4/∣w∣≤δ0<|x|\le4/|w|\le\deltaである。1x+12=w2\frac1x+\frac12=\frac w2、c2x=2c2w+2c2w(w−1)c_2x=\frac{2c_2}w+\frac{2c_2}{w(w-1)}であるから

G(w)−w2=Bw+ϵ(w),B:=−(16+2E),∣ϵ(w)∣≤K2∣w∣2,K2:=4∣c2∣+32MV∣r1∣G(w)-\frac w2=\frac Bw+\epsilon(w),\qquad B:=-\Bigl(\frac16+2E\Bigr),\qquad|\epsilon(w)|\le\frac{K_2}{|w|^2},\quad K_2:=4|c_2|+\frac{32M_V}{|r_1|}

である。

w∈Pw\in P、∣w∣≥R0|w|\ge R_0ならば∣G(w)−w2∣≤(∣B∣+K2/R0)/∣w∣|G(w)-\frac w2|\le(|B|+K_2/R_0)/|w|であるから、補題 5.1をa=12a=\frac12として適用すると、PP上でRe⁡(G(w)−w2)≥0\operatorname{Re}(G(w)-\frac w2)\ge0である。

ρ~,σ~\tilde\rho,\tilde\sigmaは実係数であるから、実数ξ≥R0\xi\ge R_0についてG(ξ)−ξ2G(\xi)-\frac\xi2は実数であり00以上である。したがってB+ξϵ(ξ)=ξ(G(ξ)−ξ2)≥0B+\xi\epsilon(\xi)=\xi(G(\xi)-\frac\xi2)\ge0であり、∣ξϵ(ξ)∣≤K2/ξ|\xi\epsilon(\xi)|\le K_2/\xiからξ→∞\xi\to\inftyとしてB≥0B\ge0を得る。B≥0B\ge0はE≤−112E\le-\frac1{12}と同値である。次数が22ならば誤差定数はEEであり、E<0E<0により∣E∣≥112|E|\ge\frac1{12}である。

(1)を示す。次数が33以上であるとすると、(2)によりC3/σ(1)≤−112C_3/\sigma(1)\le-\frac1{12}であり、命題 1.4 (2)によりC3=0C_3=0である。これらは両立しない。

1 段 Adams–Moulton 法については、命題 4.3 (2)により A 安定であり、命題 3.3により次数は22、誤差定数は−112-\frac1{12}である。▨

6 演習

問題 6.1.d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。I⊆RI\subseteq\Rを開区間、f ⁣:I×Rd→Rdf\colon I\times\R^d\to\R^dを連続写像とし、L≥0L\ge0が、任意のs∈Is\in Iとu,v∈Rdu,v\in\R^dについて∥f(s,u)−f(s,v)∥≤L∥u−v∥\|f(s,u)-f(s,v)\|\le L\|u-v\|を満たすとする。t0<Tt_0<Tとし、JJを[t0,T]⊆J⊆I[t_0,T]\subseteq J\subseteq Iを満たす開区間、y∈C3(J;Rd)y\in C^3(J;\R^d)を、任意のτ∈J\tau\in Jについてy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たす写像とする。整数N≥2N\ge2についてh:=(T−t0)/Nh:=(T-t_0)/N、tn:=t0+nht_n:=t_0+nhと置き、開始値をy0:=y(t0)y_0:=y(t_0)、y1:=y0+hf(t0,y0)y_1:=y_0+hf(t_0,y_0)として 2 段 Adams–Bashforth 法の近似値y0,…,yNy_0,\dots,y_Nを定める。C≥0C\ge0が存在して、任意の整数N≥2N\ge2についてmax⁡0≤n≤N∥yn−y(tn)∥≤Ch2\max_{0\le n\le N}\|y_n-y(t_n)\|\le Ch^2が成り立つことを示せ。

解答.

2 段 Adams–Bashforth 法は命題 3.3によりゼロ安定であり、次数は22である。Γ\Gammaをa0:=0a_0:=0、a1:=−1a_1:=-1とした補題 2.2 (1)の定数とする。β2=0\beta_2=0であるから、任意の整数N≥2N\ge2について2ΓhL∣β2∣≤12\Gamma hL|\beta_2|\le1であり、定理 2.5の仮定が成り立つ。したがって近似値y0,…,yNy_0,\dots,y_Nはすべて定まる。en:=yn−y(tn)e_n:=y_n-y(t_n)と置く。

e0=0e_0=0である。y0=y(t0)y_0=y(t_0)とy′(t0)=f(t0,y(t0))y'(t_0)=f(t_0,y(t_0))により

e1=−(y(t0+h)−y(t0)−hf(t0,y(t0)))e_1=-\bigl(y(t_0+h)-y(t_0)-hf(t_0,y(t_0))\bigr)

であり、右辺の括弧内は前進 Euler 法の、解yyの点t0t_0、刻みhhにおける局所打切り誤差δy(t0,h)\delta_y(t_0,h)である。I×RdI\times\R^dは開集合であり、y∈C2(J;Rd)y\in C^2(J;\R^d)、0<h≤T−t00<h\le T-t_0であるから、§E20.27 命題 1.6により、M2:=max⁡τ∈[t0,T]∥y′′(τ)∥M_2:=\max_{\tau\in[t_0,T]}\|y''(\tau)\|について∥e1∥≤M2h2/2\|e_1\|\le M_2h^2/2である。したがってmax⁡0≤j<2∥ej∥≤M2h2/2\max_{0\le j<2}\|e_j\|\le M_2h^2/2である。

y∈C3(J;Rd)y\in C^3(J;\R^d)であるから、系 2.7をp=2p=2、C′=M2/2C'=M_2/2として適用すると、M3:=max⁡τ∈[t0,T]∥y′′′(τ)∥M_3:=\max_{\tau\in[t_0,T]}\|y'''(\tau)\|について

max⁡0≤n≤N∥en∥≤2Γe4ΓL(T−t0)(M22+(T−t0)K2M3)h2\max_{0\le n\le N}\|e_n\|\le2\Gamma e^{4\Gamma L(T-t_0)}\Bigl(\frac{M_2}2+(T-t_0)K_2M_3\Bigr)h^2

である。ここでβˉ=32+12=2\bar\beta=\frac32+\frac12=2を用い、K2K_2は補題 1.3 (1)の定数である。Γ\GammaとK2K_2は 2 段 Adams–Bashforth 法の係数だけから定まりNNによらないので、右辺のh2h^2の係数をCCとすればよい。▨

問題 6.2. 2 段 Adams–Moulton 法の係数β0,β1,β2\beta_0,\beta_1,\beta_2を求めよ。この方法の次数が33であり、ゼロ安定であって、A 安定でないことを示せ。

解答.

ℓ0(2)(s)=(s−1)(s−2)2\ell^{(2)}_0(s)=\frac{(s-1)(s-2)}2、ℓ1(2)(s)=−s(s−2)\ell^{(2)}_1(s)=-s(s-2)、ℓ2(2)(s)=s(s−1)2\ell^{(2)}_2(s)=\frac{s(s-1)}2を[1,2][1,2]で積分すると

β0=12[s33−3s22+2s]12=−112,β1=[−s33+s2]12=23,β2=12[s33−s22]12=512\beta_0=\frac12\Bigl[\frac{s^3}3-\frac{3s^2}2+2s\Bigr]_1^2=-\frac1{12},\qquad\beta_1=\Bigl[-\frac{s^3}3+s^2\Bigr]_1^2=\frac23,\qquad\beta_2=\frac12\Bigl[\frac{s^3}3-\frac{s^2}2\Bigr]_1^2=\frac5{12}

であり、方法はyn+2−yn+1=h12(5fn+2+8fn+1−fn)y_{n+2}-y_{n+1}=\frac h{12}(5f_{n+2}+8f_{n+1}-f_n)である。命題 3.2 (2)により次数は33以上であり、

C4=16−124−16(23+8⋅512)=1524−1624=−124≠0C_4=\frac{16-1}{24}-\frac16\Bigl(\frac23+8\cdot\frac5{12}\Bigr)=\frac{15}{24}-\frac{16}{24}=-\frac1{24}\ne0

であるから、命題 1.4 (2)により次数は33である。命題 3.2 (4)によりゼロ安定である。次数が11以上であるから命題 1.4 (1)により整合的であり、A 安定ならば定理 5.2 (1)により次数は33以上でない。したがってこの方法は A 安定でない。▨

問題 6.3. 実数x<−12x<-\frac12について、2 段後退差分公式の安定多項式πx\pi_xの二つの根の絶対値がともに(3−2x)−1/2(3-2x)^{-1/2}であることを示せ。x=−100x=-100について、この値と 1 段 Adams–Moulton 法の安定多項式πx\pi_xの根を求めよ。

解答.

x<−12x<-\frac12ならばx≠32x\ne\frac32であり、命題 4.3 (3)により根はζ±=2±1+2x3−2x\zeta_\pm=\frac{2\pm\sqrt{1+2x}}{3-2x}である。1+2x<01+2x<0であるから1+2x=±i−1−2x\sqrt{1+2x}=\pm i\sqrt{-1-2x}であり、∣2±i−1−2x∣2=4−1−2x=3−2x|2\pm i\sqrt{-1-2x}|^2=4-1-2x=3-2xである。3−2x>03-2x>0により∣ζ±∣2=(3−2x)/(3−2x)2=1/(3−2x)|\zeta_\pm|^2=(3-2x)/(3-2x)^2=1/(3-2x)である。x=−100x=-100では∣ζ±∣=203−1/2=0.070186…|\zeta_\pm|=203^{-1/2}=0.070186\ldotsである。1 段 Adams–Moulton 法の根は、命題 4.3 (2)によりR1/2(−100)=1−501+50=−4951=−0.960784…R_{1/2}(-100)=\frac{1-50}{1+50}=-\frac{49}{51}=-0.960784\ldotsである。▨

前提記事