§E24.28一段法の安定性と収束

最終更新

初期値問題y′=−100yy'=-100y、y(0)=1y(0)=1の解e−100te^{-100t}は、ttとともに減衰する。この問題に前進 Euler 法を刻みh=0.021h=0.021で用いると、一歩ごとに近似値に1−100h=−1.11-100h=-1.1が掛かるので、第nn歩の近似値は(−1.1)n(-1.1)^nとなる。その絶対値は歩ごとに1.11.1倍になり、解とは逆に増大する。

一段法の近似値については、上の例から二つの問いが生じる。一つは、固定した有限区間の上で刻みを小さくしたとき、格子点での近似値が一様に解へ近づくかという問いである。もう一つは、刻みを固定して歩を重ねたとき、近似値がどのように振る舞うかという問いである。上の例は、この二つ目の問いに属する。本記事では、一段法の格子点での誤差の評価と、刻みを固定したときの近似値の振る舞いを扱う。

1 大域誤差と収束

定義 1.1.d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。Ω⊆R×Rd\Omega\subseteq\R\times\R^dを開集合、f ⁣:Ω→Rdf\colon\Omega\to\R^dを連続写像とし、一段法がffに対応させる増分関数の定義域をDD、一歩写像をΨh\Psi_hとする。t0<Tt_0<Tとし、J⊇[t0,T]J\supseteq[t_0,T]を開区間、y ⁣:J→Rdy\colon J\to\R^dを、任意のτ∈J\tau\in Jについて(τ,y(τ))∈Ω(\tau,y(\tau))\in\Omegaとy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たすC1C^1級写像とする。[t0,T][t_0,T]の格子t0<t1<⋯<tN=Tt_0<t_1<\dots<t_N=Tに対してhn:=tn+1−tnh_n:=t_{n+1}-t_n、H:=max⁡nhnH:=\max_nh_nと置く。

  1. y0,…,yN∈Rdy_0,\dots,y_N\in\R^dに対して、en:=yn−y(tn)e_n:=y_n-y(t_n)を格子点tnt_nにおける 大域誤差 (global error) といい、E:=max⁡0≤n≤N∥en∥E:=\max_{0\le n\le N}\|e_n\|を 大域誤差の大きさ (magnitude of the global error) という。
  2. 次の条件を満たすとき、一段法は解yyに対して 収束する (convergent) という。H1>0H_1>0が存在して、H≤H1H\le H_1を満たす[t0,T][t_0,T]の任意の格子について、初期値y0=y(t0)y_0=y(t_0)からyn+1:=Ψhn(tn,yn)y_{n+1}:=\Psi_{h_n}(t_n,y_n)で定まる近似値y0,…,yNy_0,\dots,y_Nがすべて定まり、任意のε>0\varepsilon>0に対してH2∈(0,H1]H_2\in(0,H_1]が存在して、H≤H2H\le H_2を満たす任意の格子についてE≤εE\le\varepsilonが成り立つ。
  3. p∈N≥1p\in\NNとする。C≥0C\ge0とH1>0H_1>0が存在して、H≤H1H\le H_1を満たす[t0,T][t_0,T]の任意の格子について、初期値y0=y(t0)y_0=y(t_0)から近似値y0,…,yNy_0,\dots,y_Nがすべて定まりE≤CHpE\le CH^pが成り立つとき、一段法は解yyに対して pp次で収束する (convergent of orderpp) という。

補題 1.2 (離散 Grönwall 不等式).Λ≥0\Lambda\ge0、N∈N≥1N\in\NNとし、実数t0<t1<⋯<tNt_0<t_1<\dots<t_Nに対してhn:=tn+1−tnh_n:=t_{n+1}-t_nと置く。E0≥0E_0\ge0とa0,…,aN−1≥0a_0,\dots,a_{N-1}\ge0に対して

Bn:=eΛ(tn−t0)E0+∑j=0n−1eΛ(tn−tj+1)aj(0≤n≤N)B_n:=e^{\Lambda(t_n-t_0)}E_0+\sum_{j=0}^{n-1}e^{\Lambda(t_n-t_{j+1})}a_j\qquad(0\le n\le N)

と置く。

  1. m∈{0,1,…,N}m\in\{0,1,\dots,N\}とし、E1,…,Em≥0E_1,\dots,E_m\ge0が0≤n<m0\le n<mについてEn+1≤(1+Λhn)En+anE_{n+1}\le(1+\Lambda h_n)E_n+a_nを満たすならば、0≤n≤m0\le n\le mについてEn≤BnE_n\le B_nが成り立つ。
  2. B0≤B1≤⋯≤BNB_0\le B_1\le\dots\le B_Nである。
  3. α≥0\alpha\ge0が0≤j<N0\le j<Nについてaj≤αhja_j\le\alpha h_jを満たすならば、0≤n≤N0\le n\le Nについて Bn≤eΛ(tn−t0)E0+α φΛ(tn−t0)B_n\le e^{\Lambda(t_n-t_0)}E_0+\alpha\,\varphi_\Lambda(t_n-t_0) が成り立つ。ここでτ≥0\tau\ge0に対して、Λ>0\Lambda>0ならばφΛ(τ):=(eΛτ−1)/Λ\varphi_\Lambda(\tau):=(e^{\Lambda\tau}-1)/\Lambda、Λ=0\Lambda=0ならばφ0(τ):=τ\varphi_0(\tau):=\tauである。

証明.0≤n<N0\le n<NについてBn+1=eΛhnBn+anB_{n+1}=e^{\Lambda h_n}B_n+a_nである。

主張 (1)を示す。B0=E0B_0=E_0である。n<mn<mについてEn≤BnE_n\le B_nであるとすると、1+Λhn≤eΛhn1+\Lambda h_n\le e^{\Lambda h_n}により

En+1≤(1+Λhn)Bn+an≤eΛhnBn+an=Bn+1E_{n+1}\le(1+\Lambda h_n)B_n+a_n\le e^{\Lambda h_n}B_n+a_n=B_{n+1}

である。nnについての帰納法により主張を得る。

主張 (2)はeΛhn≥1e^{\Lambda h_n}\ge1、Bn≥0B_n\ge0、an≥0a_n\ge0とBn+1=eΛhnBn+anB_{n+1}=e^{\Lambda h_n}B_n+a_nによる。

主張 (3)を示す。j<nj<nとs∈[tj,tj+1]s\in[t_j,t_{j+1}]についてeΛ(tn−tj+1)≤eΛ(tn−s)e^{\Lambda(t_n-t_{j+1})}\le e^{\Lambda(t_n-s)}であるから

∑j=0n−1eΛ(tn−tj+1)hj≤∑j=0n−1∫tjtj+1eΛ(tn−s) ds=∫t0tneΛ(tn−s) ds=φΛ(tn−t0)\sum_{j=0}^{n-1}e^{\Lambda(t_n-t_{j+1})}h_j\le\sum_{j=0}^{n-1}\int_{t_j}^{t_{j+1}}e^{\Lambda(t_n-s)}\,ds=\int_{t_0}^{t_n}e^{\Lambda(t_n-s)}\,ds=\varphi_\Lambda(t_n-t_0)

であり、aj≤αhja_j\le\alpha h_jを代入すると主張を得る。▨

定理 1.3.d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。Ω⊆R×Rd\Omega\subseteq\R\times\R^dを開集合、f ⁣:Ω→Rdf\colon\Omega\to\R^dを連続写像、Φ ⁣:D→Rd\Phi\colon D\to\R^dをffに対する増分関数、Ψh\Psi_hをその一歩写像とする。t0<Tt_0<Tとし、J⊇[t0,T]J\supseteq[t_0,T]を開区間、y ⁣:J→Rdy\colon J\to\R^dを、任意のτ∈J\tau\in Jについて(τ,y(τ))∈Ω(\tau,y(\tau))\in\Omegaとy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たすC1C^1級写像とする。[t0,T][t_0,T]の格子t0<⋯<tN=Tt_0<\dots<t_N=Tに対してhn:=tn+1−tnh_n:=t_{n+1}-t_nと置く。ρ∈(0,∞]\rho\in(0,\infty]とΛ≥0\Lambda\ge0を取り、各n∈{0,…,N−1}n\in\{0,\dots,N-1\}と、∥u−y(tn)∥≤ρ\|u-y(t_n)\|\le\rhoを満たす任意のu∈Rdu\in\R^d（ρ=∞\rho=\inftyのときは任意のu∈Rdu\in\R^d）について

(tn,u,hn)∈D,∥Ψhn(tn,u)−Ψhn(tn,y(tn))∥≤(1+Λhn)∥u−y(tn)∥(t_n,u,h_n)\in D,\qquad\|\Psi_{h_n}(t_n,u)-\Psi_{h_n}(t_n,y(t_n))\|\le(1+\Lambda h_n)\|u-y(t_n)\|

が成り立つとする。dn:=y(tn+1)−Ψhn(tn,y(tn))d_n:=y(t_{n+1})-\Psi_{h_n}(t_n,y(t_n))と置く。y0∈Rdy_0\in\R^dとr0,…,rN−1∈Rdr_0,\dots,r_{N-1}\in\R^dを取り、(tn,yn,hn)∈D(t_n,y_n,h_n)\in Dである限りyn+1:=Ψhn(tn,yn)+rny_{n+1}:=\Psi_{h_n}(t_n,y_n)+r_nと定め、En:=∥yn−y(tn)∥E_n:=\|y_n-y(t_n)\|と

Bn:=eΛ(tn−t0)E0+∑j=0n−1eΛ(tn−tj+1)(∥dj∥+∥rj∥)(0≤n≤N)B_n:=e^{\Lambda(t_n-t_0)}E_0+\sum_{j=0}^{n-1}e^{\Lambda(t_n-t_{j+1})}\bigl(\|d_j\|+\|r_j\|\bigr)\qquad(0\le n\le N)

と置く。BN≤ρB_N\le\rhoならば、y0,…,yNy_0,\dots,y_Nはすべて定まり、0≤n≤N0\le n\le NについてEn≤BnE_n\le B_nが成り立つ。

証明.u=y(tn)u=y(t_n)に仮定を用いると(tn,y(tn),hn)∈D(t_n,y(t_n),h_n)\in Dであるから、dnd_nは定まり、y(tn+1)=Ψhn(tn,y(tn))+dny(t_{n+1})=\Psi_{h_n}(t_n,y(t_n))+d_nである。m∈{0,…,N}m\in\{0,\dots,N\}について、y0,…,ymy_0,\dots,y_mが定まり、0≤n<m0\le n<mについて

En+1≤(1+Λhn)En+∥dn∥+∥rn∥E_{n+1}\le(1+\Lambda h_n)E_n+\|d_n\|+\|r_n\|

が成り立つという条件をP(m)P(m)とする。y0y_0は与えられているからP(0)P(0)が成り立つ。m<Nm<NについてP(m)P(m)が成り立つとする。補題 1.2 (1)をaj:=∥dj∥+∥rj∥a_j:=\|d_j\|+\|r_j\|として適用するとEm≤BmE_m\le B_mであり、補題 1.2 (2)によりBm≤BN≤ρB_m\le B_N\le\rhoである。仮定により(tm,ym,hm)∈D(t_m,y_m,h_m)\in Dであるからym+1y_{m+1}は定まり、

ym+1−y(tm+1)=Ψhm(tm,ym)−Ψhm(tm,y(tm))+rm−dmy_{m+1}-y(t_{m+1})=\Psi_{h_m}(t_m,y_m)-\Psi_{h_m}(t_m,y(t_m))+r_m-d_m

の両辺のノルムを取って仮定を用いると、n=mn=mの不等式を得る。したがってP(m+1)P(m+1)が成り立ち、帰納法によりP(N)P(N)が成り立つ。P(N)P(N)に補題 1.2 (1)を適用すると、0≤n≤N0\le n\le NについてEn≤BnE_n\le B_nである。▨

定義 1.4.d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。Ω⊆R×Rd\Omega\subseteq\R\times\R^dを開集合、f ⁣:Ω→Rdf\colon\Omega\to\R^dを連続写像とし、一段法がffに対応させる増分関数の定義域をDD、一歩写像をΨh\Psi_hとする。t0<Tt_0<Tとし、J⊇[t0,T]J\supseteq[t_0,T]を開区間、y ⁣:J→Rdy\colon J\to\R^dを、任意のτ∈J\tau\in Jについて(τ,y(τ))∈Ω(\tau,y(\tau))\in\Omegaとy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たすC1C^1級写像とする。

  1. t∈[t0,T)t\in[t_0,T)と0<h≤T−t0<h\le T-tが(t,y(t),h)∈D(t,y(t),h)\in Dを満たすとき、局所打切り誤差δy(t,h)=y(t+h)−Ψh(t,y(t))\delta_y(t,h)=y(t+h)-\Psi_h(t,y(t))をhhで割ったδy(t,h)/h\delta_y(t,h)/hを 正規化局所打切り誤差 (normalized local truncation error) という。
  2. H1>0H_1>0が存在して、t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{H1,T−t}0<h\le\min\{H_1,T-t\}を満たす任意のt,ht,hについて(t,y(t),h)∈D(t,y(t),h)\in Dであり、0<H≤H10<H\le H_1に対して τy(H):=sup⁡{∥δy(t,h)∥h ∣ t∈[t0,T), 0<h≤min⁡{H,T−t}}\tau_y(H):=\sup\Bigl\{\frac{\|\delta_y(t,h)\|}{h}\ \Big|\ t\in[t_0,T),\ 0<h\le\min\{H,T-t\}\Bigr\} と置くとlim⁡H→+0τy(H)=0\lim_{H\to+0}\tau_y(H)=0が成り立つとき、一段法は解yyと整合的であるという。すなわち、t∈[t0,T)t\in[t_0,T)について一様にδy(t,h)=o(h)\delta_y(t,h)=o(h)（h→+0h\to+0）であることをいう。任意のd∈N≥1d\in\NN、Rd\R^dの任意のノルム、任意の開集合Ω⊆R×Rd\Omega\subseteq\R\times\R^d、一段法のクラスに属する任意のf∈C1(Ω;Rd)f\in C^1(\Omega;\R^d)と、上の条件を満たす任意のt0<Tt_0<T、JJ、yyについて、一段法がyyと整合的であるとき、一段法は 整合的 (consistent) であるという。
  3. p∈N≥1p\in\NNとする。一段法の局所次数がpp以上であるとき、一段法の 整合性の次数 (order of consistency) はpp以上であるともいう。局所次数の定義の評価∥δy(t,h)∥≤Chp+1\|\delta_y(t,h)\|\le Ch^{p+1}は、正規化局所打切り誤差について∥δy(t,h)/h∥≤Chp\|\delta_y(t,h)/h\|\le Ch^pである。

注意 1.5. 局所次数が11以上の一段法は整合的である。実際、f∈C1f\in C^1の解yyについて、局所次数の定義のCCとh0h_0を取ると、t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{h0,T−t}0<h\le\min\{h_0,T-t\}について(t,y(t),h)∈D(t,y(t),h)\in Dであり、0<H≤h00<H\le h_0で定義 1.4 (2)の上限はτy(H)≤CH\tau_y(H)\le CHを満たす。

c∈Rc\in\Rとし、各ffにD:=Ω×(0,∞)D:=\Omega\times(0,\infty)上の増分関数Φ(t,u,h):=cf(t,u)\Phi(t,u,h):=cf(t,u)を対応させる一段法を考える。d=1d=1、Ω=R2\Omega=\R^2、f=0f=0の解は定数であり、任意のccについてδy(t,h)=0\delta_y(t,h)=0である。f(t,u)=1f(t,u)=1と解y(τ)=τy(\tau)=\tauについてはδy(t,h)=(1−c)h\delta_y(t,h)=(1-c)hであり、正規化局所打切り誤差は1−c1-cである。したがってc≠1c\ne1のときこの一段法は整合的でない。c=1c=1の一段法は前進 Euler 法であり、§E24.27 系 5.2 (1)により局所次数が11であるから整合的である。

系 1.6.d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。Ω⊆R×Rd\Omega\subseteq\R\times\R^dを開集合、f ⁣:Ω→Rdf\colon\Omega\to\R^dを連続写像とし、一段法がffに対応させる増分関数の定義域をDD、一歩写像をΨh\Psi_hとする。t0<Tt_0<Tとし、J⊇[t0,T]J\supseteq[t_0,T]を開区間、y ⁣:J→Rdy\colon J\to\R^dを、任意のτ∈J\tau\in Jについて(τ,y(τ))∈Ω(\tau,y(\tau))\in\Omegaとy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たすC1C^1級写像とする。ρ∈(0,∞]\rho\in(0,\infty]、Λ≥0\Lambda\ge0、H0>0H_0>0を取り、t∈[t0,T)t\in[t_0,T)、0<h≤min⁡{H0,T−t}0<h\le\min\{H_0,T-t\}と、∥u−y(t)∥≤ρ\|u-y(t)\|\le\rhoを満たす任意のu∈Rdu\in\R^dについて

(t,u,h)∈D,∥Ψh(t,u)−Ψh(t,y(t))∥≤(1+Λh)∥u−y(t)∥(t,u,h)\in D,\qquad\|\Psi_h(t,u)-\Psi_h(t,y(t))\|\le(1+\Lambda h)\|u-y(t)\|

が成り立つとする。Λ>0\Lambda>0ならばφ:=(eΛ(T−t0)−1)/Λ\varphi:=(e^{\Lambda(T-t_0)}-1)/\Lambda、Λ=0\Lambda=0ならばφ:=T−t0\varphi:=T-t_0と置き、0<H≤H00<H\le H_0に対してτy(H)\tau_y(H)を定義 1.4 (2)の上限とする。

  1. [t0,T][t_0,T]の格子t0<⋯<tN=Tt_0<\dots<t_N=TがH:=max⁡nhn≤H0H:=\max_nh_n\le H_0を満たすとする。y0∈Rdy_0\in\R^d、r0,…,rN−1∈Rdr_0,\dots,r_{N-1}\in\R^dを取り、(tn,yn,hn)∈D(t_n,y_n,h_n)\in Dである限りyn+1:=Ψhn(tn,yn)+rny_{n+1}:=\Psi_{h_n}(t_n,y_n)+r_nと定め、En:=∥yn−y(tn)∥E_n:=\|y_n-y(t_n)\|、ε:=max⁡n∥rn∥/hn\varepsilon:=\max_n\|r_n\|/h_nと置く。 eΛ(T−t0)E0+φ(τy(H)+ε)≤ρe^{\Lambda(T-t_0)}E_0+\varphi\bigl(\tau_y(H)+\varepsilon\bigr)\le\rho ならば、y0,…,yNy_0,\dots,y_Nはすべて定まり、max⁡nEn\max_nE_nは左辺以下である。
  2. 一段法が解yyと整合的ならば、一段法はyyに対して収束する。

証明.主張 (1)を示す。dn:=δy(tn,hn)d_n:=\delta_y(t_n,h_n)は∥dn∥≤τy(H)hn\|d_n\|\le\tau_y(H)h_nを満たし、∥rn∥≤εhn\|r_n\|\le\varepsilon h_nである。定理 1.3のBnB_nについて、補題 1.2 (3)をα:=τy(H)+ε\alpha:=\tau_y(H)+\varepsilonとして適用すると、eΛ(tN−t0)=eΛ(T−t0)e^{\Lambda(t_N-t_0)}=e^{\Lambda(T-t_0)}とφΛ(T−t0)=φ\varphi_\Lambda(T-t_0)=\varphiによりBNB_Nは仮定の左辺以下であり、したがってBN≤ρB_N\le\rhoである。定理 1.3と補題 1.2 (2)によりy0,…,yNy_0,\dots,y_Nはすべて定まり、En≤Bn≤BNE_n\le B_n\le B_Nである。

主張 (2)を示す。τy\tau_yはHHについて単調非減少であり、τy(H)→0\tau_y(H)\to0であるから、H1∈(0,H0]H_1\in(0,H_0]が存在してH≤H1H\le H_1ならばφτy(H)≤ρ\varphi\tau_y(H)\le\rhoである。H≤H1H\le H_1を満たす格子について、y0=y(t0)y_0=y(t_0)、rn=0r_n=0として主張 (1)を適用すると、近似値はすべて定まりE≤φτy(H)E\le\varphi\tau_y(H)である。β>0\beta>0に対してφτy(H2)≤β\varphi\tau_y(H_2)\le\betaを満たすH2∈(0,H1]H_2\in(0,H_1]を取ると、H≤H2H\le H_2を満たす任意の格子についてE≤βE\le\betaである。▨

系 1.7.系 1.6の仮定に加えて、p∈N≥1p\in\NNとC≥0C\ge0が存在し、t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{H0,T−t}0<h\le\min\{H_0,T-t\}を満たす任意のt,ht,hについて∥δy(t,h)∥≤Chp+1\|\delta_y(t,h)\|\le Ch^{p+1}が成り立つとする。φ\varphiを系 1.6の定数とする。H≤H0H\le H_0を満たす[t0,T][t_0,T]の格子、y0y_0、rnr_n、EnE_n、ε\varepsilonを系 1.6 (1)のとおりとすると、

eΛ(T−t0)E0+φ(CHp+ε)≤ρe^{\Lambda(T-t_0)}E_0+\varphi\bigl(CH^p+\varepsilon\bigr)\le\rho

ならばy0,…,yNy_0,\dots,y_Nはすべて定まり、max⁡nEn\max_nE_nは左辺以下である。特に、y0=y(t0)y_0=y(t_0)、rn=0r_n=0、φCHp≤ρ\varphi CH^p\le\rhoならば

E≤C eΛ(T−t0)−1Λ Hp(Λ=0 のときは E≤C(T−t0)Hp)E\le C\,\frac{e^{\Lambda(T-t_0)}-1}{\Lambda}\,H^p\qquad(\Lambda=0\ \text{のときは}\ E\le C(T-t_0)H^p)

であり、一段法はyyに対してpp次で収束する。

証明.0<H≤H00<H\le H_0についてτy(H)≤CHp\tau_y(H)\le CH^pであるから、系 1.6 (1)により第一の主張を得る。C>0C>0かつρ<∞\rho<\inftyならばH1:=min⁡{H0,(ρ/(Cφ))1/p}H_1:=\min\{H_0,(\rho/(C\varphi))^{1/p}\}、そうでなければH1:=H0H_1:=H_0と置くと、H≤H1H\le H_1を満たす格子でφCHp≤ρ\varphi CH^p\le\rhoである。y0=y(t0)y_0=y(t_0)、rn=0r_n=0として第一の主張を適用するとE≤CφHpE\le C\varphi H^pである。▨

注意 1.8. 等間隔の格子hn=h=(T−t0)/Nh_n=h=(T-t_0)/Nで∥rn∥≤ε′\|r_n\|\le\varepsilon'とすると、系 1.7のε\varepsilonはε′/h=Nε′/(T−t0)\varepsilon'/h=N\varepsilon'/(T-t_0)以下であり、上界eΛ(T−t0)E0+φ(Chp+ε′/h)e^{\Lambda(T-t_0)}E_0+\varphi(Ch^p+\varepsilon'/h)の最後の項は歩数NNに比例する。d=1d=1、f=0f=0、解y=0y=0に前進 Euler 法を用いるとΨh(t,u)=u\Psi_h(t,u)=uであり、y0=0y_0=0、rn=ε′r_n=\varepsilon'からyN−y(T)=Nε′y_N-y(T)=N\varepsilon'を得る。C>0C>0とε′>0\varepsilon'>0について、h↦Chp+ε′/hh\mapsto Ch^p+\varepsilon'/hはh>0h>0上でh=(ε′/(pC))1/(p+1)h=(\varepsilon'/(pC))^{1/(p+1)}において最小値をとる。

2 θ 法

定義 2.1.θ∈[0,1]\theta\in[0,1]とする。d∈N≥1d\in\NNとし、Rd\R^dにノルムを固定する。Ω⊆R×Rd\Omega\subseteq\R\times\R^dを開集合、f ⁣:Ω→Rdf\colon\Omega\to\R^dを連続写像、(t,u)∈Ω(t,u)\in\Omega、h>0h>0とし、v∈Rdv\in\R^dについての方程式

v=u+h((1−θ)f(t,u)+θf(t+h,v)),(t+h,v)∈Ωv=u+h\bigl((1-\theta)f(t,u)+\theta f(t+h,v)\bigr),\qquad(t+h,v)\in\Omega

がただ一つの解をもつとき、その解をΨh(t,u)\Psi_h(t,u)とする。Ψh(t,u)\Psi_h(t,u)が定まる(t,u,h)(t,u,h)の全体をDDとし、増分関数Φ(t,u,h):=(Ψh(t,u)−u)/h\Phi(t,u,h):=(\Psi_h(t,u)-u)/hを対応させる一段法を θ 法 (theta method) という。θ=12\theta=\frac12の θ 法、すなわち方程式

v=u+h2(f(t,u)+f(t+h,v))v=u+\frac h2\bigl(f(t,u)+f(t+h,v)\bigr)

による一段法を 台形法 (trapezoidal method) という。θ=0\theta=0の θ 法の一歩写像は、定まる点で前進 Euler 法の一歩写像に一致し、θ=1\theta=1の θ 法の一歩写像は、定まる点で§E24.27 定義 3.1 (1)により後退 Euler 法の一歩写像に一致する。

補題 2.2.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について、任意のσ∈I\sigma\in Iとv,v′∈Rdv,v'\in\R^dに対して∥f(σ,v)−f(σ,v′)∥≤L∥v−v′∥\|f(\sigma,v)-f(\sigma,v')\|\le L\|v-v'\|が成り立つとする。θ∈[0,1]\theta\in[0,1]、t∈It\in I、h>0h>0、t+h∈It+h\in I、θhL<1\theta hL<1とし、Ψh\Psi_hを θ 法の一歩写像とする。

  1. 任意のu∈Rdu\in\R^dについて、方程式v=u+h((1−θ)f(t,u)+θf(t+h,v))v=u+h((1-\theta)f(t,u)+\theta f(t+h,v))はRd\R^dにただ一つの解Ψh(t,u)\Psi_h(t,u)をもつ。任意のu,w∈Rdu,w\in\R^dについて ∥Ψh(t,u)−Ψh(t,w)∥≤(1+hL1−θhL)∥u−w∥\|\Psi_h(t,u)-\Psi_h(t,w)\|\le\Bigl(1+\frac{hL}{1-\theta hL}\Bigr)\|u-w\| が成り立つ。
  2. yyを、ttとt+ht+hを含む開区間上でy′(σ)=f(σ,y(σ))y'(\sigma)=f(\sigma,y(\sigma))を満たすC1C^1級写像とし、η:=y(t+h)−y(t)−h((1−θ)y′(t)+θy′(t+h))\eta:=y(t+h)-y(t)-h((1-\theta)y'(t)+\theta y'(t+h))と置く。このとき ∥y(t+h)−Ψh(t,y(t))∥≤∥η∥1−θhL,∥y(t+h)−Ψh(t,y(t))−η∥≤θhL∥η∥1−θhL\|y(t+h)-\Psi_h(t,y(t))\|\le\frac{\|\eta\|}{1-\theta hL},\qquad\|y(t+h)-\Psi_h(t,y(t))-\eta\|\le\frac{\theta hL\|\eta\|}{1-\theta hL} が成り立つ。

証明.主張 (1)を示す。u∈Rdu\in\R^dを取り、T(v):=u+h(1−θ)f(t,u)+hθf(t+h,v)T(v):=u+h(1-\theta)f(t,u)+h\theta f(t+h,v)と置く。TTはRd\R^dからそれ自身への写像であり、∥T(v)−T(v′)∥≤θhL∥v−v′∥\|T(v)-T(v')\|\le\theta hL\|v-v'\|を満たす。Rd\R^dは空でない完備距離空間でありθhL<1\theta hL<1であるから、§E2.7 定理 2.2によりTTはただ一つの不動点をもつ。(t+h,v)∈I×Rd(t+h,v)\in I\times\R^dは常に成り立つから、方程式の解はTTの不動点と一致し、ただ一つである。v:=Ψh(t,u)v:=\Psi_h(t,u)、v′:=Ψh(t,w)v':=\Psi_h(t,w)と置くと

v−v′=u−w+h(1−θ)(f(t,u)−f(t,w))+hθ(f(t+h,v)−f(t+h,v′))v-v'=u-w+h(1-\theta)\bigl(f(t,u)-f(t,w)\bigr)+h\theta\bigl(f(t+h,v)-f(t+h,v')\bigr)

であるから∥v−v′∥≤(1+(1−θ)hL)∥u−w∥+θhL∥v−v′∥\|v-v'\|\le(1+(1-\theta)hL)\|u-w\|+\theta hL\|v-v'\|である。移項すると∥v−v′∥≤1+(1−θ)hL1−θhL∥u−w∥\|v-v'\|\le\frac{1+(1-\theta)hL}{1-\theta hL}\|u-w\|であり、1+(1−θ)hL1−θhL=1+hL1−θhL\frac{1+(1-\theta)hL}{1-\theta hL}=1+\frac{hL}{1-\theta hL}である。

主張 (2)を示す。v:=Ψh(t,y(t))v:=\Psi_h(t,y(t))と置く。y′(t)=f(t,y(t))y'(t)=f(t,y(t))、y′(t+h)=f(t+h,y(t+h))y'(t+h)=f(t+h,y(t+h))とvvの方程式から

y(t+h)−v=η+hθ(f(t+h,y(t+h))−f(t+h,v))y(t+h)-v=\eta+h\theta\bigl(f(t+h,y(t+h))-f(t+h,v)\bigr)

である。したがって∥y(t+h)−v∥≤∥η∥+θhL∥y(t+h)−v∥\|y(t+h)-v\|\le\|\eta\|+\theta hL\|y(t+h)-v\|であり、移項すると第一の評価を得る。∥y(t+h)−v−η∥≤θhL∥y(t+h)−v∥\|y(t+h)-v-\eta\|\le\theta hL\|y(t+h)-v\|に第一の評価を代入すると第二の評価を得る。▨

補題 2.3.m∈N≥1m\in\NNとし、Rm\R^mにノルム∥⋅∥\|\cdot\|を固定する。J⊆RJ\subseteq\Rを開区間、θ∈[0,1]\theta\in[0,1]、t∈Jt\in J、h>0h>0、t+h∈Jt+h\in Jとし、C1C^1級写像x ⁣:J→Rmx\colon J\to\R^mに対して

ηθ:=x(t+h)−x(t)−h((1−θ)x′(t)+θx′(t+h))\eta_\theta:=x(t+h)-x(t)-h\bigl((1-\theta)x'(t)+\theta x'(t+h)\bigr)

と置く。k∈N≥0k\in\Nについてx∈Ck(J;Rm)x\in C^k(J;\R^m)のときMk:=max⁡τ∈[t,t+h]∥x(k)(τ)∥M_k:=\max_{\tau\in[t,t+h]}\|x^{(k)}(\tau)\|と置く。

  1. x∈C2(J;Rm)x\in C^2(J;\R^m)ならば∥ηθ∥≤h22M2\|\eta_\theta\|\le\frac{h^2}2M_2である。
  2. x∈C3(J;Rm)x\in C^3(J;\R^m)ならば ∥ηθ−(12−θ)h2x′′(t)∥≤(16+θ2)h3M3\Bigl\|\eta_\theta-\Bigl(\frac12-\theta\Bigr)h^2x''(t)\Bigr\|\le\Bigl(\frac16+\frac\theta2\Bigr)h^3M_3 である。特に∥η1/2∥≤512h3M3\|\eta_{1/2}\|\le\frac5{12}h^3M_3である。
  3. x∈C4(J;Rm)x\in C^4(J;\R^m)ならば ∥η1/2+h312x′′′(t)∥≤h48M4\Bigl\|\eta_{1/2}+\frac{h^3}{12}x'''(t)\Bigr\|\le\frac{h^4}8M_4 である。

証明.主張 (1)を示す。

ηθ=(1−θ)(x(t+h)−x(t)−hx′(t))−θ(x(t)−x(t+h)−(−h)x′(t+h))\eta_\theta=(1-\theta)\bigl(x(t+h)-x(t)-hx'(t)\bigr)-\theta\bigl(x(t)-x(t+h)-(-h)x'(t+h)\bigr)

である。§E24.27 補題 1.5 (1)をr=1r=1として、展開点tt、増分hhの場合と、展開点t+ht+h、増分−h-hの場合に適用すると、二つの括弧のノルムはともにh22M2\frac{h^2}2M_2以下であり、(1−θ)+θ=1(1-\theta)+\theta=1により主張を得る。

主張 (2)を示す。

ηθ−(12−θ)h2x′′(t)=(x(t+h)−x(t)−hx′(t)−h22x′′(t))−θh(x′(t+h)−x′(t)−hx′′(t))\eta_\theta-\Bigl(\frac12-\theta\Bigr)h^2x''(t)=\Bigl(x(t+h)-x(t)-hx'(t)-\frac{h^2}2x''(t)\Bigr)-\theta h\bigl(x'(t+h)-x'(t)-hx''(t)\bigr)

である。§E24.27 補題 1.5 (1)をxxとr=2r=2に適用すると第一の括弧のノルムはh36M3\frac{h^3}6M_3以下であり、x′x'とr=1r=1に適用すると第二の括弧のノルムはh22M3\frac{h^2}2M_3以下である。θ=12\theta=\frac12では左辺はη1/2\eta_{1/2}であり、16+14=512\frac16+\frac14=\frac5{12}である。

主張 (3)を示す。

R1:=x(t+h)−∑j=03hjj!x(j)(t),R2:=x′(t+h)−∑j=02hjj!x(j+1)(t)R_1:=x(t+h)-\sum_{j=0}^3\frac{h^j}{j!}x^{(j)}(t),\qquad R_2:=x'(t+h)-\sum_{j=0}^2\frac{h^j}{j!}x^{(j+1)}(t)

と置くと、

∑j=13hjj!x(j)(t)−h2(x′(t)+∑j=02hjj!x(j+1)(t))=(16−14)h3x′′′(t)=−h312x′′′(t)\sum_{j=1}^3\frac{h^j}{j!}x^{(j)}(t)-\frac h2\Bigl(x'(t)+\sum_{j=0}^2\frac{h^j}{j!}x^{(j+1)}(t)\Bigr)=\Bigl(\frac16-\frac14\Bigr)h^3x'''(t)=-\frac{h^3}{12}x'''(t)

によりη1/2=R1−h2R2−h312x′′′(t)\eta_{1/2}=R_1-\frac h2R_2-\frac{h^3}{12}x'''(t)である。§E24.27 補題 1.5 (1)により∥R1∥≤h424M4\|R_1\|\le\frac{h^4}{24}M_4、∥R2∥≤h36M4\|R_2\|\le\frac{h^3}6M_4であり、124+112=18\frac1{24}+\frac1{12}=\frac18である。▨

定理 2.4.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について、任意のσ∈I\sigma\in Iとv,v′∈Rdv,v'\in\R^dに対して∥f(σ,v)−f(σ,v′)∥≤L∥v−v′∥\|f(\sigma,v)-f(\sigma,v')\|\le L\|v-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をy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たすC1C^1級写像とする。θ∈[0,1]\theta\in[0,1]と、θh0L<1\theta h_0L<1を満たすh0>0h_0>0を取り、

Λ:=L1−θh0L,φ:=eΛ(T−t0)−1Λ(Λ=0 のときは φ:=T−t0)\Lambda:=\frac L{1-\theta h_0L},\qquad\varphi:=\frac{e^{\Lambda(T-t_0)}-1}\Lambda\quad(\Lambda=0\ \text{のときは}\ \varphi:=T-t_0)

と置く。[t0,T][t_0,T]の格子がH≤h0H\le h_0を満たすとし、y0∈Rdy_0\in\R^d、r0,…,rN−1∈Rdr_0,\dots,r_{N-1}\in\R^dに対して、θ 法の一歩写像Ψh\Psi_hによりyn+1:=Ψhn(tn,yn)+rny_{n+1}:=\Psi_{h_n}(t_n,y_n)+r_nと定め、En:=∥yn−y(tn)∥E_n:=\|y_n-y(t_n)\|、ε:=max⁡n∥rn∥/hn\varepsilon:=\max_n\|r_n\|/h_nと置く。このときy0,…,yNy_0,\dots,y_Nはすべて定まり、次が成り立つ。

  1. y∈C2(J;Rd)y\in C^2(J;\R^d)ならば、M2:=max⁡τ∈[t0,T]∥y′′(τ)∥M_2:=\max_{\tau\in[t_0,T]}\|y''(\tau)\|について max⁡nEn≤eΛ(T−t0)E0+φ(M22(1−θh0L)H+ε)\max_nE_n\le e^{\Lambda(T-t_0)}E_0+\varphi\Bigl(\frac{M_2}{2(1-\theta h_0L)}H+\varepsilon\Bigr) である。
  2. θ=12\theta=\frac12かつy∈C3(J;Rd)y\in C^3(J;\R^d)ならば、M3:=max⁡τ∈[t0,T]∥y′′′(τ)∥M_3:=\max_{\tau\in[t_0,T]}\|y'''(\tau)\|について max⁡nEn≤eΛ(T−t0)E0+φ(5M312(1−h0L/2)H2+ε)\max_nE_n\le e^{\Lambda(T-t_0)}E_0+\varphi\Bigl(\frac{5M_3}{12(1-h_0L/2)}H^2+\varepsilon\Bigr) である。

特に、θ 法はC2C^2級の解yyに対して11次で収束し、台形法はC3C^3級の解yyに対して22次で収束する。

証明.t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{h0,T−t}0<h\le\min\{h_0,T-t\}を取る。t,t+h∈J⊆It,t+h\in J\subseteq IでありθhL≤θh0L<1\theta hL\le\theta h_0L<1であるから、補題 2.2 (1)により任意のu∈Rdu\in\R^dについて(t,u,h)∈D(t,u,h)\in Dであり、

∥Ψh(t,u)−Ψh(t,y(t))∥≤(1+hL1−θhL)∥u−y(t)∥≤(1+Λh)∥u−y(t)∥\|\Psi_h(t,u)-\Psi_h(t,y(t))\|\le\Bigl(1+\frac{hL}{1-\theta hL}\Bigr)\|u-y(t)\|\le(1+\Lambda h)\|u-y(t)\|

である。したがって系 1.6の仮定はρ=∞\rho=\infty、H0:=h0H_0:=h_0と上のΛ\Lambdaについて成り立つ。補題 2.2 (2)と補題 2.3 (1)により

∥δy(t,h)∥≤11−θh0L⋅h22M2\|\delta_y(t,h)\|\le\frac{1}{1-\theta h_0L}\cdot\frac{h^2}2M_2

であり、θ=12\theta=\frac12かつy∈C3y\in C^3ならば補題 2.2 (2)と補題 2.3 (2)により∥δy(t,h)∥≤5M312(1−h0L/2)h3\|\delta_y(t,h)\|\le\frac{5M_3}{12(1-h_0L/2)}h^3である。系 1.7をそれぞれp=1p=1とp=2p=2で適用すると、ρ=∞\rho=\inftyであるから条件は常に満たされ、二つの評価と収束の次数を得る。▨

注意 2.5.定理 2.4でθ=1\theta=1とし、h0L<1h_0L<1、ε′>0\varepsilon'>0とする。各nnでv(0):=ynv^{(0)}:=y_nから始めた反復v(k+1):=yn+hnf(tn+1,v(k))v^{(k+1)}:=y_n+h_nf(t_{n+1},v^{(k)})を∥v(k+1)−v(k)∥≤(1−hnL)hnε′\|v^{(k+1)}-v^{(k)}\|\le(1-h_nL)h_n\varepsilon'を満たすkkで止めてyn+1:=v(k)y_{n+1}:=v^{(k)}とすると、§E24.27 命題 3.4 (1)によりrn:=v(k)−Ψhn(tn,yn)r_n:=v^{(k)}-\Psi_{h_n}(t_n,y_n)は∥rn∥≤hnε′\|r_n\|\le h_n\varepsilon'を満たす。したがって定理 2.4の評価がε=ε′\varepsilon=\varepsilon'で成り立つ。

注意 2.6.d=1d=1、I=RI=\R、f(t,u)=λuf(t,u)=\lambda uとすると、λ=100\lambda=100とλ=−100\lambda=-100のどちらでもL=100L=100である。定理 2.4をθ=1\theta=1で適用するにはh0<1/100h_0<1/100が必要であり、評価にはeΛ(T−t0)≥e100(T−t0)e^{\Lambda(T-t_0)}\ge e^{100(T-t_0)}が現れる。この評価ではλ\lambdaの符号は用いられないので、解がe100te^{100t}で増大するλ=100\lambda=100の場合も含まれる。λ=−100\lambda=-100、y(0)=1y(0)=1、h=0.1h=0.1ではhL=10hL=10であり、定理 2.4は適用されない。一方、後退 Euler 法の方程式v=u−10vv=u-10vの解はv=u/11v=u/11であり、yn=11−ny_n=11^{-n}とy(tn)=e−10ny(t_n)=e^{-10n}はともに正であるから∣yn−y(tn)∣<11−n|y_n-y(t_n)|<11^{-n}である。λ\lambdaの符号を用いる評価は命題 5.3により与えられる。

命題 2.7.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について、任意のσ∈I\sigma\in Iとv,v′∈Rdv,v'\in\R^dに対して∥f(σ,v)−f(σ,v′)∥≤L∥v−v′∥\|f(\sigma,v)-f(\sigma,v')\|\le L\|v-v'\|が成り立つとする。t0<Tt_0<Tとし、JJを[t0,T]⊆J⊆I[t_0,T]\subseteq J\subseteq Iを満たす開区間、y∈C2(J;Rd)y\in C^2(J;\R^d)をy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たす写像、M2:=max⁡τ∈[t0,T]∥y′′(τ)∥M_2:=\max_{\tau\in[t_0,T]}\|y''(\tau)\|とする。[t0,T][t_0,T]の任意の格子について、初期値y0=y(t0)y_0=y(t_0)から前進 Euler 法で定まる近似値の大域誤差の大きさEEは、L>0L>0ならば

E≤M2H2⋅eL(T−t0)−1L,E\le\frac{M_2H}{2}\cdot\frac{e^{L(T-t_0)}-1}{L},

L=0L=0ならばE≤M2(T−t0)H/2E\le M_2(T-t_0)H/2を満たす。

証明.θ=0\theta=0の θ 法の一歩写像はΨh(t,u)=u+hf(t,u)\Psi_h(t,u)=u+hf(t,u)（t,t+h∈It,t+h\in I）であり、格子上で前進 Euler 法の一歩写像に一致する。h0:=T−t0h_0:=T-t_0は任意の格子についてH≤h0H\le h_0とθh0L=0<1\theta h_0L=0<1を満たし、Λ=L\Lambda=Lである。定理 2.4 (1)をθ=0\theta=0、E0=0E_0=0、rn=0r_n=0として適用すると主張を得る。▨

注意 2.8. 等間隔の格子hn=h=(T−t0)/Nh_n=h=(T-t_0)/Nと系 1.7の仮定の下で、N=(T−t0)/hN=(T-t_0)/h個の局所打切り誤差の和は∑j<N∥dj∥≤NChp+1=C(T−t0)hp\sum_{j<N}\|d_j\|\le NCh^{p+1}=C(T-t_0)h^pを満たす。各djd_jに掛かる係数eΛ(tn−tj+1)e^{\Lambda(t_n-t_{j+1})}はeΛ(T−t0)e^{\Lambda(T-t_0)}以下であり、hhによらない。一歩ごとの係数1+Λh1+\Lambda hを掛け合わせた和∑j<n(1+Λh)j\sum_{j<n}(1+\Lambda h)^jはnn以上であって、n=Nn=Nではh→+0h\to+0で有界でない。有界であるのはこれにhhを掛けた量であり、Λ>0\Lambda>0ならば

h∑j=0n−1(1+Λh)j=(1+Λh)n−1Λ≤eΛ(tn−t0)−1Λh\sum_{j=0}^{n-1}(1+\Lambda h)^j=\frac{(1+\Lambda h)^n-1}{\Lambda}\le\frac{e^{\Lambda(t_n-t_0)}-1}{\Lambda}

である。

例 2.9.d=1d=1、I=RI=\R、f(t,u)=uf(t,u)=u（L=1L=1）、解y(τ)=eτy(\tau)=e^\tau、t0=0t_0=0、T=1T=1、N∈N≥1N\in\NN、h=1/Nh=1/Nとする。前進 Euler 法の近似値はyn=(1+h)ny_n=(1+h)^nであり、終点ではyN=(1+h)Ny_N=(1+h)^Nである。a:=eha:=e^h、b:=1+hb:=1+hと置くと

e−yN=aN−bN=(a−b)∑k=0N−1aN−1−kbke-y_N=a^N-b^N=(a-b)\sum_{k=0}^{N-1}a^{N-1-k}b^k

である。§E24.27 例 1.7によりh22≤a−b≤h22eh\frac{h^2}2\le a-b\le\frac{h^2}2e^hであり、1≤b≤a1\le b\le aによりN≤∑k=0N−1aN−1−kbk≤NaN−1=Ne1−hN\le\sum_{k=0}^{N-1}a^{N-1-k}b^k\le Na^{N-1}=Ne^{1-h}である。Nh=1Nh=1を用いると

h2≤e−yN≤e2h\frac h2\le e-y_N\le\frac e2h

を得る。命題 2.7によりM2=eM_2=eとしてE≤e(e−1)2hE\le\frac{e(e-1)}2hが与えられ、e(e−1)2=2.335387135…\frac{e(e-1)}2=2.335387135\ldots、e2=1.359140914…\frac e2=1.359140914\ldotsである。値を丸めると次のとおりである。

NN yNy_N e−yNe-y_N (e−yN)/h(e-y_N)/h 一つ上の行のe−yNe-y_Nとの比
22 2.2500000000002.250000000000 0.4682818284590.468281828459 0.9365636570.936563657 —
44 2.4414062500002.441406250000 0.2768755784590.276875578459 1.1075023141.107502314 1.6913078111.691307811
88 2.5657845139502.565784513950 0.1524973145090.152497314509 1.2199785161.219978516 1.8156095361.815609536
1616 2.6379284973672.637928497367 0.0803533310920.080353331092 1.2856532971.285653297 1.8978343831.897834383
3232 2.6769901293782.676990129378 0.0412916990810.041291699081 1.3213343711.321334371 1.9459923641.945992364

3 陽的 Runge–Kutta 法

補題 3.1.s∈N≥1s\in\NNとし、A=(aij)∈Rs×sA=(a_{ij})\in\R^{s\times s}を狭義下三角行列、b∈Rsb\in\R^s、c:=A1c:=A\mathbf1とする。d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。Ω⊆R×Rd\Omega\subseteq\R\times\R^dを開集合、f ⁣:Ω→Rdf\colon\Omega\to\R^dを連続写像とし、Φ\PhiとΨh\Psi_hを Butcher 配列(A,b)(A,b)の陽的 Runge–Kutta 法の増分関数と一歩写像とする。L≥0L\ge0とh0>0h_0>0に対して

L1:=L,Li:=L(1+h0∑j<i∣aij∣Lj)(2≤i≤s),Λ:=∑i=1s∣bi∣LiL_1:=L,\qquad L_i:=L\Bigl(1+h_0\sum_{j<i}|a_{ij}|L_j\Bigr)\quad(2\le i\le s),\qquad\Lambda:=\sum_{i=1}^s|b_i|L_i

と置く。B⊆ΩB\subseteq\Omegaが、(σ,w),(σ,w′)∈B(\sigma,w),(\sigma,w')\in Bを満たす任意のσ,w,w′\sigma,w,w'について∥f(σ,w)−f(σ,w′)∥≤L∥w−w′∥\|f(\sigma,w)-f(\sigma,w')\|\le L\|w-w'\|を満たすとする。0<h≤h00<h\le h_0、t∈Rt\in\R、u,v∈Rdu,v\in\R^dとし、k1:=f(t,u)k_1:=f(t,u)、ki:=f(t+cih, u+h∑j<iaijkj)k_i:=f(t+c_ih,\ u+h\sum_{j<i}a_{ij}k_j)で順に定める値の引数と、uuをvvに替えて定める値ki′k'_iの引数が、すべて定まりBBに属するとする。このとき(t,u,h),(t,v,h)∈D(t,u,h),(t,v,h)\in Dであり、

∥Φ(t,u,h)−Φ(t,v,h)∥≤Λ∥u−v∥,∥Ψh(t,u)−Ψh(t,v)∥≤(1+Λh)∥u−v∥\|\Phi(t,u,h)-\Phi(t,v,h)\|\le\Lambda\|u-v\|,\qquad\|\Psi_h(t,u)-\Psi_h(t,v)\|\le(1+\Lambda h)\|u-v\|

が成り立つ。

証明. 引数はすべてB⊆ΩB\subseteq\Omegaに属するから、§E24.27 定義 2.2 (3)により(t,u,h),(t,v,h)∈D(t,u,h),(t,v,h)\in Dである。k1k_1とk1′k'_1の引数(t,u)(t,u)と(t,v)(t,v)は時刻が等しくBBに属するから、∥k1−k1′∥≤L∥u−v∥=L1∥u−v∥\|k_1-k'_1\|\le L\|u-v\|=L_1\|u-v\|である。i≥2i\ge2とし、j<ij<iについて∥kj−kj′∥≤Lj∥u−v∥\|k_j-k'_j\|\le L_j\|u-v\|であるとする。kik_iとki′k'_iの引数は時刻t+ciht+c_ihが等しくBBに属するから

∥ki−ki′∥≤L∥u−v+h∑j<iaij(kj−kj′)∥≤L(1+h0∑j<i∣aij∣Lj)∥u−v∥=Li∥u−v∥\|k_i-k'_i\|\le L\Bigl\|u-v+h\sum_{j<i}a_{ij}(k_j-k'_j)\Bigr\|\le L\Bigl(1+h_0\sum_{j<i}|a_{ij}|L_j\Bigr)\|u-v\|=L_i\|u-v\|

である。帰納法により1≤i≤s1\le i\le sについて∥ki−ki′∥≤Li∥u−v∥\|k_i-k'_i\|\le L_i\|u-v\|である。Φ(t,u,h)−Φ(t,v,h)=∑ibi(ki−ki′)\Phi(t,u,h)-\Phi(t,v,h)=\sum_ib_i(k_i-k'_i)であるから第一の評価を得て、Ψh(t,u)=u+hΦ(t,u,h)\Psi_h(t,u)=u+h\Phi(t,u,h)から第二の評価を得る。▨

系 3.2.s,p∈N≥1s,p\in\NNとし、Butcher 配列(A,b)(A,b)の陽的 Runge–Kutta 法の局所次数がpp以上であるとする。d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。Ω⊆R×Rd\Omega\subseteq\R\times\R^dを開集合、f∈Cp(Ω;Rd)f\in C^p(\Omega;\R^d)、t0<Tt_0<Tとし、J⊇[t0,T]J\supseteq[t_0,T]を開区間、y ⁣:J→Rdy\colon J\to\R^dを、任意のτ∈J\tau\in Jについて(τ,y(τ))∈Ω(\tau,y(\tau))\in\Omegaとy′(τ)=f(τ,y(τ))y'(\tau)=f(\tau,y(\tau))を満たすC1C^1級写像とする。このとき、ρ>0\rho>0、Λ≥0\Lambda\ge0、H0>0H_0>0が存在して系 1.6の仮定が成り立ち、この Runge–Kutta 法はyyに対してpp次で収束する。特に、§E24.27 系 5.2により、解yyに対して、前進 Euler 法はf∈C1(Ω;Rd)f\in C^1(\Omega;\R^d)ならば11次で、中点法と Heun 法はf∈C2(Ω;Rd)f\in C^2(\Omega;\R^d)ならば22次で、古典的 Runge–Kutta 法はf∈C4(Ω;Rd)f\in C^4(\Omega;\R^d)ならば44次で収束する。

証明.R×Rd\R\times\R^dにノルムmax⁡{∣t∣,∥u∥}\max\{|t|,\|u\|\}を入れる。Γ:={(τ,y(τ))∣τ∈[t0,T]}\Gamma:=\{(\tau,y(\tau))\mid\tau\in[t_0,T]\}は空でないコンパクト集合であり、§E24.27 補題 1.4によりΓρ′⊆Ω\Gamma^{\rho'}\subseteq\Omegaを満たすρ′>0\rho'>0が存在し、Γρ′\Gamma^{\rho'}はコンパクトである。M:=max⁡Γρ′∥f∥M:=\max_{\Gamma^{\rho'}}\|f\|、α:=max⁡i∑j∣aij∣\alpha:=\max_i\sum_j|a_{ij}|、γ:=max⁡i∣ci∣\gamma:=\max_i|c_i|と置く。p≥1p\ge1によりf∈C1(Ω;Rd)f\in C^1(\Omega;\R^d)であるから、§E24.27 補題 2.1 (1)をK=ΓK=\Gamma、ρ=ρ′\rho=\rho'に適用してL≥0L\ge0を得る。局所次数の定義により、C≥0C\ge0とh1>0h_1>0が存在して、t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{h1,T−t}0<h\le\min\{h_1,T-t\}について∥δy(t,h)∥≤Chp+1\|\delta_y(t,h)\|\le Ch^{p+1}である。H0∈(0,h1]H_0\in(0,h_1]をH0γ≤ρ′H_0\gamma\le\rho'とH0αM≤ρ′/2H_0\alpha M\le\rho'/2を満たすように取り、ρ:=ρ′/2\rho:=\rho'/2と置く。

t∈[t0,T)t\in[t_0,T)、0<h≤min⁡{H0,T−t}0<h\le\min\{H_0,T-t\}、∥u−y(t)∥≤ρ\|u-y(t)\|\le\rhoとし、Bt:={(σ,w)∣∣σ−t∣≤ρ′, ∥w−y(t)∥≤ρ′}B_t:=\{(\sigma,w)\mid|\sigma-t|\le\rho',\ \|w-y(t)\|\le\rho'\}と置くとBt⊆Γρ′⊆ΩB_t\subseteq\Gamma^{\rho'}\subseteq\Omegaである。k1,…,ki−1k_1,\dots,k_{i-1}が定まり∥kj∥≤M\|k_j\|\le Mを満たすとすると、∣cih∣≤γH0≤ρ′|c_ih|\le\gamma H_0\le\rho'と

∥u+h∑j<iaijkj−y(t)∥≤ρ+hαM≤ρ′\Bigl\|u+h\sum_{j<i}a_{ij}k_j-y(t)\Bigr\|\le\rho+h\alpha M\le\rho'

によりkik_iの引数はBtB_tに属し、kik_iは定まって∥ki∥≤M\|k_i\|\le Mである。iiについての帰納法により、k1,…,ksk_1,\dots,k_sの引数はすべてBtB_tに属する。u=y(t)u=y(t)の場合も同じである。§E24.27 補題 2.1 (1)により、(σ,w),(σ,w′)∈Bt(\sigma,w),(\sigma,w')\in B_tならば∥f(σ,w)−f(σ,w′)∥≤L∥w−w′∥\|f(\sigma,w)-f(\sigma,w')\|\le L\|w-w'\|である。補題 3.1をB=BtB=B_t、h0=H0h_0=H_0、v=y(t)v=y(t)として適用すると、(t,u,h)∈D(t,u,h)\in Dと∥Ψh(t,u)−Ψh(t,y(t))∥≤(1+Λh)∥u−y(t)∥\|\Psi_h(t,u)-\Psi_h(t,y(t))\|\le(1+\Lambda h)\|u-y(t)\|を得る。ここでΛ\LambdaはLL、H0H_0、(A,b)(A,b)だけで定まる。系 1.7により、この一段法はyyに対してpp次で収束する。▨

例 3.3.d=1d=1、Ω=R2\Omega=\R^2、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とする。古典的 Runge–Kutta 法の段はk1=uk_1=u、k2=(1+h2)uk_2=(1+\frac h2)u、k3=(1+h2+h24)uk_3=(1+\frac h2+\frac{h^2}4)u、k4=(1+h+h22+h34)uk_4=(1+h+\frac{h^2}2+\frac{h^3}4)uであり、

Ψh(t,u)=R(h)u,R(h):=1+h+h22+h36+h424\Psi_h(t,u)=R(h)u,\qquad R(h):=1+h+\frac{h^2}2+\frac{h^3}6+\frac{h^4}{24}

である。したがってyN=R(h)Ny_N=R(h)^Nである。a:=eha:=e^h、b:=R(h)b:=R(h)と置くとa−b=∑j≥5hj/j!≥h5/120a-b=\sum_{j\ge5}h^j/j!\ge h^5/120であり、§E24.27 補題 1.5 (1)をexp⁡\expとr=4r=4に適用するとa−b≤h5eh/120a-b\le h^5e^h/120である。1≤b≤a1\le b\le aと

e−yN=(a−b)∑k=0N−1aN−1−kbke-y_N=(a-b)\sum_{k=0}^{N-1}a^{N-1-k}b^k

からh4120≤e−yN≤e120h4\frac{h^4}{120}\le e-y_N\le\frac{e}{120}h^4を得る。1120=0.008333333…\frac1{120}=0.008333333\ldots、e120=0.022652348570…\frac e{120}=0.022652348570\ldotsである。値を丸めると次のとおりである。

NN R(h)R(h) yNy_N e−yNe-y_N (e−yN)/h4(e-y_N)/h^4 一つ上の行のe−yNe-y_Nとの比
22 1.6484375000001.648437500000 2.7173461914062.717346191406 0.0009356370530.000935637053 0.0149701930.014970193 —
44 1.2840169270831.284016927083 2.7182099392012.718209939201 0.0000718892580.000071889258 0.0184036500.018403650 13.01497713.014977
88 1.1331481933591.133148193359 2.7182768444172.718276844417 0.0000049840420.000004984042 0.0204146370.020414637 14.42388614.423886
1616 1.0644944508871.064494450887 2.7182815003412.718281500341 0.0000003281180.000000328118 0.0215035710.021503571 15.18976515.189765

4 絶対連続な補間の欠陥

定理 4.1.d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。t0<Tt_0<T、S⊆[t0,T]×RdS\subseteq[t_0,T]\times\R^dとし、f ⁣:S→Rdf\colon S\to\R^dを写像とする。非負関数ℓ∈L1([t0,T])\ell\in L^1([t_0,T])と Lebesgue 零集合N0⊆[t0,T]N_0\subseteq[t_0,T]が存在して、任意のt∈[t0,T]∖N0t\in[t_0,T]\setminus N_0と、(t,u),(t,v)∈S(t,u),(t,v)\in Sを満たす任意のu,vu,vについて∥f(t,u)−f(t,v)∥≤ℓ(t)∥u−v∥\|f(t,u)-f(t,v)\|\le\ell(t)\|u-v\|が成り立つとする。y,y~ ⁣:[t0,T]→Rdy,\tilde y\colon[t_0,T]\to\R^dを絶対連続写像とし、それらのグラフはSSに含まれ、ほとんど至る所でy′(t)=f(t,y(t))y'(t)=f(t,y(t))であるとする。y~\tilde yが微分可能な点でδ(t):=y~′(t)−f(t,y~(t))\delta(t):=\tilde y'(t)-f(t,\tilde y(t))と置き、δ\deltaが[t0,T][t_0,T]上で可積分であるとする。このときt∈[t0,T]t\in[t_0,T]について

∥y~(t)−y(t)∥≤∥y~(t0)−y(t0)∥exp⁡(∫t0tℓ(s) ds)+∫t0texp⁡(∫stℓ(σ) dσ)∥δ(s)∥ ds\|\tilde y(t)-y(t)\|\le\|\tilde y(t_0)-y(t_0)\|\exp\Bigl(\int_{t_0}^t\ell(s)\,ds\Bigr)+\int_{t_0}^t\exp\Bigl(\int_s^t\ell(\sigma)\,d\sigma\Bigr)\|\delta(s)\|\,ds

が成り立つ。

証明.y~\tilde yが微分可能な点ではδˉ(t):=δ(t)\bar\delta(t):=\delta(t)、それ以外の点ではδˉ(t):=0\bar\delta(t):=0と置き、g(t,u):=f(t,u)+δˉ(t)g(t,u):=f(t,u)+\bar\delta(t)（(t,u)∈S(t,u)\in S）と置く。ほとんど至る所でy~′(t)=g(t,y~(t))\tilde y'(t)=g(t,\tilde y(t))であり、任意のttについて∥f(t,y~(t))−g(t,y~(t))∥=∥δˉ(t)∥\|f(t,\tilde y(t))-g(t,\tilde y(t))\|=\|\bar\delta(t)\|である。§E13.7 定理 2.1を、区間[t0,T][t_0,T]、基準の時刻t0t_0、右辺ggの解y~\tilde y、右辺ffの解yy、ffの Lipschitz 係数ℓ\ellと零集合N0N_0、右辺の差の上界r:=∥δˉ∥∈L1([t0,T])r:=\|\bar\delta\|\in L^1([t_0,T])として適用すると、主張の評価を得る。▨

例 4.2.d∈N≥1d\in\NNとし、Rd\R^dにノルム∥⋅∥\|\cdot\|を固定する。t0<Tt_0<T、f ⁣:[t0,T]×Rd→Rdf\colon[t_0,T]\times\R^d\to\R^dとし、L,K≥0L,K\ge0について、任意のt,s∈[t0,T]t,s\in[t_0,T]とu,v∈Rdu,v\in\R^dに対して∥f(t,u)−f(t,v)∥≤L∥u−v∥\|f(t,u)-f(t,v)\|\le L\|u-v\|と∥f(t,u)−f(s,u)∥≤K∣t−s∣\|f(t,u)-f(s,u)\|\le K|t-s|が成り立つとする。y ⁣:[t0,T]→Rdy\colon[t_0,T]\to\R^dをy′=f(t,y)y'=f(t,y)を満たすC1C^1級写像とする。[t0,T][t_0,T]の格子とy0∈Rdy_0\in\R^dから前進 Euler 法の近似値yny_nを定め、t∈[tn,tn+1]t\in[t_n,t_{n+1}]についてy~(t):=yn+(t−tn)f(tn,yn)\tilde y(t):=y_n+(t-t_n)f(t_n,y_n)と置く。t=tn+1t=t_{n+1}ではyn+hnf(tn,yn)=yn+1y_n+h_nf(t_n,y_n)=y_{n+1}であるからy~\tilde yは連続であり、区分的に一次であるから絶対連続である。t∈(tn,tn+1)t\in(t_n,t_{n+1})ではy~′(t)=f(tn,yn)\tilde y'(t)=f(t_n,y_n)であり、

∥δ(t)∥=∥f(tn,yn)−f(t,y~(t))∥≤K(t−tn)+L∥yn−y~(t)∥=(K+L∥f(tn,yn)∥)(t−tn)\|\delta(t)\|=\|f(t_n,y_n)-f(t,\tilde y(t))\|\le K(t-t_n)+L\|y_n-\tilde y(t)\|=\bigl(K+L\|f(t_n,y_n)\|\bigr)(t-t_n)

である。δ\deltaは各(tn,tn+1)(t_n,t_{n+1})上で連続かつ有界であるから可積分であり、∫tntn+1∥δ(s)∥ ds≤(K+L∥f(tn,yn)∥)hn2/2\int_{t_n}^{t_{n+1}}\|\delta(s)\|\,ds\le(K+L\|f(t_n,y_n)\|)h_n^2/2である。定理 4.1をS=[t0,T]×RdS=[t_0,T]\times\R^d、ℓ=L\ell=Lとして適用し、eL(tn−s)≤eL(tn−t0)e^{L(t_n-s)}\le e^{L(t_n-t_0)}を用いると

∥yn−y(tn)∥≤eL(tn−t0)(∥y0−y(t0)∥+∑j=0n−1(K+L∥f(tj,yj)∥)hj22)\|y_n-y(t_n)\|\le e^{L(t_n-t_0)}\Bigl(\|y_0-y(t_0)\|+\sum_{j=0}^{n-1}\bigl(K+L\|f(t_j,y_j)\|\bigr)\frac{h_j^2}2\Bigr)

を得る。右辺はLL、KK、格子と計算した近似値だけで定まる。t→tn+0t\to t_n+0で∥δ(t)∥→0\|\delta(t)\|\to0であるから、各格子点でのδ\deltaの右極限は00である。一方、d=1d=1、f(t,u)=uf(t,u)=uではt∈(tn,tn+1)t\in(t_n,t_{n+1})でδ(t)=yn−y~(t)=−(t−tn)yn\delta(t)=y_n-\tilde y(t)=-(t-t_n)y_nであり、yn≠0y_n\ne0ならば∫tntn+1∣δ(s)∣ ds=hn2∣yn∣/2>0\int_{t_n}^{t_{n+1}}|\delta(s)|\,ds=h_n^2|y_n|/2>0である。

5 絶対安定性

定義 5.1.λ∈C\lambda\in\Cとする。C\Cを実22次元の線形空間R2\R^2と同一視し、ノルムとして絶対値∣⋅∣|\cdot|を用いる。f(t,u):=λuf(t,u):=\lambda u（(t,u)∈R×C(t,u)\in\R\times\C）による初期値問題y′=λyy'=\lambda y、y(0)=1y(0)=1を テスト方程式 (test equation) という。その解はy(t)=eλty(t)=e^{\lambda t}（t∈Rt\in\R）である。

定義 5.2. 一段法と、集合DR⊆CD_R\subseteq\Cと写像R ⁣:DR→CR\colon D_R\to\Cが次を満たすとき、RRを一段法の 安定関数 (stability function) という。任意のλ∈C\lambda\in\Cとh>0h>0について、テスト方程式の右辺f(t,u)=λuf(t,u)=\lambda uに対する増分関数の定義域をDD、一歩写像をΨh\Psi_hとすると、hλ∈DRh\lambda\in D_Rならば任意の(t,u)∈R×C(t,u)\in\R\times\Cについて(t,u,h)∈D(t,u,h)\in DかつΨh(t,u)=R(hλ)u\Psi_h(t,u)=R(h\lambda)uであり、hλ∉DRh\lambda\notin D_Rならば(t,u,h)∉D(t,u,h)\notin Dを満たす(t,u)(t,u)が存在する。DRD_RとRRは一段法から一意に定まる。実際、DRD_Rはh=1h=1、λ=z\lambda=zの場合の条件で定まり、R(z)=Ψ1(0,1)R(z)=\Psi_1(0,1)である。

S:={z∈DR∣∣R(z)∣≤1}S:=\{z\in D_R\mid|R(z)|\le1\}

を一段法の 絶対安定領域 (region of absolute stability) という。

命題 5.3. 一段法が安定関数R ⁣:DR→CR\colon D_R\to\Cと絶対安定領域SSをもつとする。λ∈C\lambda\in\C、h>0h>0、z:=hλ∈DRz:=h\lambda\in D_Rとし、テスト方程式の等間隔の格子tn=nht_n=nh上で、w0∈Cw_0\in\Cとr0,r1,…∈Cr_0,r_1,\ldots\in\Cからwn+1:=Ψh(tn,wn)+rnw_{n+1}:=\Psi_h(t_n,w_n)+r_nと定める。ε≥0\varepsilon\ge0が任意のnnについて∣rn∣≤ε|r_n|\le\varepsilonを満たすとする。

  1. z∈Sz\in Sであることとsup⁡n∣R(z)n∣<∞\sup_n|R(z)^n|<\inftyは同値であり、∣R(z)∣<1|R(z)|<1であることとR(z)n→0R(z)^n\to0（n→∞n\to\infty）は同値である。
  2. z∈Sz\in Sならば∣wn−R(z)nw0∣≤nε|w_n-R(z)^nw_0|\le n\varepsilonである。R(z)=1R(z)=1かつ任意のnnについてrn=εr_n=\varepsilonならば等号が成り立つ。
  3. ∣R(z)∣<1|R(z)|<1ならば∣wn−R(z)nw0∣≤ε/(1−∣R(z)∣)|w_n-R(z)^nw_0|\le\varepsilon/(1-|R(z)|)である。
  4. ∣R(z)∣>1|R(z)|>1、r0=εr_0=\varepsilon、rn=0r_n=0（n≥1n\ge1）ならば∣wn−R(z)nw0∣=∣R(z)∣n−1ε|w_n-R(z)^nw_0|=|R(z)|^{n-1}\varepsilon（n≥1n\ge1）である。

証明.R:=R(z)R:=R(z)と書く。Ψh(tn,w)=Rw\Psi_h(t_n,w)=Rwであるから、nnについての帰納法によりwn−Rnw0=∑j=0n−1Rn−1−jrjw_n-R^nw_0=\sum_{j=0}^{n-1}R^{n-1-j}r_jである。主張 (1)は∣Rn∣=∣R∣n|R^n|=|R|^nによる。主張 (2)では∣R∣n−1−j≤1|R|^{n-1-j}\le1により和の絶対値はnεn\varepsilon以下であり、R=1R=1、rj=εr_j=\varepsilonならば和はnεn\varepsilonである。主張 (3)では和の絶対値はε∑k≥0∣R∣k=ε/(1−∣R∣)\varepsilon\sum_{k\ge0}|R|^k=\varepsilon/(1-|R|)以下である。主張 (4)では和はRn−1εR^{n-1}\varepsilonである。▨

定義 5.4. 安定関数R ⁣:DR→CR\colon D_R\to\Cと絶対安定領域SSをもつ一段法が A 安定 (A-stable) であるとは、{z∈C∣ℜz≤0}⊆S\{z\in\C\mid\Re z\le0\}\subseteq Sであることをいう。

命題 5.5.θ∈[0,1]\theta\in[0,1]とする。θ 法は安定関数

Rθ(z)=1+(1−θ)z1−θzR_\theta(z)=\frac{1+(1-\theta)z}{1-\theta z}

をもち、その定義域はθ>0\theta>0ならばC∖{1/θ}\C\setminus\{1/\theta\}、θ=0\theta=0ならばC\Cである。

  1. z∈DRz\in D_Rについて、∣Rθ(z)∣≤1|R_\theta(z)|\le1であることと2ℜz+(1−2θ)∣z∣2≤02\Re z+(1-2\theta)|z|^2\le0は同値である。
  2. θ 法が A 安定であることとθ≥12\theta\ge\frac12は同値である。θ≥12\theta\ge\frac12かつℜz<0\Re z<0ならば∣Rθ(z)∣<1|R_\theta(z)|<1であり、θ=12\theta=\frac12かつℜz=0\Re z=0ならば∣R1/2(z)∣=1|R_{1/2}(z)|=1である。
  3. θ>0\theta>0ならば、実数x→−∞x\to-\inftyのときRθ(x)→−(1−θ)/θR_\theta(x)\to-(1-\theta)/\thetaである。特にR1(x)→0R_1(x)\to0、R1/2(x)→−1R_{1/2}(x)\to-1である。

証明.λ∈C\lambda\in\C、h>0h>0、z:=hλz:=h\lambdaとする。テスト方程式ではΩ=R×C\Omega=\R\times\Cであり、θ 法の方程式は(1−θz)v=(1+(1−θ)z)u(1-\theta z)v=(1+(1-\theta)z)uである。θz≠1\theta z\ne1ならば、任意の(t,u)(t,u)についてただ一つの解v=Rθ(z)uv=R_\theta(z)uが存在する。θ>0\theta>0かつz=1/θz=1/\thetaならば、方程式は0=u/θ0=u/\thetaとなり、u=1u=1では解が存在しない。したがってRθR_\thetaは定義域C∖{1/θ}\C\setminus\{1/\theta\}（θ=0\theta=0ではC\C）の安定関数である。

主張 (1)を示す。z∈DRz\in D_Rでは∣1−θz∣>0|1-\theta z|>0であり、

∣1+(1−θ)z∣2−∣1−θz∣2=2(1−θ)ℜz+(1−θ)2∣z∣2+2θℜz−θ2∣z∣2=2ℜz+(1−2θ)∣z∣2|1+(1-\theta)z|^2-|1-\theta z|^2=2(1-\theta)\Re z+(1-\theta)^2|z|^2+2\theta\Re z-\theta^2|z|^2=2\Re z+(1-2\theta)|z|^2

である。

主張 (2)を示す。θ≥12\theta\ge\frac12とし、ℜz≤0\Re z\le0とする。1/θ1/\thetaは正の実数であるからz∈DRz\in D_Rであり、2ℜz+(1−2θ)∣z∣2≤2ℜz≤02\Re z+(1-2\theta)|z|^2\le2\Re z\le0であるから主張 (1)によりz∈Sz\in Sである。ℜz<0\Re z<0ならば左辺は負であり∣Rθ(z)∣<1|R_\theta(z)|<1である。θ=12\theta=\frac12では左辺は2ℜz2\Re zであり、ℜz=0\Re z=0ならば∣R1/2(z)∣=1|R_{1/2}(z)|=1である。θ<12\theta<\frac12とし、x>2/(1−2θ)x>2/(1-2\theta)とすると、z=−xz=-xについて2ℜz+(1−2θ)∣z∣2=x((1−2θ)x−2)>02\Re z+(1-2\theta)|z|^2=x((1-2\theta)x-2)>0であるから−x∉S-x\notin Sであり、θ 法は A 安定でない。

主張 (3)は、x<0x<0についてRθ(x)=(x−1+1−θ)/(x−1−θ)R_\theta(x)=(x^{-1}+1-\theta)/(x^{-1}-\theta)であることによる。▨

系 5.6. 前進 Euler 法の安定関数はR(z)=1+zR(z)=1+z（DR=CD_R=\C）であり、絶対安定領域はS={z∈C∣∣1+z∣≤1}S=\{z\in\C\mid|1+z|\le1\}、S∩R=[−2,0]S\cap\R=[-2,0]である。前進 Euler 法は A 安定でない。後退 Euler 法の安定関数はR(z)=1/(1−z)R(z)=1/(1-z)（DR=C∖{1}D_R=\C\setminus\{1\}）であり、絶対安定領域はS={z∈C∣∣1−z∣≥1}S=\{z\in\C\mid|1-z|\ge1\}である。後退 Euler 法は A 安定である。

証明. テスト方程式ではΩ=R×C\Omega=\R\times\Cであり、前進 Euler 法の増分関数の定義域はΩ×(0,∞)\Omega\times(0,\infty)、一歩写像はu+hλu=(1+hλ)uu+h\lambda u=(1+h\lambda)uである。したがってR(z)=1+zR(z)=1+z、DR=CD_R=\Cである。実数xxについて∣1+x∣≤1|1+x|\le1と−2≤x≤0-2\le x\le0は同値であり、−3∉S-3\notin Sであるから前進 Euler 法は A 安定でない。

後退 Euler 法について、z:=hλ≠1z:=h\lambda\ne1ならば方程式v=u+zvv=u+zvはただ一つの解u/(1−z)u/(1-z)をもち、§E24.27 定義 3.1 (1)によりそれが一歩写像の値である。z=1z=1としu=1u=1とすると方程式v=1+vv=1+vは解をもたない。§E24.27 定義 3.1 (2)の値も方程式の解であるから、(t,1,h)∉D(t,1,h)\notin Dであり、安定関数はR(z)=1/(1−z)R(z)=1/(1-z)、DR=C∖{1}D_R=\C\setminus\{1\}である。z=1z=1は∣1−z∣≥1|1-z|\ge1を満たさないからS={z∈C∣∣1−z∣≥1}S=\{z\in\C\mid|1-z|\ge1\}である。この安定関数はR1R_1に等しいから、命題 5.5 (2)により後退 Euler 法は A 安定である。▨

命題 5.7.s∈N≥1s\in\NNとし、A∈Rs×sA\in\R^{s\times s}を狭義下三角行列、b∈Rsb\in\R^sとする。Butcher 配列(A,b)(A,b)の陽的 Runge–Kutta 法は、定義域C\Cの安定関数

R(z)=1+∑j=1s(bTAj−11)zjR(z)=1+\sum_{j=1}^s\bigl(b^{\mathsf T}A^{j-1}\mathbf1\bigr)z^j

をもつ。bT1≠0b^{\mathsf T}\mathbf1\ne0ならば絶対安定領域は有界であり、この方法は A 安定でない。古典的 Runge–Kutta 法の安定関数はR(z)=1+z+z22+z36+z424R(z)=1+z+\frac{z^2}2+\frac{z^3}6+\frac{z^4}{24}である。

証明.λ∈C\lambda\in\C、h>0h>0、z:=hλz:=h\lambda、(t,u)∈R×C(t,u)\in\R\times\Cとする。テスト方程式ではΩ=R×C\Omega=\R\times\Cであるから、§E24.27 定義 2.2 (3)の段の値はすべて定まり、k:=(k1,…,ks)T∈Csk:=(k_1,\dots,k_s)^{\mathsf T}\in\C^sはk=λu1+zAkk=\lambda u\mathbf1+zAkを満たす。AAは狭義下三角行列であるからAs=0A^s=0であり、(I−zA)−1=∑j=0s−1zjAj(I-zA)^{-1}=\sum_{j=0}^{s-1}z^jA^jである。したがってk=λu∑j=0s−1zjAj1k=\lambda u\sum_{j=0}^{s-1}z^jA^j\mathbf1であり、

Ψh(t,u)=u+hbTk=(1+∑j=0s−1(bTAj1)zj+1)u=R(z)u\Psi_h(t,u)=u+hb^{\mathsf T}k=\Bigl(1+\sum_{j=0}^{s-1}\bigl(b^{\mathsf T}A^j\mathbf1\bigr)z^{j+1}\Bigr)u=R(z)u

を得る。bT1≠0b^{\mathsf T}\mathbf1\ne0ならばRRは次数q≥1q\ge1の多項式であり、最高次の係数をβq≠0\beta_q\ne0、他の係数をβ0,…,βq−1\beta_0,\ldots,\beta_{q-1}とすると∣R(z)∣≥∣βq∣∣z∣q−∑k<q∣βk∣∣z∣k→∞|R(z)|\ge|\beta_q||z|^q-\sum_{k<q}|\beta_k||z|^k\to\infty（∣z∣→∞|z|\to\infty）である。したがってSSは有界であり、閉左半平面を含まない。古典的 Runge–Kutta 法ではc=A1=(0,12,12,1)Tc=A\mathbf1=(0,\frac12,\frac12,1)^{\mathsf T}、Ac=(0,0,14,12)TAc=(0,0,\frac14,\frac12)^{\mathsf T}、A2c=(0,0,0,14)TA^2c=(0,0,0,\frac14)^{\mathsf T}、b=16(1,2,2,1)Tb=\frac16(1,2,2,1)^{\mathsf T}であるから、bT1=1b^{\mathsf T}\mathbf1=1、bTc=12b^{\mathsf T}c=\frac12、bTAc=16b^{\mathsf T}Ac=\frac16、bTA2c=124b^{\mathsf T}A^2c=\frac1{24}である。▨

例 5.8. 台形法、後退 Euler 法の安定関数R1/2R_{1/2}、R1R_1とeze^zの値は次のとおりである。小数は小数第6位に丸めた。

zz R1/2(z)R_{1/2}(z) R1(z)R_1(z) eze^z
−10-10 −2/3-2/3 1/11=0.0909091/11=0.090909 0.0000450.000045
−2-2 00 1/31/3 0.1353350.135335
5i5i (−21+20i)/29(-21+20i)/29 (1+5i)/26(1+5i)/26 0.283662−0.958924i0.283662-0.958924i

z=5iz=5iでは∣R1/2(z)∣=1=∣ez∣|R_{1/2}(z)|=1=|e^z|、∣R1(z)∣=1/26=0.196116|R_1(z)|=1/\sqrt{26}=0.196116である。R1/2(5i)R_{1/2}(5i)の偏角は2arctan⁡52=2.3805802\arctan\frac52=2.380580であり、e5ie^{5i}の偏角5−2π=−1.2831855-2\pi=-1.283185と異なる。

命題 5.9.m∈N≥1m\in\NNとし、B∈Cm×mB\in\C^{m\times m}と正則行列V∈Cm×mV\in\C^{m\times m}がV−1BV=diag⁡(λ1,…,λm)V^{-1}BV=\operatorname{diag}(\lambda_1,\dots,\lambda_m)を満たすとする。一段法を θ 法（θ∈[0,1]\theta\in[0,1]）または陽的 Runge–Kutta 法とし、その安定関数をR ⁣:DR→CR\colon D_R\to\Cとする。Cm\C^mをR2m\R^{2m}と同一視し、f(t,u):=Buf(t,u):=Bu（(t,u)∈R×Cm(t,u)\in\R\times\C^m）とする。h>0h>0がhλ1,…,hλm∈DRh\lambda_1,\dots,h\lambda_m\in D_Rを満たすとする。

  1. 任意の(t,u)(t,u)について(t,u,h)∈D(t,u,h)\in Dであり、Ψh(t,u)=Vdiag⁡(R(hλ1),…,R(hλm))V−1u\Psi_h(t,u)=V\operatorname{diag}(R(h\lambda_1),\dots,R(h\lambda_m))V^{-1}uである。
  2. 1≤q≤∞1\le q\le\inftyとし、Cm\C^mにqqノルム∥⋅∥q\|\cdot\|_qと、それによる作用素ノルムを入れる。u0∈Cmu_0\in\C^mからun+1:=Ψh(tn,un)u_{n+1}:=\Psi_h(t_n,u_n)と定めると ∥un∥q≤∥V∥q∥V−1∥qmax⁡i∣R(hλi)∣n∥u0∥q\|u_n\|_q\le\|V\|_q\|V^{-1}\|_q\max_i|R(h\lambda_i)|^n\|u_0\|_q である。q=2q=2でVVがユニタリ行列ならば∥V∥2∥V−1∥2=1\|V\|_2\|V^{-1}\|_2=1である。

証明.主張 (1)を示す。x:=V−1ux:=V^{-1}uと置く。θ 法では、v∈Cmv\in\C^mが方程式v=u+h((1−θ)Bu+θBv)v=u+h((1-\theta)Bu+\theta Bv)を満たすことと、w:=V−1vw:=V^{-1}vの各成分がwl=xl+h((1−θ)λlxl+θλlwl)w_l=x_l+h((1-\theta)\lambda_lx_l+\theta\lambda_lw_l)を満たすことは同値である。各成分の方程式はλ=λl\lambda=\lambda_lのテスト方程式に対する θ 法の方程式であり、hλl∈DRh\lambda_l\in D_Rであるから、ただ一つの解wl=R(hλl)xlw_l=R(h\lambda_l)x_lをもつ。陽的 Runge–Kutta 法では、段の値ki∈Cmk_i\in\C^mに対してk~i:=V−1ki\tilde k_i:=V^{-1}k_iと置くとk~i=diag⁡(λ1,…,λm)(x+h∑j<iaijk~j)\tilde k_i=\operatorname{diag}(\lambda_1,\ldots,\lambda_m)(x+h\sum_{j<i}a_{ij}\tilde k_j)であり、各成分はλ=λl\lambda=\lambda_lのテスト方程式に対する段の値である。したがってV−1Ψh(t,u)V^{-1}\Psi_h(t,u)の第ll成分はR(hλl)xlR(h\lambda_l)x_lである。

主張 (2)を示す。μ∈Cm\mu\in\C^mについて∥diag⁡(μ)y∥q≤max⁡i∣μi∣∥y∥q\|\operatorname{diag}(\mu)y\|_q\le\max_i|\mu_i|\|y\|_qである。主張 (1)によりun=Vdiag⁡(R(hλi)n)V−1u0u_n=V\operatorname{diag}(R(h\lambda_i)^n)V^{-1}u_0であり、作用素ノルムの劣乗法性から評価を得る。ユニタリ行列VVは∥Vy∥2=∥y∥2\|Vy\|_2=\|y\|_2を満たすから∥V∥2=∥V−1∥2=1\|V\|_2=\|V^{-1}\|_2=1である。▨

例 5.10.

B=(−11000−2),V=(1−10001),V−1=(110001)B=\begin{pmatrix}-1&100\\0&-2\end{pmatrix},\qquad V=\begin{pmatrix}1&-100\\0&1\end{pmatrix},\qquad V^{-1}=\begin{pmatrix}1&100\\0&1\end{pmatrix}

とするとV−1BV=diag⁡(−1,−2)V^{-1}BV=\operatorname{diag}(-1,-2)である。前進 Euler 法をh=12h=\frac12で用いるとhλ1=−12h\lambda_1=-\frac12、hλ2=−1h\lambda_2=-1はともにS∩R=[−2,0]S\cap\R=[-2,0]に属し、∣R(hλ1)∣=12|R(h\lambda_1)|=\frac12、R(hλ2)=0R(h\lambda_2)=0である。一方、Ψh(t,u)=(I+12B)u\Psi_h(t,u)=(I+\frac12B)u、I+12B=(1/25000)I+\frac12B=\begin{pmatrix}1/2&50\\0&0\end{pmatrix}であり、u0=(0,1)Tu_0=(0,1)^{\mathsf T}からu1=(50,0)Tu_1=(50,0)^{\mathsf T}、un=(50⋅21−n,0)Tu_n=(50\cdot2^{1-n},0)^{\mathsf T}（n≥1n\ge1）を得る。したがって∥u1∥∞=50∥u0∥∞\|u_1\|_\infty=50\|u_0\|_\inftyである。命題 5.9 (2)の係数は∥V∥∞∥V−1∥∞=1012=10201\|V\|_\infty\|V^{-1}\|_\infty=101^2=10201である。

6 硬い問題

定義 6.1.m∈N≥1m\in\NNとし、B∈Cm×mB\in\C^{m\times m}を対角化可能な行列、λ1,…,λm\lambda_1,\dots,\lambda_mを重複を込めた固有値とし、任意のiiについてℜλi<0\Re\lambda_i<0であるとする。

max⁡i∣ℜλi∣min⁡i∣ℜλi∣\frac{\max_i|\Re\lambda_i|}{\min_i|\Re\lambda_i|}

をBBの 硬さ比 (stiffness ratio) という。安定関数と絶対安定領域SSをもつ一段法に対して、

hS:=sup⁡{h>0 ∣ 0<h′≤h を満たす任意の h′ と任意の i について h′λi∈S}(sup⁡∅:=0)h_S:=\sup\bigl\{h>0\ \big|\ 0<h'\le h\ \text{を満たす任意の}\ h'\ \text{と任意の}\ i\ \text{について}\ h'\lambda_i\in S\bigr\}\qquad(\sup\varnothing:=0)

と置く。初期値問題u′=Buu'=Buが区間[t0,T][t_0,T]とこの一段法について 硬い (stiff) であるとは、hSh_Sが、局所打切り誤差の大きさから[t0,T][t_0,T]上の精度の要求が許す刻みより著しく小さいことをいう。この語は二つの刻みの大小を比べて述べる語であり、閾値を定めない。

例 6.2.

  1. λ=−100\lambda=-100のテスト方程式に前進 Euler 法を用いる。系 5.6により、hλ∈Sh\lambda\in Sであることと0<h≤0.020<h\le0.02は同値である。h=0.019h=0.019ではR(hλ)=−0.9R(h\lambda)=-0.9、h=0.021h=0.021ではR(hλ)=−1.1R(h\lambda)=-1.1であり、近似値はyn=R(hλ)ny_n=R(h\lambda)^nである。後退 Euler 法をh=0.1h=0.1で用いるとyn=11−ny_n=11^{-n}である。y(tn)=e−100tny(t_n)=e^{-100t_n}と比べると次のとおりである。指数表記の値は有効数字7桁に丸めた。

    nn 前進 Euler 法h=0.019h=0.019 e−1.9ne^{-1.9n} 前進 Euler 法h=0.021h=0.021 e−2.1ne^{-2.1n} 後退 Euler 法h=0.1h=0.1 e−10ne^{-10n}
    11 −0.9-0.9 1.495686×10−11.495686\times10^{-1} −1.1-1.1 1.224564×10−11.224564\times10^{-1} 9.090909×10−29.090909\times10^{-2} 4.539993×10−54.539993\times10^{-5}
    22 0.810.81 2.237077×10−22.237077\times10^{-2} 1.211.21 1.499558×10−21.499558\times10^{-2} 8.264463×10−38.264463\times10^{-3} 2.061154×10−92.061154\times10^{-9}
    33 −0.729-0.729 3.345965×10−33.345965\times10^{-3} −1.331-1.331 1.836305×10−31.836305\times10^{-3} 7.513148×10−47.513148\times10^{-4} 9.357623×10−149.357623\times10^{-14}
    44 0.65610.6561 5.004514×10−45.004514\times10^{-4} 1.46411.4641 2.248673×10−42.248673\times10^{-4} 6.830135×10−56.830135\times10^{-5} 4.248354×10−184.248354\times10^{-18}
    55 −0.59049-0.59049 7.485183×10−57.485183\times10^{-5} −1.61051-1.61051 2.753645×10−52.753645\times10^{-5} 6.209213×10−66.209213\times10^{-6} 1.928750×10−221.928750\times10^{-22}
    66 0.5314410.531441 1.119548×10−51.119548\times10^{-5} 1.7715611.771561 3.372015×10−63.372015\times10^{-6} 5.644739×10−75.644739\times10^{-7} 8.756511×10−278.756511\times10^{-27}

    同じ最終時刻T=0.399=21×0.019=19×0.021T=0.399=21\times0.019=19\times0.021では、h=0.019h=0.019の近似値は(−0.9)21=−0.109418989131512359209(-0.9)^{21}=-0.109418989131512359209、h=0.021h=0.021の近似値は(−1.1)19=−6.115909…(-1.1)^{19}=-6.115909\ldotsであり、y(T)=e−39.9=4.695158×10−18y(T)=e^{-39.9}=4.695158\times10^{-18}である。

  2. B=diag⁡(−1,−1000)B=\operatorname{diag}(-1,-1000)の硬さ比は10001000である。前進 Euler 法では、∣1−h∣≤1|1-h|\le1と∣1−1000h∣≤1|1-1000h|\le1がともに成り立つことと0<h≤0.0020<h\le0.002は同値であり、hS=0.002h_S=0.002である。h=0.002h=0.002、u0=(1,1)Tu_0=(1,1)^{\mathsf T}ではun=(0.998n,(−1)n)Tu_n=(0.998^n,(-1)^n)^{\mathsf T}であり、第2成分の絶対値は任意のnnで11であるが、厳密解の第2成分はe−1000tn=e−2ne^{-1000t_n}=e^{-2n}である。後退 Euler 法をh=0.1h=0.1で用いると、二つの成分の増幅係数はR1(−0.1)=1/1.1R_1(-0.1)=1/1.1とR1(−100)=1/101R_1(-100)=1/101である。

  3. B=diag⁡(−1,−1000)B=\operatorname{diag}(-1,-1000)、区間[0,1][0,1]、初期値u(0)=(1,10−8)Tu(0)=(1,10^{-8})^{\mathsf T}とし、R2\R^2に最大値ノルム∥⋅∥∞\|\cdot\|_\inftyを入れて前進 Euler 法をh=0.01h=0.01で用いる。解はu(t)=(e−t,10−8e−1000t)Tu(t)=(e^{-t},10^{-8}e^{-1000t})^{\mathsf T}、一歩写像はΨh(t,v)=(I+hB)v\Psi_h(t,v)=(I+hB)vであるから、0≤t≤1−h0\le t\le1-hについて局所打切り誤差は

    δu(t,h)=(e−t(e−h−1+h), 10−8e−1000t(e−1000h−1+1000h))T\delta_u(t,h)=\bigl(e^{-t}(e^{-h}-1+h),\ 10^{-8}e^{-1000t}(e^{-1000h}-1+1000h)\bigr)^{\mathsf T}

    である。§E24.27 補題 1.5 (1)をτ↦e−τ\tau\mapsto e^{-\tau}とr=1r=1に適用すると∣e−0.01−1+0.01∣≤0.012/2=5×10−5|e^{-0.01}-1+0.01|\le0.01^2/2=5\times10^{-5}であり、0<e−10<10<e^{-10}<1から0<e−10−1+10<100<e^{-10}-1+10<10である。したがって0≤t≤0.990\le t\le0.99について∥δu(t,0.01)∥∞≤max⁡{5×10−5,10−7}<10−4\|\delta_u(t,0.01)\|_\infty\le\max\{5\times10^{-5},10^{-7}\}<10^{-4}であり、各一歩の局所打切り誤差の最大値ノルムを10−410^{-4}未満とする要求を刻み0.010.01は満たす。一方、(2)により同じ前進 Euler 法のhSh_Sは0.0020.002であり、刻み0.010.01はその55倍である。h=0.01h=0.01では第2成分の増幅係数はR(−10)=−9R(-10)=-9であり、u0=u(0)u_0=u(0)から定まる近似値はun=(0.99n, 10−8(−9)n)Tu_n=(0.99^n,\ 10^{-8}(-9)^n)^{\mathsf T}である。∥u9∥∞=99×10−8>3\|u_9\|_\infty=9^9\times10^{-8}>3であるが、0≤t≤10\le t\le1について∥u(t)∥∞≤1\|u(t)\|_\infty\le1である。

  4. m=1m=1、B=(−1000)B=(-1000)の硬さ比は11である。前進 Euler 法ではhS=0.002h_S=0.002である。区間[0,10][0,10]で後退 Euler 法をh=1h=1で用いるとyn=1001−ny_n=1001^{-n}であり、yny_nとy(tn)=e−1000ny(t_n)=e^{-1000n}はともに正であるから、n≥1n\ge1について∣yn−y(tn)∣<1001−n<10−3|y_n-y(t_n)|<1001^{-n}<10^{-3}である。精度の要求を格子点での絶対誤差10−310^{-3}以下とすると、後退 Euler 法は刻み11でこの要求を満たし、前進 Euler 法のhS=0.002h_S=0.002はその1/5001/500である。

7 演習

問題 7.1.B=diag⁡(−1,−1000)B=\operatorname{diag}(-1,-1000)、u(0)=(1,1)Tu(0)=(1,1)^{\mathsf T}とし、台形法と後退 Euler 法をh=0.1h=0.1で1010歩用いる。それぞれのu10u_{10}の各成分と、u(1)=(e−1,e−1000)Tu(1)=(e^{-1},e^{-1000})^{\mathsf T}との差を求めよ。

解答.

命題 5.9により、各成分は安定関数をhλlh\lambda_lで評価した値の冪である。台形法ではR1/2(−0.1)=0.95/1.05=19/21R_{1/2}(-0.1)=0.95/1.05=19/21、R1/2(−100)=(1−50)/(1+50)=−49/51R_{1/2}(-100)=(1-50)/(1+50)=-49/51であり、

u10=((19/21)10, (49/51)10)T=(0.367572542…, 0.670284288…)Tu_{10}=\bigl((19/21)^{10},\ (49/51)^{10}\bigr)^{\mathsf T}=(0.367572542\ldots,\ 0.670284288\ldots)^{\mathsf T}

である。後退 Euler 法ではR1(−0.1)=10/11R_1(-0.1)=10/11、R1(−100)=1/101R_1(-100)=1/101であり、

u10=((10/11)10, 101−10)T=(0.385543289…, 9.0528695469298328727…×10−21)Tu_{10}=\bigl((10/11)^{10},\ 101^{-10}\bigr)^{\mathsf T}=(0.385543289\ldots,\ 9.0528695469298328727\ldots\times10^{-21})^{\mathsf T}

である。e−1=0.367879441…e^{-1}=0.367879441\ldotsであるから、第1成分の差の絶対値は台形法で0.000306898788573…0.000306898788573\ldots、後退 Euler 法で0.017663848…0.017663848\ldotsである。第2成分の差は台形法で(49/51)10−e−1000(49/51)^{10}-e^{-1000}、後退 Euler 法で101−10−e−1000101^{-10}-e^{-1000}である。0<e−1000<10−4340<e^{-1000}<10^{-434}であるから、これらはそれぞれ(49/51)10(49/51)^{10}と101−10101^{-10}との差が10−43410^{-434}未満の正の数である。▨

問題 7.2. 中点法と Heun 法の安定関数がR(z)=1+z+z22R(z)=1+z+\frac{z^2}2であり、S∩R=[−2,0]S\cap\R=[-2,0]であることを示せ。λ=−100\lambda=-100のテスト方程式でh=0.019h=0.019とh=0.021h=0.021に対するR(hλ)R(h\lambda)を求めよ。

解答.

中点法ではA=(00120)A=\begin{pmatrix}0&0\\\frac12&0\end{pmatrix}、b=(0,1)Tb=(0,1)^{\mathsf T}、Heun 法ではA=(0010)A=\begin{pmatrix}0&0\\1&0\end{pmatrix}、b=(12,12)Tb=(\frac12,\frac12)^{\mathsf T}であり、どちらもbT1=1b^{\mathsf T}\mathbf1=1、bTA1=12b^{\mathsf T}A\mathbf1=\frac12である。命題 5.7によりR(z)=1+z+z22R(z)=1+z+\frac{z^2}2である。実数xxについてR(x)+1=(x+1)2+32>0R(x)+1=\frac{(x+1)^2+3}2>0であるから、∣R(x)∣≤1|R(x)|\le1はR(x)−1=x(1+x2)≤0R(x)-1=x(1+\frac x2)\le0と同値であり、これは−2≤x≤0-2\le x\le0と同値である。h=0.019h=0.019ではR(−1.9)=1−1.9+1.805=0.905R(-1.9)=1-1.9+1.805=0.905、h=0.021h=0.021ではR(−2.1)=1−2.1+2.205=1.105R(-2.1)=1-2.1+2.205=1.105である。▨

前提記事