§E20.32二点境界値問題

最終更新

二階の微分方程式を一段法で扱うときは、解とその導関数を成分とする一階の系に直し、区間の一方の端での解の値と導関数の値の組を初期値として、解を他方の端へ進める。有限閉区間[a,b][a,b]上の実連続関数c≥0c\ge0、ffについて、方程式−u′′+cu=f-u''+cu=fに両端の値u(a)=αu(a)=\alpha、u(b)=βu(b)=\betaを課した問題は、C2C^2級の解をただ一つもつ。しかし左端で与えられるのは解の値だけであり、一段法で解を進めるのに必要な傾きu′(a)u'(a)は与えられていない。

傾きssを仮定して初期値問題を解き、右端での値がβ\betaに一致するようにssを定める方法を射撃法という。この問題では右端での値とβ\betaの差はssの一次関数であり、s=0s=0とs=1s=1の二つの初期値問題の解から正しい傾きが定まる。ところが[a,b]=[0,1][a,b]=[0,1]、c=402c=40^2、f=0f=0、α=1\alpha=1、β=0\beta=0の場合、解は00と11の間にとどまるのに対して、左端の値を10−1610^{-16}だけ変えた初期値問題の解は右端で1111以上変わる。

もう一つの方法は、幅hhの等間隔格子の内部の点で二階導関数を二階中心差分に置き換え、両端の値を直接課すことである。こうして作る中心差分方程式は、内部の格子点での値についての連立一次方程式と同値であり、その係数行列である差分行列は、既出の三重対角行列AnA_nのh−2h^{-2}倍と、ccの格子点での値を並べた対角行列との和である。差分行列の条件数は、ccと区間を固定して格子を細かくするとh−2h^{-2}の程度で増大する。それでも上と同じ係数で左端の値をϑ≥0\vartheta\ge0だけ変えたとき、離散解の変化は離散最大原理によってすべての格子点でϑ\vartheta以下である。

本記事は、この境界値問題に射撃法と中心差分法を構成し、解がC4C^4級であるとき、中心差分法の離散解と解の格子点での差が最大ノルムでh2h^2の定数倍以下であることを示す。

1 連続問題と射撃法

命題 1.1.a<ba<bとし、c,f ⁣:[a,b]→Rc,f\colon[a,b]\to\Rを連続関数、c≥0c\ge0、α,β∈R\alpha,\beta\in\Rとする。このとき

−u′′+cu=f,u(a)=α,u(b)=β-u''+cu=f,\qquad u(a)=\alpha,\qquad u(b)=\beta

を満たすu∈C2([a,b];R)u\in C^2([a,b];\R)がただ一つ存在する。

証明.w∈C2([a,b];R)w\in C^2([a,b];\R)がw′′−cw=0w''-cw=0、w(a)=w(b)=0w(a)=w(b)=0を満たすとする。部分積分により

0=∫ab(w′′−cw)w dx=[w′w]ab−∫ab(w′)2 dx−∫abcw2 dx=−∫ab(w′)2 dx−∫abcw2 dx0=\int_a^b(w''-cw)w\,dx=\bigl[w'w\bigr]_a^b-\int_a^b(w')^2\,dx-\int_a^bcw^2\,dx=-\int_a^b(w')^2\,dx-\int_a^bcw^2\,dx

であり、二つの被積分関数は連続かつ非負であるからw′=0w'=0である。wwは定数であり、w(a)=0w(a)=0からw=0w=0である。したがって、P=0P=0、Q=−cQ=-c、境界係数(αa,βa)=(αb,βb)=(1,0)(\alpha_a,\beta_a)=(\alpha_b,\beta_b)=(1,0)の斉次境界値問題は自明解だけをもつ。

ℓ(x):=α+(β−α)(x−a)/(b−a)\ell(x):=\alpha+(\beta-\alpha)(x-a)/(b-a)と置く。−f+cℓ-f+c\ellは連続であるから、§E10.13 系 3.6により

v′′−cv=−f+cℓ,v(a)=v(b)=0v''-cv=-f+c\ell,\qquad v(a)=v(b)=0

を満たすv∈C2([a,b];R)v\in C^2([a,b];\R)がただ一つ存在する。ℓ′′=0\ell''=0であるから、u:=v+ℓu:=v+\ellはu′′−cu=v′′−cv−cℓ=−fu''-cu=v''-cv-c\ell=-fを満たし、u(a)=αu(a)=\alpha、u(b)=βu(b)=\betaである。二つの解の差はw′′−cw=0w''-cw=0、w(a)=w(b)=0w(a)=w(b)=0を満たすC2C^2級関数であるから、第一段により00である。▨

命題 1.2.a<ba<b、L:=b−aL:=b-aとし、c,f ⁣:[a,b]→Rc,f\colon[a,b]\to\Rを連続関数、c≥0c\ge0、α,β∈R\alpha,\beta\in\Rとし、uuを命題 1.1の解とする。

  1. q′′=cqq''=cq、q(a)=0q(a)=0、q′(a)=1q'(a)=1を満たすq∈C2([a,b];R)q\in C^2([a,b];\R)がただ一つ存在する。qqは[a,b][a,b]上で狭義増加し、q(b)≥Lq(b)\ge Lである。
  2. 任意のs∈Rs\in\Rについて、y′′=cy−fy''=cy-f、y(a)=αy(a)=\alpha、y′(a)=sy'(a)=sを満たすy∈C2([a,b];R)y\in C^2([a,b];\R)がただ一つ存在し、それはys:=u+(s−u′(a))qy_s:=u+(s-u'(a))qである。
  3. s∈Rs\in\Rに対してR(s):=ys(b)−βR(s):=y_s(b)-\betaと置くと、任意のs∈Rs\in\RについてR(s)=R(0)+sq(b)R(s)=R(0)+sq(b)である。RRの零点はs∗:=u′(a)s_*:=u'(a)だけであり、s∗=−R(0)/(R(1)−R(0))s_*=-R(0)/(R(1)-R(0))、ys∗=uy_{s_*}=uである。
  4. s∈Rs\in\R、ε,η≥0\varepsilon,\eta\ge0とし、実数y^\hat yが∣y^−ys(b)∣≤ε|\hat y-y_s(b)|\le\varepsilonと∣y^−β∣≤η|\hat y-\beta|\le\etaを満たすとする。このとき ∣s−s∗∣≤η+εq(b)≤η+εL|s-s_*|\le\frac{\eta+\varepsilon}{q(b)}\le\frac{\eta+\varepsilon}{L} であり、任意のx∈[a,b]x\in[a,b]について∣ys(x)−u(x)∣≤η+ε|y_s(x)-u(x)|\le\eta+\varepsilonである。

証明.(1)を示す。§E10.13 補題 3.2と§E10.13 補題 3.3をP=0P=0、Q=−cQ=-c、(αa,βa)=(−1,0)(\alpha_a,\beta_a)=(-1,0)に適用すると、q′′=cqq''=cq、q(a)=0q(a)=0、q′(a)=1q'(a)=1を満たすqqはC1C^1級で導関数が絶対連続な関数の中でただ一つ存在し、C2C^2級である。q(a)=0q(a)=0、q′(a)=1q'(a)=1であるから、δ∈(0,L)\delta\in(0,L)が存在してqqは(a,a+δ](a,a+\delta]上で正である。Z:={x∈[a+δ,b]∣q(x)≤0}Z:=\{x\in[a+\delta,b]\mid q(x)\le0\}が空でないと仮定し、x1:=min⁡Zx_1:=\min Zと置く。(a,x1)(a,x_1)上でq′′=cq≥0q''=cq\ge0であるからq′q'は[a,x1][a,x_1]上で単調非減少であり、q′≥q′(a)=1q'\ge q'(a)=1である。平均値の定理によりq(x1)≥x1−a>0q(x_1)\ge x_1-a>0であり、これはq(x1)≤0q(x_1)\le0と両立しない。したがってqqは(a,b](a,b]上で正であり、[a,b][a,b]上でq′′=cq≥0q''=cq\ge0、q′≥1q'\ge1である。平均値の定理によりqqは狭義増加し、q(b)≥q(a)+L=Lq(b)\ge q(a)+L=Lである。

(2)を示す。u′′=cu−fu''=cu-fとq′′=cqq''=cqからys′′=cys−fy_s''=cy_s-fであり、ys(a)=αy_s(a)=\alpha、ys′(a)=u′(a)+s−u′(a)=sy_s'(a)=u'(a)+s-u'(a)=sである。yyが同じ条件を満たすとすると、w:=y−ysw:=y-y_sはw′′=cww''=cw、w(a)=w′(a)=0w(a)=w'(a)=0を満たし、q+wq+wは(1)の条件を満たす。(1)の一意性によりq+w=qq+w=qであり、y=ysy=y_sである。

(3)を示す。u(b)=βu(b)=\betaからR(s)=(s−u′(a))q(b)R(s)=(s-u'(a))q(b)であり、R(0)=−u′(a)q(b)R(0)=-u'(a)q(b)であるからR(s)=R(0)+sq(b)R(s)=R(0)+sq(b)である。(1)によりq(b)>0q(b)>0であるから、RRの零点はu′(a)u'(a)だけである。R(1)−R(0)=q(b)R(1)-R(0)=q(b)から−R(0)/(R(1)−R(0))=u′(a)-R(0)/(R(1)-R(0))=u'(a)であり、yu′(a)=uy_{u'(a)}=uである。

(4)を示す。∣R(s)∣≤∣ys(b)−y^∣+∣y^−β∣≤ε+η|R(s)|\le|y_s(b)-\hat y|+|\hat y-\beta|\le\varepsilon+\etaであり、(3)により∣R(s)∣=∣s−s∗∣q(b)|R(s)|=|s-s_*|q(b)である。q(b)≥Lq(b)\ge Lから第一の評価を得る。ys−u=(s−s∗)qy_s-u=(s-s_*)qであり、(1)により任意のx∈[a,b]x\in[a,b]について0≤q(x)≤q(b)0\le q(x)\le q(b)であるから、∣ys(x)−u(x)∣≤∣s−s∗∣q(b)≤η+ε|y_s(x)-u(x)|\le|s-s_*|q(b)\le\eta+\varepsilonである。▨

注意 1.3.命題 1.2 (2)の初期値問題は、(y,v)(a)=(α,s)(y,v)(a)=(\alpha,s)を初期値とする一階の系y′=vy'=v、v′=c(x)y−f(x)v'=c(x)y-f(x)と同値である。射撃法は、この系をs=0s=0とs=1s=1について一段法でaaからbbまで積分し、終端値の第一成分y^0\hat y_0、y^1\hat y_1がy^1≠y^0\hat y_1\ne\hat y_0を満たすときs^:=−(y^0−β)/(y^1−y^0)\hat s:=-(\hat y_0-\beta)/(\hat y_1-\hat y_0)を作る。s^\hat sの計算にuuとu′(a)u'(a)は現れない。s^\hat sについて再び積分した終端値の第一成分をy^\hat yとすると、∣y^−β∣|\hat y-\beta|は計算された値であり、これを命題 1.2 (4)のη\etaに取る。

ε\varepsilonは計算された値ではない。R2\R^2に最大ノルムを入れ、c,fc,fが[a,b][a,b]を含む開区間JJ上の連続関数の制限であり、(ys,ys′)(y_s,y_s')がJJ上で上の系を満たすC1C^1級写像に延長されるとする。延長した解について§E20.28 系 1.7の仮定がρ\rho、Λ\Lambda、H0H_0、pp、CCで成り立ち、[a,b][a,b]の格子の最大刻みHHがH≤H0H\le H_0を満たすとする。初期誤差E0E_0と各歩で加わる計算誤差rnr_nについてσ:=max⁡n∥rn∥/hn\sigma:=\max_n\|r_n\|/h_n、φ:=(eΛL−1)/Λ\varphi:=(e^{\Lambda L}-1)/\Lambda（Λ=0\Lambda=0ではφ:=L\varphi:=L）と置いてeΛLE0+φ(CHp+σ)≤ρe^{\Lambda L}E_0+\varphi(CH^p+\sigma)\le\rhoならば、この左辺をε\varepsilonに取ることができる。局所打切り誤差がChp+1Ch^{p+1}以下であることは同系の仮定であり、c,fc,fの連続性からは従わない。

刻み制御で積分する場合、§E20.30 命題 2.3により手続きはt=bt=bまたはh<hmin⁡h<h_{\min}で停止し、後者の停止では終端bbでの値は得られない。t=bt=bで停止し、受理された各刻みの局所打切り誤差の上界τn\tau_nが与えられたときは、§E20.30 系 2.4の評価をε\varepsilonに取ることができる。受理は誤差の推定値による判定であり、τn\tau_nを与えない。

2 中心差分法と差分行列

定義 2.1.a<ba<b、L:=b−aL:=b-aとし、c,f ⁣:[a,b]→Rc,f\colon[a,b]\to\Rを連続関数、α,β∈R\alpha,\beta\in\R、n∈N≥1n\in\NNとする。h:=L/(n+1)h:=L/(n+1)、xi:=a+ihx_i:=a+ih（0≤i≤n+10\le i\le n+1）、ci:=c(xi)c_i:=c(x_i)、fi:=f(xi)f_i:=f(x_i)と置き、集合{x0,…,xn+1}\{x_0,\dots,x_{n+1}\}上の実数値関数VVを値の列(V0,…,Vn+1)(V_0,\dots,V_{n+1})、Vi:=V(xi)V_i:=V(x_i)と同一視する。

U0=α,Un+1=β,−Dh2U(xi)+ciUi=fi(1≤i≤n)U_0=\alpha,\qquad U_{n+1}=\beta,\qquad -D^2_hU(x_i)+c_iU_i=f_i\quad(1\le i\le n)

を境界値問題−u′′+cu=f-u''+cu=f、u(a)=αu(a)=\alpha、u(b)=βu(b)=\betaの格子幅hhの 中心差分方程式 (central difference equation) といい、その解(U0,…,Un+1)(U_0,\dots,U_{n+1})を 離散解 (discrete solution) という。AnA_nを§E20.10 命題 5.1の行列とし、

Bh:=h−2An+diag⁡(c1,…,cn)B_h:=h^{-2}A_n+\operatorname{diag}(c_1,\dots,c_n)

を 差分行列 (difference matrix) という。e1,…,ene_1,\dots,e_nをRn\R^nの標準基底としてF:=(f1,…,fn)T+h−2(αe1+βen)F:=(f_1,\dots,f_n)^{\mathsf T}+h^{-2}(\alpha e_1+\beta e_n)と置く。n=1n=1ではF=(f1+(α+β)/h2)F=(f_1+(\alpha+\beta)/h^2)である。

命題 2.2.定義 2.1の設定でc≥0c\ge0とする。

  1. BhB_hは実対称かつ正定値であり、任意のv∈Rnv\in\R^nについてvTBhv=h−2vTAnv+∑i=1ncivi2≥h−2vTAnvv^{\mathsf T}B_hv=h^{-2}v^{\mathsf T}A_nv+\sum_{i=1}^nc_iv_i^2\ge h^{-2}v^{\mathsf T}A_nvである。
  2. v∈Rnv\in\R^nにv0:=vn+1:=0v_0:=v_{n+1}:=0を添えると、1≤i≤n1\le i\le nについて(Bhv)i=−Dh2v(xi)+civi(B_hv)_i=-D^2_hv(x_i)+c_iv_iである。
  3. 中心差分方程式の離散解はただ一つ存在し、(U1,…,Un)T(U_1,\dots,U_n)^{\mathsf T}はBhU^=FB_h\hat U=Fのただ一つの解U^\hat Uである。
  4. BhB_hのピボット選択なしの消去を厳密算術で実行すると、どの段でも失敗せず、すべてのピボットは正である。§E20.5 定理 6.1 (1)の漸化式による消去、前進代入と後退代入はBhU^=FB_h\hat U=Fの解を8n−78n-7回の四則演算で与える。

証明.(1)を示す。AnA_nと対角行列は対称であるからBhB_hは対称であり、等式はBhB_hの定義から従う。ci≥0c_i\ge0から不等式が成り立つ。§E20.10 命題 5.1 (1)によりv≠0v\ne0ならばvTAnv>0v^{\mathsf T}A_nv>0であるから、BhB_hは正定値である。

(2)を示す。AnA_nの第ii行から(Anv)i=2vi−vi−1−vi+1(A_nv)_i=2v_i-v_{i-1}-v_{i+1}であり、i=1i=1とi=ni=nで現れない項はv0=vn+1=0v_0=v_{n+1}=0により00である。−Dh2v(xi)=(2vi−vi−1−vi+1)/h2-D^2_hv(x_i)=(2v_i-v_{i-1}-v_{i+1})/h^2から等式を得る。

(3)を示す。(U0,…,Un+1)(U_0,\dots,U_{n+1})を実数の列とし、U^:=(U1,…,Un)T\hat U:=(U_1,\dots,U_n)^{\mathsf T}と置く。1≤i≤n1\le i\le nについて、−Dh2U(xi)+ciUi-D^2_hU(x_i)+c_iU_iは、U^\hat Uに00を添えた列に(2)を適用した値(BhU^)i(B_h\hat U)_iから、i=1i=1のときU0/h2U_0/h^2を、i=ni=nのときUn+1/h2U_{n+1}/h^2を引いた値である。したがってU0=αU_0=\alpha、Un+1=βU_{n+1}=\betaの下で、中心差分方程式はBhU^=FB_h\hat U=Fと同値である。(1)によりBhB_hは正則であるから、解はただ一つである。

(4)は、(1)と三重対角なBhB_hに§E20.5 定理 6.1 (3)と§E20.5 定理 6.1 (1)を適用して得られる。▨

命題 2.3.定義 2.1の設定でc≥0c\ge0とし、C:=max⁡[a,b]cC:=\max_{[a,b]}c、1≤k≤n1\le k\le nについてμk:=4h−2sin⁡2(kπ/(2(n+1)))\mu_k:=4h^{-2}\sin^2\bigl(k\pi/(2(n+1))\bigr)と置き、ν1≤⋯≤νn\nu_1\le\dots\le\nu_nをBhB_hの固有値を重複度を込めて並べたものとする。

  1. 1≤k≤n1\le k\le nについてμk≤νk≤μk+C\mu_k\le\nu_k\le\mu_k+Cである。
  2. 4/L2≤μ1≤π2/L24/L^2\le\mu_1\le\pi^2/L^2、2h−2≤μn≤4h−22h^{-2}\le\mu_n\le4h^{-2}である。
  3. 2L2(π2+CL2)h2≤μnμ1+C≤κ2(Bh)≤μn+Cμ1≤L2h2+CL24\frac{2L^2}{(\pi^2+CL^2)h^2}\le\frac{\mu_n}{\mu_1+C}\le\kappa_2(B_h)\le\frac{\mu_n+C}{\mu_1}\le\frac{L^2}{h^2}+\frac{CL^2}4 である。特にccとLLを固定してn→∞n\to\inftyとするとκ2(Bh)=Θ(h−2)\kappa_2(B_h)=\Theta(h^{-2})である。
  4. c=0c=0ならばκ2(Bh)=cot⁡2(π/(2(n+1)))\kappa_2(B_h)=\cot^2\bigl(\pi/(2(n+1))\bigr)であり、n→∞n\to\inftyのときπ2h2κ2(Bh)/(4L2)→1\pi^2h^2\kappa_2(B_h)/(4L^2)\to1である。

証明.(1)を示す。§E20.10 命題 5.1 (2)によりh−2Anh^{-2}A_nの固有値を非減少に並べた列はμ1<⋯<μn\mu_1<\dots<\mu_nである。D:=diag⁡(c1,…,cn)D:=\operatorname{diag}(c_1,\dots,c_n)と置くと、∥x∥2=1\|x\|_2=1を満たすx∈Rnx\in\R^nについて0≤xTDx≤C0\le x^{\mathsf T}Dx\le Cであるから

xT(h−2An)x≤xTBhx≤xT(h−2An)x+Cx^{\mathsf T}(h^{-2}A_n)x\le x^{\mathsf T}B_hx\le x^{\mathsf T}(h^{-2}A_n)x+C

である。dim⁡S=k\dim S=kを満たす各部分空間SSで単位ベクトルx∈Sx\in Sについての最大値をとり、さらにSSについての最小値をとると、§E20.7 定理 4.2をBhB_hとh−2Anh^{-2}A_nに適用してμk≤νk≤μk+C\mu_k\le\nu_k\le\mu_k+Cを得る。

(2)を示す。θ:=π/(2(n+1))=πh/(2L)\theta:=\pi/(2(n+1))=\pi h/(2L)と置くと0<θ≤π/40<\theta\le\pi/4である。sin⁡\sinは[0,π/2][0,\pi/2]上で凹でありsin⁡t≤t\sin t\le tであるから、2θ/π≤sin⁡θ≤θ2\theta/\pi\le\sin\theta\le\thetaであり、μ1=4h−2sin⁡2θ\mu_1=4h^{-2}\sin^2\thetaから第一の評価を得る。nπ/(2(n+1))=π/2−θn\pi/(2(n+1))=\pi/2-\thetaからμn=4h−2cos⁡2θ\mu_n=4h^{-2}\cos^2\thetaであり、1/2≤cos⁡2θ≤11/2\le\cos^2\theta\le1から第二の評価を得る。

(3)を示す。命題 2.2 (1)によりBhB_hは実対称かつ正定値であるから、§E20.5 命題 4.4 (5)によりκ2(Bh)=νn/ν1\kappa_2(B_h)=\nu_n/\nu_1である。(1)から内側の二つの不等式を得る。(2)によりμn/(μ1+C)≥2h−2/(π2/L2+C)\mu_n/(\mu_1+C)\ge2h^{-2}/(\pi^2/L^2+C)、(μn+C)/μ1≤(4h−2+C)L2/4(\mu_n+C)/\mu_1\le(4h^{-2}+C)L^2/4である。

(4)を示す。c=0c=0ではBh=h−2AnB_h=h^{-2}A_nであり、κ2(Bh)=μn/μ1=cos⁡2θ/sin⁡2θ=cot⁡2θ\kappa_2(B_h)=\mu_n/\mu_1=\cos^2\theta/\sin^2\theta=\cot^2\thetaである。n→∞n\to\inftyのときθ→0\theta\to0であり、θ2cot⁡2θ=(θ/sin⁡θ)2cos⁡2θ→1\theta^2\cot^2\theta=(\theta/\sin\theta)^2\cos^2\theta\to1、θ2=π2h2/(4L2)\theta^2=\pi^2h^2/(4L^2)である。▨

例 2.4.c=0c=0、a=0a=0、b=Lb=Lとする。§E12.13 例 4.2により、−y′′=λy-y''=\lambda y、y(0)=y(L)=0y(0)=y(L)=0の固有値はλk=(kπ/L)2\lambda_k=(k\pi/L)^2、固有関数はsin⁡(kπx/L)\sin(k\pi x/L)の定数倍である（k∈N≥1k\in\NN）。xj=jh=jL/(n+1)x_j=jh=jL/(n+1)であるから、§E20.10 命題 5.1 (2)の固有ベクトルs(k)s^{(k)}の成分sin⁡(jkπ/(n+1))\sin(jk\pi/(n+1))はsin⁡(kπx/L)\sin(k\pi x/L)のxjx_jでの値であり、h−2Ans(k)=μks(k)h^{-2}A_ns^{(k)}=\mu_ks^{(k)}である。1≤k≤n1\le k\le nについてθk:=kπh/(2L)∈(0,π/2)\theta_k:=k\pi h/(2L)\in(0,\pi/2)と置くと

μk=λk(sin⁡θkθk)2,4π2λk<μk<λk\mu_k=\lambda_k\Bigl(\frac{\sin\theta_k}{\theta_k}\Bigr)^2,\qquad\frac4{\pi^2}\lambda_k<\mu_k<\lambda_k

である。kkを固定してn→∞n\to\inftyとするとμk→λk\mu_k\to\lambda_kであり、μn/λn=(sin⁡θn/θn)2→4/π2\mu_n/\lambda_n=(\sin\theta_n/\theta_n)^2\to4/\pi^2である。

3 離散最大原理と収束

補題 3.1.n∈N≥1n\in\NN、h>0h>0、c1,…,cn≥0c_1,\dots,c_n\ge0とし、実数v0,…,vn+1v_0,\dots,v_{n+1}がv0≥0v_0\ge0、vn+1≥0v_{n+1}\ge0と

−vi−1−2vi+vi+1h2+civi≥0(1≤i≤n)-\frac{v_{i-1}-2v_i+v_{i+1}}{h^2}+c_iv_i\ge0\qquad(1\le i\le n)

を満たすとする。このとき、0≤i≤n+10\le i\le n+1を満たすすべてのiiについてvi≥0v_i\ge0である。

証明.m:=min⁡0≤i≤n+1vi<0m:=\min_{0\le i\le n+1}v_i<0と仮定し、vj=mv_j=mを満たす最小の添字をjjとする。v0,vn+1≥0>mv_0,v_{n+1}\ge0>mであるから1≤j≤n1\le j\le nである。jjの最小性とv0≥0v_0\ge0からvj−1>mv_{j-1}>mであり、vj+1≥mv_{j+1}\ge mである。したがって

−vj−1+2vj−vj+1=(vj−vj−1)+(vj−vj+1)<0-v_{j-1}+2v_j-v_{j+1}=(v_j-v_{j-1})+(v_j-v_{j+1})<0

であり、cj≥0c_j\ge0、vj<0v_j<0からcjvj≤0c_jv_j\le0である。二つを合わせると−(vj−1−2vj+vj+1)/h2+cjvj<0-(v_{j-1}-2v_j+v_{j+1})/h^2+c_jv_j<0であり、これはi=ji=jの仮定の不等式と両立しない。▨

命題 3.2.定義 2.1の設定でc≥0c\ge0とし、1≤i≤n1\le i\le nについてwi:=(xi−a)(b−xi)/2w_i:=(x_i-a)(b-x_i)/2、w:=(w1,…,wn)Tw:=(w_1,\dots,w_n)^{\mathsf T}と置く。

  1. r∈Rnr\in\R^nの成分がすべて非負ならば、Bh−1rB_h^{-1}rの成分もすべて非負である。
  2. 1≤i≤n1\le i\le nについて(Bhw)i=1+ciwi≥1(B_hw)_i=1+c_iw_i\ge1である。
  3. 任意のr∈Rnr\in\R^nと1≤i≤n1\le i\le nについて∣(Bh−1r)i∣≤∥r∥∞wi≤L2∥r∥∞/8|(B_h^{-1}r)_i|\le\|r\|_\infty w_i\le L^2\|r\|_\infty/8である。特に∥Bh−1∥∞≤L2/8\|B_h^{-1}\|_\infty\le L^2/8である。

証明.(1)を示す。v:=Bh−1rv:=B_h^{-1}rにv0:=vn+1:=0v_0:=v_{n+1}:=0を添えると、命題 2.2 (2)により1≤i≤n1\le i\le nについて−Dh2v(xi)+civi=ri≥0-D^2_hv(x_i)+c_iv_i=r_i\ge0である。補題 3.1によりvi≥0v_i\ge0である。

(2)を示す。W(x):=(x−a)(b−x)/2W(x):=(x-a)(b-x)/2は二次多項式であり、W′′=−1W''=-1、W(4)=0W^{(4)}=0である。§E20.19 定理 1.2 (3)により1≤i≤n1\le i\le nについてDh2W(xi)=W′′(xi)=−1D^2_hW(x_i)=W''(x_i)=-1である。W(x0)=W(xn+1)=0W(x_0)=W(x_{n+1})=0であるから、命題 2.2 (2)により(Bhw)i=1+ciwi(B_hw)_i=1+c_iw_iであり、ci≥0c_i\ge0、wi≥0w_i\ge0から(Bhw)i≥1(B_hw)_i\ge1である。

(3)を示す。∥r∥∞=:M\|r\|_\infty=:Mと置く。(2)によりBh(Mw±Bh−1r)=MBhw±rB_h(Mw\pm B_h^{-1}r)=MB_hw\pm rの各成分はM±ri≥0M\pm r_i\ge0以上である。(1)によりMw±Bh−1rMw\pm B_h^{-1}rの成分はすべて非負であり、∣(Bh−1r)i∣≤Mwi|(B_h^{-1}r)_i|\le Mw_iである。(x−a)(b−x)≤L2/4(x-a)(b-x)\le L^2/4からwi≤L2/8w_i\le L^2/8である。▨

例 3.3.κ>0\kappa>0、a=0a=0、b=1b=1、c=κ2c=\kappa^2、f=0f=0、α=1\alpha=1、β=0\beta=0とする。命題 1.1の解はu(x)=sinh⁡(κ(1−x))/sinh⁡κu(x)=\sinh(\kappa(1-x))/\sinh\kappaであり、0≤u≤10\le u\le1である。命題 1.2 (1)の関数はq(x)=sinh⁡(κx)/κq(x)=\sinh(\kappa x)/\kappaであり、命題 1.2 (2)により任意のs,δ∈Rs,\delta\in\Rについてys+δ−ys=δsinh⁡(κx)/κy_{s+\delta}-y_s=\delta\sinh(\kappa x)/\kappaである。命題 1.2 (3)により∣R(s)∣=∣s−s∗∣sinh⁡κ/κ|R(s)|=|s-s_*|\sinh\kappa/\kappaである。κ=40\kappa=40ではsinh⁡40/40=2.94…×1015\sinh40/40=2.94\ldots\times10^{15}であり、∣R(s)∣≤1|R(s)|\le1は∣s−s∗∣≤κ/sinh⁡κ=3.39…×10−16|s-s_*|\le\kappa/\sinh\kappa=3.39\ldots\times10^{-16}と同値である。

初期値をα+ϑ\alpha+\varthetaに替えると、y′′=κ2yy''=\kappa^2y、y(0)=1+ϑy(0)=1+\vartheta、y′(0)=sy'(0)=sの解はys+ϑcosh⁡(κx)y_s+\vartheta\cosh(\kappa x)である。この関数は方程式と初期条件を満たし、命題 1.2 (2)をα+ϑ\alpha+\varthetaに適用した一意性により解はこれに限る。κ=40\kappa=40、ϑ=10−16\vartheta=10^{-16}では終端値の変化は10−16cosh⁡40=11.76…10^{-16}\cosh40=11.76\ldotsである。

同じ係数の中心差分方程式で、境界値α\alphaをα+ϑ\alpha+\vartheta（ϑ≥0\vartheta\ge0）に替えた離散解から元の離散解を引いた差をVVとする。V0=ϑV_0=\vartheta、Vn+1=0V_{n+1}=0であり、1≤i≤n1\le i\le nについて−Dh2V(xi)+κ2Vi=0-D^2_hV(x_i)+\kappa^2V_i=0、−Dh2(ϑ−V)(xi)+κ2(ϑ−Vi)=κ2ϑ≥0-D^2_h(\vartheta-V)(x_i)+\kappa^2(\vartheta-V_i)=\kappa^2\vartheta\ge0である。補題 3.1をVVとϑ−V\vartheta-Vに適用すると、すべてのiiについて0≤Vi≤ϑ0\le V_i\le\varthetaである。命題 3.2 (3)の定数L2/8=1/8L^2/8=1/8もκ\kappaによらない。

定理 3.4.定義 2.1の設定でc≥0c\ge0とし、uuを命題 1.1の解、(U0,…,Un+1)(U_0,\dots,U_{n+1})を離散解とする。u∈C4([a,b];R)u\in C^4([a,b];\R)ならば、M4:=max⁡[a,b]∣u(4)∣M_4:=\max_{[a,b]}|u^{(4)}|と置くと、0≤i≤n+10\le i\le n+1について

∣Ui−u(xi)∣≤h2M424(xi−a)(b−xi)≤L2h2M496|U_i-u(x_i)|\le\frac{h^2M_4}{24}(x_i-a)(b-x_i)\le\frac{L^2h^2M_4}{96}

であり、任意のV^∈Rn\hat V\in\R^nについて

max⁡1≤i≤n∣V^i−u(xi)∣≤L2h2M496+L28∥F−BhV^∥∞\max_{1\le i\le n}|\hat V_i-u(x_i)|\le\frac{L^2h^2M_4}{96}+\frac{L^2}8\|F-B_h\hat V\|_\infty

である。c,f∈C2([a,b];R)c,f\in C^2([a,b];\R)ならばu∈C4([a,b];R)u\in C^4([a,b];\R)である。

証明.1≤i≤n1\le i\le nとする。[xi−h,xi+h]=[xi−1,xi+1]⊆[a,b][x_i-h,x_i+h]=[x_{i-1},x_{i+1}]\subseteq[a,b]であるから、§E20.19 定理 1.2 (3)によりξi∈[xi−1,xi+1]\xi_i\in[x_{i-1},x_{i+1}]が存在してDh2u(xi)=u′′(xi)+h2u(4)(ξi)/12D^2_hu(x_i)=u''(x_i)+h^2u^{(4)}(\xi_i)/12である。−u′′(xi)+ciu(xi)=fi-u''(x_i)+c_iu(x_i)=f_iから

−Dh2u(xi)+ciu(xi)=fi−h212u(4)(ξi)-D^2_hu(x_i)+c_iu(x_i)=f_i-\frac{h^2}{12}u^{(4)}(\xi_i)

である。Ei:=Ui−u(xi)E_i:=U_i-u(x_i)（0≤i≤n+10\le i\le n+1）と置くとE0=En+1=0E_0=E_{n+1}=0であり、τi:=h2u(4)(ξi)/12\tau_i:=h^2u^{(4)}(\xi_i)/12と置くと、中心差分方程式との差から−Dh2E(xi)+ciEi=τi-D^2_hE(x_i)+c_iE_i=\tau_iである。命題 2.2 (2)によりBh(E1,…,En)T=τB_h(E_1,\dots,E_n)^{\mathsf T}=\tauであり、∥τ∥∞≤h2M4/12\|\tau\|_\infty\le h^2M_4/12である。命題 3.2 (3)により∣Ei∣≤∥τ∥∞wi≤h2M4(xi−a)(b−xi)/24|E_i|\le\|\tau\|_\infty w_i\le h^2M_4(x_i-a)(b-x_i)/24であり、(xi−a)(b−xi)≤L2/4(x_i-a)(b-x_i)\le L^2/4から第一の評価を得る。

U^:=(U1,…,Un)T\hat U:=(U_1,\dots,U_n)^{\mathsf T}と置く。命題 2.2 (3)によりBhU^=FB_h\hat U=FであるからV^−U^=−Bh−1(F−BhV^)\hat V-\hat U=-B_h^{-1}(F-B_h\hat V)であり、命題 3.2 (3)により∥V^−U^∥∞≤L2∥F−BhV^∥∞/8\|\hat V-\hat U\|_\infty\le L^2\|F-B_h\hat V\|_\infty/8である。∣V^i−u(xi)∣≤∣V^i−Ui∣+∣Ei∣|\hat V_i-u(x_i)|\le|\hat V_i-U_i|+|E_i|から第二の評価を得る。

c,f∈C2c,f\in C^2ならば、u∈C2u\in C^2からu′′=cu−fu''=cu-fはC2C^2級であり、u∈C4u\in C^4である。▨

例 3.5.a=0a=0、b=1b=1、c=0c=0、u(x)=x4u(x)=x^4、f(x)=−12x2f(x)=-12x^2、α=0\alpha=0、β=1\beta=1とすると、uuは命題 1.1の解であり、M4=24M_4=24である。§E20.19 定理 1.2 (3)により1≤i≤n1\le i\le nについてDh2u(xi)−u′′(xi)=2h2D^2_hu(x_i)-u''(x_i)=2h^2であり、離散解との差Ei:=Ui−u(xi)E_i:=U_i-u(x_i)はE0=En+1=0E_0=E_{n+1}=0、−Dh2E(xi)=2h2-D^2_hE(x_i)=2h^2を満たす。命題 2.2 (2)によりBh(E1,…,En)T=2h2(1,…,1)TB_h(E_1,\dots,E_n)^{\mathsf T}=2h^2(1,\dots,1)^{\mathsf T}であり、c=0c=0では命題 3.2 (2)によりBhw=(1,…,1)TB_hw=(1,\dots,1)^{\mathsf T}である。BhB_hは正則であるから(E1,…,En)T=2h2w(E_1,\dots,E_n)^{\mathsf T}=2h^2w、すなわちEi=h2xi(1−xi)E_i=h^2x_i(1-x_i)である。nnが奇数ならばx(n+1)/2=1/2x_{(n+1)/2}=1/2であり、max⁡i∣Ei∣=h2/4=L2h2M4/96\max_i|E_i|=h^2/4=L^2h^2M_4/96である。したがって定理 3.4の第一の評価は等号で成り立つ。h=1/2,1/4,1/8,1/16h=1/2,1/4,1/8,1/16でのmax⁡i∣Ei∣\max_i|E_i|は1/16,1/64,1/256,1/10241/16,1/64,1/256,1/1024であり、格子幅を半分にするごとに1/41/4倍になる。

4 非線形の差分方程式と Newton 法

例 4.1.a<ba<b、L:=b−aL:=b-a、n∈N≥1n\in\NN、h:=L/(n+1)h:=L/(n+1)、xi:=a+ihx_i:=a+ihとし、u(x):=(x−a)(b−x)u(x):=(x-a)(b-x)、uh:=(u(x1),…,u(xn))Tu_h:=(u(x_1),\dots,u(x_n))^{\mathsf T}、gi:=2+u(xi)3g_i:=2+u(x_i)^3と置く。G ⁣:Rn→RnG\colon\R^n\to\R^nを

G(V):=h−2AnV+(V13,…,Vn3)T−(g1,…,gn)TG(V):=h^{-2}A_nV+(V_1^3,\dots,V_n^3)^{\mathsf T}-(g_1,\dots,g_n)^{\mathsf T}

で定める。命題 2.2 (2)をc=0c=0として用いると、V∈RnV\in\R^nにV0:=Vn+1:=0V_0:=V_{n+1}:=0を添えた列について、G(V)=0G(V)=0は−Dh2V(xi)+Vi3=gi-D^2_hV(x_i)+V_i^3=g_i（1≤i≤n1\le i\le n）と同値であり、これは−v′′+v3=2+u3-v''+v^3=2+u^3、v(a)=v(b)=0v(a)=v(b)=0の二階導関数を二階中心差分で置き換えた方程式である。

u′′=−2u''=-2、u(4)=0u^{(4)}=0、u(x0)=u(xn+1)=0u(x_0)=u(x_{n+1})=0であるから、§E20.19 定理 1.2 (3)と命題 2.2 (2)（c=0c=0）によりh−2Anuh=(2,…,2)Th^{-2}A_nu_h=(2,\dots,2)^{\mathsf T}であり、G(uh)=0G(u_h)=0である。GGはC1C^1級であり、DG(V)=h−2An+3diag⁡(V12,…,Vn2)DG(V)=h^{-2}A_n+3\operatorname{diag}(V_1^2,\dots,V_n^2)である。μ1:=4h−2sin⁡2(π/(2(n+1)))\mu_1:=4h^{-2}\sin^2(\pi/(2(n+1)))と置くと、§E20.7 定理 4.2をk=1k=1でh−2Anh^{-2}A_nに適用して、任意のx∈Rnx\in\R^nについてxTDG(V)x≥xTh−2Anx≥μ1∥x∥22x^{\mathsf T}DG(V)x\ge x^{\mathsf T}h^{-2}A_nx\ge\mu_1\|x\|_2^2である。したがってDG(V)DG(V)は実対称かつ正定値であり、Newton 法の各段の連立一次方程式は§E20.5 定理 6.1 (3)によりピボット選択なしの三重対角消去で解かれる。Cauchy–Schwarz の不等式から∥DG(V)x∥2≥μ1∥x∥2\|DG(V)x\|_2\ge\mu_1\|x\|_2であり、命題 2.3 (2)（c=0c=0）によりβ:=∥DG(uh)−1∥2≤1/μ1≤L2/4\beta:=\|DG(u_h)^{-1}\|_2\le1/\mu_1\le L^2/4である。

r:=L2/4r:=L^2/4とし、x,y∈Rnx,y\in\R^nが∥x−uh∥2≤r\|x-u_h\|_2\le r、∥y−uh∥2≤r\|y-u_h\|_2\le rを満たすとする。0≤u≤L2/40\le u\le L^2/4であるから∣xi∣,∣yi∣≤L2/2|x_i|,|y_i|\le L^2/2であり、

∥DG(x)−DG(y)∥2≤3max⁡i∣xi+yi∣∣xi−yi∣≤3L2∥x−y∥2\|DG(x)-DG(y)\|_2\le3\max_i|x_i+y_i||x_i-y_i|\le3L^2\|x-y\|_2

である。したがってuhu_hはGGの正則零点であり、半径rr、定数γ:=3L2\gamma:=3L^2の Lipschitz 条件を満たし、1/(2βγ)≥2/(3L4)1/(2\beta\gamma)\ge2/(3L^4)である。§E20.23 定理 2.5により、∥V0−uh∥2≤min⁡{L2/4, 2/(3L4)}\|V^0-u_h\|_2\le\min\{L^2/4,\ 2/(3L^4)\}を満たす任意のV0V^0から Newton 法の反復列(Vk)(V^k)が定まり、∥Vk+1−uh∥2≤βγ∥Vk−uh∥22≤3L44∥Vk−uh∥22\|V^{k+1}-u_h\|_2\le\beta\gamma\|V^k-u_h\|_2^2\le\frac{3L^4}4\|V^k-u_h\|_2^2を満たしてuhu_hに収束する。

a=0a=0、b=1b=1、n=7n=7とすると、この半径は1/41/4である。V0=0V^0=0は∥V0−uh∥2=0.516…\|V^0-u_h\|_2=0.516\ldotsを満たし、この半径の球に属さない。60桁の十進演算で計算すると、∥Vk−uh∥2\|V^k-u_h\|_2はk=1,2,3,4k=1,2,3,4で2.55×10−32.55\times10^{-3}、1.94×10−71.94\times10^{-7}、1.12×10−151.12\times10^{-15}、3.73×10−323.73\times10^{-32}である。V1V^1はこの半径の球に属し、k≥1k\ge1の値は評価∥Vk+1−uh∥2≤34∥Vk−uh∥22\|V^{k+1}-u_h\|_2\le\frac34\|V^k-u_h\|_2^2を満たす。

前提記事

11 本の記事・単元を表示