コンテンツにスキップ

有限要素法 I:Galerkin 近似・Céa の補題・1 次元の誤差評価

対象読者: 楕円型方程式の弱解(Lax–Milgram の定理)とヒルベルト空間の射影定理を学び、本テーマの差分法の章を読んだ修士課程の人。有限要素法は未習でかまいません。

前提となる記事:楕円型方程式の弱解と変分法:Lax–Milgram の定理・Gårding 不等式・Dirichlet 原理・直接法、Hilbert 空間の基礎:射影定理・リースの表現定理・正規直交基底、楕円型・放物型方程式の差分法:離散最大値原理と θ 法

生 Markdown

曲がった形の物体の熱の伝わり方や、橋の部材のたわみを計算機で求めるとき、答えの関数をそのまま扱うことはできません。そこで答えを「折れ線のような単純な関数のつなぎ合わせ」で近似し、その係数だけを連立一次方程式で求めます。この記事は、その近似がどれだけ正しいか(細かく刻むと誤差がどう小さくなるか)を、1 次元の問題で最後まで確かめます。読み終えると、有限要素法の誤差が刻み幅 hh のどんな冪で減るか、その理由を説明できるようになります。

  • 変分問題「すべての v∈Vv \in V で a(u,v)=f(v)a(u,v) = f(v)」を、有限次元の部分空間 Vh⊂VV_h \subset V の上で同じ式のまま解くのが Galerkin 近似です。aa が有界で強圧的なら、近似解 uhu_h はただ 1 つ存在し、hh によらず ∥uh∥≤∥f∥/α\|u_h\| \le \|f\|/\alpha です(安定性)。
  • 誤差は VhV_h のすべての元と「直交」します(Galerkin 直交性 a(u−uh,vh)=0a(u - u_h, v_h) = 0)。ここから Céa の補題 ∥u−uh∥≤Mαinf⁡vh∥u−vh∥\|u - u_h\| \le \frac{M}{\alpha}\inf_{v_h}\|u - v_h\| が出て、誤差評価が「VhV_h で uu をどれだけ近似できるか」という近似論の問題に帰着します。
  • 1 次元の区分 1 次要素では、節点補間の誤差が ∥(u−Ihu)′∥≤hπ∥u′′∥\|(u - I_hu)'\| \le \frac{h}{\pi}\|u''\|、∥u−Ihu∥≤h2π2∥u′′∥\|u - I_hu\| \le \frac{h^2}{\pi^2}\|u''\| です。したがって有限要素解の誤差は H1H^1 で O(h)O(h) です。
  • L2L^2 での誤差は、双対問題を使う Aubin–Nitsche の論法で O(h2)O(h^2) に上がります。−u′′=f-u'' = f では、有限要素解は節点で厳密解に一致します。
  • 等間隔の分割では、剛性行列は 3 点差分の行列の hh 倍です。剛性行列の条件数(連立方程式の解きにくさの目安)は O(h−2)O(h^{-2})、質量行列の条件数は O(1)O(1) です。

差分法は、微分方程式の導関数を差分商に置き換えました。前章の収束定理は解が C4C^4 級であることを仮定し、格子は座標軸に沿っていました。領域の形が複雑なとき、係数が不連続なとき、解が滑らかでないときには、どちらも不都合です。

有限要素法は、方程式ではなく解のほうを近似します。−u′′=f-u'' = f、u(0)=u(1)=0u(0) = u(1) = 0 に試験関数 vv を掛けて部分積分すると

∫01u′v′ dx=∫01fv dx(v(0)=v(1)=0)\int_0^1u'v'\,dx = \int_0^1fv\,dx\qquad(v(0) = v(1) = 0)

となります。これが弱形式(変分問題)で、uu に 1 階の導関数しか要求しません。ここで uu と vv を「折れ線」に限れば、有限個の未知数(節点での値)の連立一次方程式になります。

この方法の利点は、誤差解析の構造が単純なことです。差分法での「整合性 + 安定性 ⇒\Rightarrow 収束」は、有限要素法では「近似性 + 強圧性 ⇒\Rightarrow 収束」になります。安定性は連続の問題の強圧性からそのまま受け継がれ(部分空間に制限しても強圧性は壊れない)、整合性にあたるものは「部分空間が解をどれだけよく近似できるか」です。しかも Céa の補題は、有限要素解が部分空間の中の最良近似と定数倍しか違わないことを保証します。

この章では、まず抽象的な Hilbert 空間で枠組みを作り、次に 1 次元の区分 1 次要素で誤差解析を最後まで行います。弱解の一般論(Sobolev 空間、Lax–Milgram の定理、存在と一意性)は 楕円型方程式の弱解 にあります。2 次元の三角形要素は次章で扱います。

定義 2.1(有界で強圧的な双線形形式と変分問題)↗

VV を実 Hilbert 空間(ノルム ∥⋅∥V\|\cdot\|_V)とする。双線形形式 a:V×V→Ra : V \times V \to \mathbb{R} が、定数 M,α>0M, \alpha > 0 について

∣a(u,v)∣≤M∥u∥V∥v∥V,a(v,v)≥α∥v∥V2(u,v∈V)|a(u,v)| \le M\|u\|_V\|v\|_V,\qquad a(v,v) \ge \alpha\|v\|_V^2\qquad(u, v \in V)

を満たすとき、aa は有界で強圧的であるという。有界線形汎関数 f∈V∗f \in V^{*} に対し

u∈V,a(u,v)=f(v)(すべての v∈V)u \in V,\qquad a(u,v) = f(v)\quad(\text{すべての}\ v \in V)

を変分問題という。aa が対称(a(u,v)=a(v,u)a(u,v) = a(v,u))のとき、∥v∥a=a(v,v)1/2\|v\|_a = a(v,v)^{1/2} をエネルギーノルムという。

Lax–Milgram の定理(定理 3.2[楕円型方程式の弱解と変分法])から、変分問題の解 uu はただ 1 つ存在し、∥u∥V≤∥f∥V∗/α\|u\|_V \le \|f\|_{V^{*}}/\alpha です。aa が対称なら、a(⋅,⋅)a(\cdot,\cdot) は VV の内積で、α∥v∥V2≤∥v∥a2≤M∥v∥V2\alpha\|v\|_V^2 \le \|v\|_a^2 \le M\|v\|_V^2 からエネルギーノルムは ∥⋅∥V\|\cdot\|_V と同値です。

定義 2.2(Galerkin 近似と Ritz 射影)↗

Vh⊂VV_h \subset V を有限次元の部分空間とする。

uh∈Vh,a(uh,vh)=f(vh)(すべての vh∈Vh)u_h \in V_h,\qquad a(u_h, v_h) = f(v_h)\quad(\text{すべての}\ v_h \in V_h)

を満たす uhu_h を変分問題の Galerkin 近似という。w∈Vw \in V に対し、a(Rhw−w,vh)=0a(R_hw - w, v_h) = 0(すべての vh∈Vhv_h \in V_h)を満たす Rhw∈VhR_hw \in V_h を ww の Ritz 射影という。VhV_h の基底 φ1,…,φn\varphi_1, \ldots, \varphi_n について

K=(a(φj,φi))i,j,b=(f(φi))iK = \bigl(a(\varphi_j, \varphi_i)\bigr)_{i,j},\qquad \boldsymbol{b} = \bigl(f(\varphi_i)\bigr)_i

を剛性行列、荷重ベクトルという。

Ritz 射影は、右辺を v↦a(w,v)v \mapsto a(w, v) に置き換えた問題を VhV_h の中で解いた Galerkin 近似にあたります。あとの補題で、uhu_h がちょうど uu の Ritz 射影であることを確かめます。

uh=∑jcjφju_h = \sum_jc_j\varphi_j と書くと、Galerkin 近似の式は vh=φiv_h = \varphi_i(i=1,…,ni = 1, \ldots, n)で確かめれば十分で、連立一次方程式 Kc=bK\boldsymbol{c} = \boldsymbol{b} になります。添字 hh は部分空間の細かさを表すパラメータで、次節で分割の幅になります。

定理 2.3(Galerkin 近似の適切性)✓

aa が有界で強圧的ならば、任意の有限次元部分空間 Vh⊂VV_h \subset V と f∈V∗f \in V^{*} について、Galerkin 近似 uhu_h はただ 1 つ存在し

∥uh∥V≤1α∥f∥V∗\|u_h\|_V \le \frac{1}{\alpha}\|f\|_{V^{*}}

を満たす。剛性行列 KK は正則で、c≠0\boldsymbol{c} \ne \boldsymbol{0} なら cTKc>0\boldsymbol{c}^{\mathsf{T}}K\boldsymbol{c} > 0 である。aa が対称なら KK は正定値対称行列である。任意の w∈Vw \in V の Ritz 射影もただ 1 つ存在する。

証明(定理 2.3)

有限次元の部分空間は閉じているので、VhV_h は VV の内積で Hilbert 空間です。aa を Vh×VhV_h \times V_h に制限したものは、同じ定数 MM、α\alpha で有界かつ強圧的です。ff の VhV_h への制限 fhf_h は有界線形汎関数で、∥fh∥Vh∗=sup⁡vh≠0∣f(vh)∣/∥vh∥V≤∥f∥V∗\|f_h\|_{V_h^{*}} = \sup_{v_h \ne 0}|f(v_h)|/\|v_h\|_V \le \|f\|_{V^{*}} です。定理 3.2[楕円型方程式の弱解と変分法] を VhV_h に適用すると、uhu_h がただ 1 つ存在して ∥uh∥V≤∥fh∥/α≤∥f∥V∗/α\|u_h\|_V \le \|f_h\|/\alpha \le \|f\|_{V^{*}}/\alpha です。

c∈Rn\boldsymbol{c} \in \mathbb{R}^n、v=∑jcjφjv = \sum_jc_j\varphi_j とすると、双線形性から cTKc=∑i,jcia(φj,φi)cj=a(v,v)≥α∥v∥V2\boldsymbol{c}^{\mathsf{T}}K\boldsymbol{c} = \sum_{i,j}c_ia(\varphi_j, \varphi_i)c_j = a(v, v) \ge \alpha\|v\|_V^2 です。基底の一次独立性から c≠0\boldsymbol{c} \ne \boldsymbol{0} なら v≠0v \ne 0 で、cTKc>0\boldsymbol{c}^{\mathsf{T}}K\boldsymbol{c} > 0 です。特に Kc=0K\boldsymbol{c} = \boldsymbol{0} なら c=0\boldsymbol{c} = \boldsymbol{0} なので KK は正則です。aa が対称なら Kij=a(φj,φi)=a(φi,φj)=KjiK_{ij} = a(\varphi_j,\varphi_i) = a(\varphi_i,\varphi_j) = K_{ji} です。

Ritz 射影は、汎関数 v↦a(w,v)v \mapsto a(w, v)(∣a(w,v)∣≤M∥w∥V∥v∥V|a(w,v)| \le M\|w\|_V\|v\|_V から有界)を右辺とする Galerkin 近似なので、ただ 1 つ存在します。

評価 ∥uh∥V≤∥f∥/α\|u_h\|_V \le \|f\|/\alpha が VhV_h によらないことが、Galerkin 法の安定性です。差分法では安定性をスキームごとに確かめる必要がありましたが、ここでは連続の問題の強圧性から自動的に従います。

補題 2.4(Galerkin 直交性)✓

uu を変分問題の解、uhu_h をその Galerkin 近似とすると、すべての vh∈Vhv_h \in V_h について a(u−uh,vh)=0a(u - u_h, v_h) = 0 である。すなわち uh=Rhuu_h = R_hu である。

証明(補題 2.4)

Vh⊂VV_h \subset V なので、変分問題の式で v=vhv = v_h とでき、a(u,vh)=f(vh)a(u, v_h) = f(v_h) です。Galerkin 近似の式 a(uh,vh)=f(vh)a(u_h, v_h) = f(v_h) を引くと a(u−uh,vh)=0a(u - u_h, v_h) = 0 です。Ritz 射影の一意性(定理 2.3)から uh=Rhuu_h = R_hu です。

Vh⊂VV_h \subset V(適合性)が効いています。差分法の言葉では、厳密解が離散の方程式を残差なしで満たす、ということで、整合性が最も強い形で成り立っています。

定理 2.5(Céa の補題)✓

aa が有界で強圧的ならば

∥u−uh∥V≤Mαinf⁡vh∈Vh∥u−vh∥V\|u - u_h\|_V \le \frac{M}{\alpha}\inf_{v_h \in V_h}\|u - v_h\|_V

が成り立つ。aa が対称ならば、uhu_h はエネルギーノルムでの最良近似

∥u−uh∥a=min⁡vh∈Vh∥u−vh∥a\|u - u_h\|_a = \min_{v_h \in V_h}\|u - v_h\|_a

で、∥u−uh∥V≤Mαinf⁡vh∈Vh∥u−vh∥V\|u - u_h\|_V \le \sqrt{\dfrac{M}{\alpha}}\displaystyle\inf_{v_h \in V_h}\|u - v_h\|_V である。

この定理は、「有限要素解の誤差は、VhV_h の中で uu に最も近い元との距離の定数倍で抑えられる」という主張です。証明の方針は、Galerkin 直交性が使えるのは VhV_h の元を入れたときだけなので、VhV_h の元 uhu_h と任意の vhv_h の差 vh−uhv_h - u_h も VhV_h の元であることを使い、誤差 u−uhu - u_h を u−vhu - v_h の誤差に置き換えることです。

証明(定理 2.5)

vh∈Vhv_h \in V_h を任意にとります。vh−uh∈Vhv_h - u_h \in V_h なので、補題 2.4 から a(u−uh,vh−uh)=0a(u - u_h, v_h - u_h) = 0 です。強圧性と有界性から

α∥u−uh∥V2≤a(u−uh,u−uh)=a(u−uh,u−vh)+a(u−uh,vh−uh)=a(u−uh,u−vh)≤M∥u−uh∥V∥u−vh∥V\alpha\|u - u_h\|_V^2 \le a(u - u_h, u - u_h) = a(u - u_h, u - v_h) + a(u - u_h, v_h - u_h) = a(u - u_h, u - v_h) \le M\|u - u_h\|_V\|u - v_h\|_V

です。∥u−uh∥V≠0\|u - u_h\|_V \ne 0 なら両辺を割って ∥u−uh∥V≤Mα∥u−vh∥V\|u - u_h\|_V \le \frac{M}{\alpha}\|u - v_h\|_V で、00 なら自明です。vhv_h について下限をとれば前半です。

aa が対称なら aa は内積で、同じ計算と Cauchy–Schwarz の不等式から ∥u−uh∥a2=a(u−uh,u−vh)≤∥u−uh∥a∥u−vh∥a\|u - u_h\|_a^2 = a(u - u_h, u - v_h) \le \|u - u_h\|_a\|u - v_h\|_a、すなわち ∥u−uh∥a≤∥u−vh∥a\|u - u_h\|_a \le \|u - v_h\|_a です。vh=uhv_h = u_h で等号なので最小です。ノルムの同値性から

α∥u−uh∥V2≤∥u−uh∥a2≤∥u−vh∥a2≤M∥u−vh∥V2\alpha\|u - u_h\|_V^2 \le \|u - u_h\|_a^2 \le \|u - v_h\|_a^2 \le M\|u - v_h\|_V^2

で、後半を得ます。

V_hu_h(Galerkin 近似)v_hu(厳密解)誤差 u − u_hu − v_h
a が対称なときの Céa の補題の図(内積 a で測った長さ)。V_h を水平な線、解 u をその上の点とすると、u_h は u から V_h へ下ろした垂線の足で、誤差 u − u_h は V_h に直交する(Galerkin 直交性)。V_h の他のどの元 v_h と比べても、u_h が u に最も近い。

対称な場合、uhu_h は内積 aa についての uu の VhV_h への直交射影です(定理 3.4[Hilbert 空間の基礎] の最近点)。Céa の補題は、誤差評価を 2 つに分けます。方程式によるのは定数 M/αM/\alpha だけで、残りは「関数 uu を空間 VhV_h の元でどれだけ近似できるか」という、方程式と無関係な問題です。

定理 2.6(Ritz の最小化原理)✓

aa が対称、有界、強圧的とし、J(v)=12a(v,v)−f(v)J(v) = \frac12a(v,v) - f(v) とおく。

  1. すべての v∈Vv \in V について J(v)−J(u)=12∥v−u∥a2J(v) - J(u) = \frac12\|v - u\|_a^2。したがって uu は VV 上で JJ をただ 1 点で最小にする。
  2. すべての vh∈Vhv_h \in V_h について J(vh)−J(uh)=12∥vh−uh∥a2J(v_h) - J(u_h) = \frac12\|v_h - u_h\|_a^2。したがって uhu_h は VhV_h 上で JJ をただ 1 点で最小にする。
  3. J(uh)−J(u)=12∥u−uh∥a2≥0J(u_h) - J(u) = \frac12\|u - u_h\|_a^2 \ge 0。
証明(定理 2.6)

WW を VV または VhV_h、zz を WW での解(a(z,w)=f(w)a(z, w) = f(w)、w∈Ww \in W)とします。v∈Wv \in W について、対称性から a(v,v)=a(v−z,v−z)+2a(z,v−z)+a(z,z)a(v,v) = a(v - z, v - z) + 2a(z, v - z) + a(z,z) なので

J(v)−J(z)=12a(v−z,v−z)+a(z,v−z)−f(v−z)=12∥v−z∥a2J(v) - J(z) = \tfrac12a(v - z, v - z) + a(z, v - z) - f(v - z) = \tfrac12\|v - z\|_a^2

です(v−z∈Wv - z \in W なので a(z,v−z)=f(v−z)a(z, v-z) = f(v-z))。W=VW = V、z=uz = u で 1、W=VhW = V_h、z=uhz = u_h で 2 です。∥v−z∥a=0\|v - z\|_a = 0 は v=zv = z のときに限るので最小点は一意です。3 は 1 で v=uhv = u_h としたものです。

JJ は物理ではポテンシャルエネルギーにあたります。−Δu=f-\Delta u = f の場合の 1 が Dirichlet 原理(定理 5.1[楕円型方程式の弱解と変分法])です。3 は、Galerkin 近似のエネルギー J(uh)J(u_h) が常に真の値より大きく、その差がちょうどエネルギーノルムの誤差の 2 乗の半分であることを示しています。

注意 2.7(強圧的でない問題)

aa が強圧的でない問題(鞍点問題、Helmholtz 方程式など)や、試験関数の空間を解の空間と変える Petrov–Galerkin 法では、安定性は自動的には従いません。強圧性の代わりに inf-sup 条件 inf⁡uhsup⁡vha(uh,vh)∥uh∥∥vh∥≥β>0\inf_{u_h}\sup_{v_h}\frac{a(u_h, v_h)}{\|u_h\|\|v_h\|} \ge \beta > 0 を VhV_h ごとに確かめる必要があります(Babuška の定理)。また、移流が拡散よりずっと強い問題 −εu′′+u′=f-\varepsilon u'' + u' = f では、強圧性の定数 α=ε\alpha = \varepsilon が小さく M/αM/\alpha が大きくなり、標準的な Galerkin 法の解は境界層の近くで振動します。等間隔の区分 1 次要素では移流項 u′u' が中心差分 (ci+1−ci−1)/2(c_{i+1} - c_{i-1})/2 になることが原因で、前々章の移流方程式で風上差分を使ったのと同じく、風上化にあたる安定化が要ります。Ern–Guermond(参考文献)がこれらを体系的に扱っています。

3. 1 次元のモデル問題と有限要素空間

Section titled “3. 1 次元のモデル問題と有限要素空間”

区間 (0,1)(0,1) の上で、関数空間を具体的に定めます。∥⋅∥\|\cdot\| は L2(0,1)L^2(0,1) のノルム、⟨⋅,⋅⟩\langle\cdot,\cdot\rangle はその内積です。

定義 3.1(区間上の H¹・H¹₀・H²)↗

[0,1][0,1] 上の絶対連続な関数 uu で u′∈L2(0,1)u' \in L^2(0,1) となるものの全体を H1(0,1)H^1(0,1)、そのうち u(0)=u(1)=0u(0) = u(1) = 0 を満たすものの全体を H01(0,1)H^1_0(0,1) と書く。C1C^1 級で u′u' が絶対連続、u′′∈L2(0,1)u'' \in L^2(0,1) となる関数の全体を H2(0,1)H^2(0,1) と書く。H01(0,1)H^1_0(0,1) には内積 ⟨u,v⟩V=⟨u′,v′⟩\langle u, v\rangle_V = \langle u', v'\rangle とノルム ∥u∥V=∥u′∥\|u\|_V = \|u'\| を入れる。

絶対連続な関数はほとんどいたるところ微分可能で、u(y)−u(x)=∫xyu′ dtu(y) - u(x) = \int_x^yu'\,dt が成り立ちます(定理 6.9[符号付き測度と Radon–Nikodym の定理])。この定義は、弱微分による定義(定義 2.1[Sobolev 空間])と 1 次元では同じ空間を与えますが、以下の証明では上の具体的な定義だけを使います。

補題 3.2(H¹₀(0,1) の基本性質)✓

  1. u,v∈H1(0,1)u, v \in H^1(0,1) なら uvuv は絶対連続で、∫xyu′v dt=[uv]xy−∫xyuv′ dt\int_x^yu'v\,dt = \bigl[uv\bigr]_x^y - \int_x^yuv'\,dt(部分積分)。
  2. u∈H01(0,1)u \in H^1_0(0,1) なら max⁡∣u∣≤∥u′∥\max|u| \le \|u'\|、∥u∥≤1π∥u′∥\|u\| \le \frac{1}{\pi}\|u'\|。
  3. H01(0,1)H^1_0(0,1) は内積 ⟨u′,v′⟩\langle u', v'\rangle について Hilbert 空間である。
証明(補題 3.2)

1. uu、vv は連続なので有界で、∣u(y)v(y)−u(x)v(x)∣≤max⁡∣u∣ ∣v(y)−v(x)∣+max⁡∣v∣ ∣u(y)−u(x)∣|u(y)v(y) - u(x)v(x)| \le \max|u|\,|v(y) - v(x)| + \max|v|\,|u(y) - u(x)| です。したがって、互いに交わらない区間の有限列についての和をとれば、uu と vv の絶対連続性から uvuv の絶対連続性が従います。uu と vv が微分可能な点(ほとんどすべての点)で (uv)′=u′v+uv′(uv)' = u'v + uv' なので、定理 6.9[符号付き測度と Radon–Nikodym の定理] から [uv]xy=∫xy(u′v+uv′)dt[uv]_x^y = \int_x^y(u'v + uv')dt です。

2. u(x)=∫0xu′ dtu(x) = \int_0^xu'\,dt と Cauchy–Schwarz の不等式から ∣u(x)∣≤x ∥u′∥≤∥u′∥|u(x)| \le \sqrt{x}\,\|u'\| \le \|u'\| です。後半は 補題 6.1[無限次元の固有値の包含] を a=0a = 0、h=1h = 1 で使ったものです。

3. 内積であることは、∥u′∥=0\|u'\| = 0 なら 2 から u=0u = 0 であることから分かります。完備性を示します。(un)(u_n) を Cauchy 列とすると、(un′)(u_n') は L2(0,1)L^2(0,1) の Cauchy 列なので、ある g∈L2g \in L^2 に収束します。u(x)=∫0xg dtu(x) = \int_0^xg\,dt とおくと、uu は絶対連続で u′=gu' = g(a.e.)、u(0)=0u(0) = 0 です。u(1)=∫01g dt=lim⁡n∫01un′ dt=lim⁡n(un(1)−un(0))=0u(1) = \int_0^1g\,dt = \lim_n\int_0^1u_n'\,dt = \lim_n(u_n(1) - u_n(0)) = 0 なので u∈H01(0,1)u \in H^1_0(0,1) で、∥un−u∥V=∥un′−g∥→0\|u_n - u\|_V = \|u_n' - g\| \to 0 です。

モデル問題. p∈C1[0,1]p \in C^1[0,1]、q∈C[0,1]q \in C[0,1]、p(x)≥p0>0p(x) \ge p_0 > 0、q(x)≥0q(x) \ge 0、f∈L2(0,1)f \in L^2(0,1) とし、

−(pu′)′+qu=f(0<x<1),u(0)=u(1)=0-(pu')' + qu = f\quad(0 < x < 1),\qquad u(0) = u(1) = 0

を考えます。V=H01(0,1)V = H^1_0(0,1)、

a(u,v)=∫01(pu′v′+quv) dx,f(v)=∫01fv dxa(u,v) = \int_0^1(pu'v' + quv)\,dx,\qquad f(v) = \int_0^1fv\,dx

とおきます。p1=max⁡pp_1 = \max p、q1=max⁡qq_1 = \max q とすると、補題 3.2 の 2 から

∣a(u,v)∣≤p1∥u′∥∥v′∥+q1∥u∥∥v∥≤M∥u∥V∥v∥V,M=p1+q1π2,a(v,v)≥p0∥v∥V2,∣f(v)∣≤∥f∥π∥v∥V|a(u,v)| \le p_1\|u'\|\|v'\| + q_1\|u\|\|v\| \le M\|u\|_V\|v\|_V,\quad M = p_1 + \frac{q_1}{\pi^2},\qquad a(v,v) \ge p_0\|v\|_V^2,\qquad|f(v)| \le \frac{\|f\|}{\pi}\|v\|_V

なので、aa は対称、有界、強圧的(α=p0\alpha = p_0)で、∥f∥V∗≤∥f∥/π\|f\|_{V^{*}} \le \|f\|/\pi です。変分問題の解 u∈Vu \in V がモデル問題の弱解です(多次元の一般の定義は 定義 3.1[楕円型方程式の弱解と変分法]、存在と一意性は 定理 4.2[楕円型方程式の弱解と変分法])。

補題 3.3(1 次元の H² 正則性)✓

モデル問題の弱解 uu は H2(0,1)H^2(0,1) に属し、ほとんどいたるところ −(pu′)′+qu=f-(pu')' + qu = f を満たす。さらに

∥u′∥≤∥f∥πp0,∥u′′∥≤CR∥f∥,CR=1p0(1+q1/π+max⁡∣p′∣πp0)\|u'\| \le \frac{\|f\|}{\pi p_0},\qquad\|u''\| \le C_R\|f\|,\qquad C_R = \frac{1}{p_0}\Bigl(1 + \frac{q_1/\pi + \max|p'|}{\pi p_0}\Bigr)

が成り立つ。p=1p = 1、q=0q = 0 のときは u′′=−fu'' = -f で ∥u′′∥=∥f∥\|u''\| = \|f\| である。

この補題は、あとの誤差評価で必要な「u′′u'' の大きさを右辺 ff で抑える」ための準備です。証明の方針は、弱形式から pu′+Gpu' + G が定数になること(GG は gg の原始関数)を示し、そこから u′u' が絶対連続で、方程式が各点で成り立つことを読み取ることです。定数になることは、試験関数 vv をうまく選んで二乗の積分を 00 にして示します。

証明(補題 3.3)

∥u′∥=∥u∥V≤∥f∥V∗/α≤∥f∥/(πp0)\|u'\| = \|u\|_V \le \|f\|_{V^{*}}/\alpha \le \|f\|/(\pi p_0) です。g=f−qu∈L2g = f - qu \in L^2、G(x)=∫0xg dtG(x) = \int_0^xg\,dt とおくと、GG は絶対連続で G′=gG' = g(a.e.)です。弱形式は ∫01pu′v′ dx=∫01gv dx\int_0^1pu'v'\,dx = \int_0^1gv\,dx で、補題 3.2 の 1 と v(0)=v(1)=0v(0) = v(1) = 0 から ∫01gv dx=−∫01Gv′ dx\int_0^1gv\,dx = -\int_0^1Gv'\,dx なので

∫01(pu′+G)v′ dx=0(すべての v∈V)\int_0^1(pu' + G)v'\,dx = 0\qquad(\text{すべての}\ v \in V)

です。w=pu′+G∈L2w = pu' + G \in L^2、c=∫01w dxc = \int_0^1w\,dx とし、v(x)=∫0x(w−c) dtv(x) = \int_0^x(w - c)\,dt とおくと、vv は絶対連続で v(0)=0v(0) = 0、v(1)=c−c=0v(1) = c - c = 0、v′=w−c∈L2v' = w - c \in L^2 なので v∈Vv \in V です。この vv を代入すると

0=∫01w(w−c) dx=∫01(w−c)2dx+c∫01(w−c) dx=∫01(w−c)2dx0 = \int_0^1w(w - c)\,dx = \int_0^1(w - c)^2dx + c\int_0^1(w - c)\,dx = \int_0^1(w - c)^2dx

で、w=cw = c(a.e.)、すなわち u′=(c−G)/pu' = (c - G)/p(a.e.)です。右辺は連続なので、u(x)=∫0xu′ dt=∫0xc−Gpdtu(x) = \int_0^xu'\,dt = \int_0^x\frac{c - G}{p}dt は C1C^1 級で、すべての点で u′=(c−G)/pu' = (c - G)/p です。1/p1/p は C1C^1 級、c−Gc - G は絶対連続なので、補題 3.2 の 1 と同じ議論で u′u' は絶対連続です。pu′=c−Gpu' = c - G を微分して (pu′)′=−g(pu')' = -g(a.e.)、すなわち −(pu′)′+qu=f-(pu')' + qu = f です。u′′=−(g+p′u′)/pu'' = -(g + p'u')/p(a.e.)なので

∥u′′∥≤1p0(∥f∥+q1∥u∥+max⁡∣p′∣ ∥u′∥)≤1p0(1+q1/π+max⁡∣p′∣πp0)∥f∥\|u''\| \le \frac{1}{p_0}\bigl(\|f\| + q_1\|u\| + \max|p'|\,\|u'\|\bigr) \le \frac{1}{p_0}\Bigl(1 + \frac{q_1/\pi + \max|p'|}{\pi p_0}\Bigr)\|f\|

です(∥u∥≤∥u′∥/π\|u\| \le \|u'\|/\pi と ∥u′∥\|u'\| の評価を使いました)。p=1p = 1、q=0q = 0 なら g=fg = f、u′′=−fu'' = -f です。

多次元では、弱解が H2H^2 に属するかどうかは領域の形によります(次章)。1 次元では、方程式から直接 u′′u'' が求まるので、常に成り立ちます。

定義 3.4(1 次元の分割と有限要素空間)↗

0=x0<x1<⋯<xN=10 = x_0 < x_1 < \cdots < x_N = 1 を [0,1][0,1] の分割、ei=[xi−1,xi]e_i = [x_{i-1}, x_i] を要素、xix_i を節点、hi=xi−xi−1h_i = x_i - x_{i-1}、h=max⁡ihih = \max_ih_i とする。k=1,2k = 1, 2 について

Vhk={v∈C[0,1]:v∣ei は k 次以下の多項式 (1≤i≤N), v(0)=v(1)=0}V_h^k = \{v \in C[0,1] : v|_{e_i}\ \text{は}\ k\ \text{次以下の多項式}\ (1 \le i \le N),\ v(0) = v(1) = 0\}

を区分 kk 次の(連続な)有限要素空間という。Vh=Vh1V_h = V_h^1 と書く。φi∈Vh\varphi_i \in V_h(1≤i≤N−11 \le i \le N-1)を φi(xj)=δij\varphi_i(x_j) = \delta_{ij} で定め、節点基底(帽子関数)という。連続関数 uu(u(0)=u(1)=0u(0) = u(1) = 0)に対し

Ihu=∑i=1N−1u(xi)φiI_hu = \sum_{i=1}^{N-1}u(x_i)\varphi_i

を節点補間という。Vh2V_h^2 では、節点に各要素の中点を加えた点で同じように節点基底と補間 Ih(2)I_h^{(2)} を定める。G=(⟨φj,φi⟩)i,jG = (\langle\varphi_j, \varphi_i\rangle)_{i,j} を質量行列という(基底の Gram 行列。有界性の定数 MM と区別するため GG と書く)。定数 γ≥1\gamma \ge 1 について min⁡ihi≥h/γ\min_ih_i \ge h/\gamma となる分割を(定数 γ\gamma で)準一様という。

x₀ = 0x₁x₂x₃x₄x₅ = 1φ₁φ₂φ₄u と I_h u
区分 1 次の節点基底(帽子関数)と、関数 u の節点補間 I_h u。φ_i は節点 x_i で 1、他の節点で 0 で、隣り合う 2 つの要素の上でだけ 0 でない。このため剛性行列は三重対角になる。

区分的に多項式の連続関数は Lipschitz 連続なので絶対連続で、導関数は要素ごとの多項式の導関数(有界)です。したがって Vhk⊂H01(0,1)V_h^k \subset H^1_0(0,1) で、Galerkin 近似の枠組みに入ります。VhV_h の元は節点での値で決まるので dim⁡Vh=N−1\dim V_h = N - 1、vh=∑ivh(xi)φiv_h = \sum_iv_h(x_i)\varphi_i です。

剛性行列・質量行列・荷重ベクトルは、要素ごとの積分の和として計算します。要素 eie_i の上で 00 でない基底は φi−1\varphi_{i-1} と φi\varphi_i の 2 つだけなので、要素ごとに 2×22\times2 の行列(要素行列)を作り、全体の行列の対応する位置に足し込みます。これを組み立て(アセンブリ)といいます。疎行列の格納形式(定義 6.1[Cholesky 分解と帯行列・疎行列の直接法])のうち COO 形式は、同じ位置への重複した登録を足し合わせる約束にできるので、組み立てに向いています(§7 のコード)。

Céa の補題の右辺 inf⁡vh∥u−vh∥V\inf_{v_h}\|u - v_h\|_V を上から抑えるには、VhV_h の元を 1 つ選んで誤差を評価すれば十分です。節点補間 IhuI_hu がその役をします。

定理 4.1(1 次元の節点補間の誤差)✓

任意の分割について次が成り立つ。

  1. u∈H2(0,1)u \in H^2(0,1)、u(0)=u(1)=0u(0) = u(1) = 0 ならば ∥(u−Ihu)′∥≤hπ∥u′′∥\|(u - I_hu)'\| \le \dfrac{h}{\pi}\|u''\|、∥u−Ihu∥≤(hπ)2∥u′′∥\|u - I_hu\| \le \Bigl(\dfrac{h}{\pi}\Bigr)^2\|u''\|。
  2. uu が C2C^2 級で u′′u'' が絶対連続、u′′′∈L2(0,1)u''' \in L^2(0,1)、u(0)=u(1)=0u(0) = u(1) = 0 ならば、区分 2 次の節点補間について ∥(u−Ih(2)u)′∥≤h2∥u′′′∥\|(u - I_h^{(2)}u)'\| \le h^2\|u'''\|、∥u−Ih(2)u∥≤h32π∥u′′′∥\|u - I_h^{(2)}u\| \le \dfrac{h^3}{2\pi}\|u'''\|。
  3. k≥1k \ge 1 とし、uu は [0,1][0,1] を含む開区間で Ck+1C^{k+1} 級とする。各要素を kk 等分した k+1k+1 個の点での値を補間する区分 kk 次の連続関数を Ih(k)uI_h^{(k)}u とすると、max⁡[0,1]∣u−Ih(k)u∣≤hk+1(k+1)!max⁡∣u(k+1)∣\max_{[0,1]}|u - I_h^{(k)}u| \le \dfrac{h^{k+1}}{(k+1)!}\max|u^{(k+1)}|。

1 の直感は次のとおりです。補間誤差 w=u−Ihuw = u - I_hu は要素の両端で 00 なので、短い要素の上では大きくなれません(Poincaré 不等式)。しかも w′′=u′′w'' = u'' なので、ww の曲がり具合は u′′u'' で決まります。この 2 点から、要素の長さ ℓ\ell が小さいほど誤差が小さい、という評価が出ます。

証明(定理 4.1)

要素 e=[xi−1,xi]e = [x_{i-1}, x_i](長さ ℓ=hi≤h\ell = h_i \le h)の上で評価し、最後に要素について足します。∥⋅∥e\|\cdot\|_e は L2(e)L^2(e) のノルムです。

1. w=u−Ihuw = u - I_hu とおくと、ee の上で ww は C1C^1 級、w′w' は絶対連続で、IhuI_hu が 1 次式なので w′′=u′′w'' = u''(a.e.)、また w(xi−1)=w(xi)=0w(x_{i-1}) = w(x_i) = 0 です。補題 6.1[無限次元の固有値の包含] から

∥w∥e≤ℓπ∥w′∥e\|w\|_e \le \frac{\ell}{\pi}\|w'\|_e

です。次に w′w' を評価します。補題 3.2 の 1 の部分積分(区間 ee で、ww と w′w' に使う)と、両端で w=0w = 0 であることから

∥w′∥e2=∫ew′⋅w′ dx=[ww′]xi−1xi−∫eww′′ dx=−∫ewu′′ dx≤∥w∥e∥u′′∥e≤ℓπ∥w′∥e∥u′′∥e\|w'\|_e^2 = \int_ew'\cdot w'\,dx = \bigl[ww'\bigr]_{x_{i-1}}^{x_i} - \int_eww''\,dx = -\int_ewu''\,dx \le \|w\|_e\|u''\|_e \le \frac{\ell}{\pi}\|w'\|_e\|u''\|_e

で、∥w′∥e≤ℓπ∥u′′∥e\|w'\|_e \le \frac{\ell}{\pi}\|u''\|_e です。したがって ∥w∥e≤(ℓπ)2∥u′′∥e\|w\|_e \le (\frac{\ell}{\pi})^2\|u''\|_e です。ℓ≤h\ell \le h を使い、2 乗して要素について足すと、∥w′∥2=∑e∥w′∥e2≤(hπ)2∥u′′∥2\|w'\|^2 = \sum_e\|w'\|_e^2 \le (\frac{h}{\pi})^2\|u''\|^2、同様に ∥w∥2≤(hπ)4∥u′′∥2\|w\|^2 \le (\frac{h}{\pi})^4\|u''\|^2 です。

2. w=u−Ih(2)uw = u - I_h^{(2)}u は ee の両端と中点 mm の 3 点で 00 です。長さ ℓ/2\ell/2 の 2 つの半区間のそれぞれに 補題 6.1[無限次元の固有値の包含] を使って足すと ∥w∥e≤ℓ2π∥w′∥e\|w\|_e \le \frac{\ell}{2\pi}\|w'\|_e です。ww は ee の上で C2C^2 級なので、Rolle の定理から w′(ξ1)=w′(ξ2)=0w'(\xi_1) = w'(\xi_2) = 0 となる ξ1∈(xi−1,m)\xi_1 \in (x_{i-1}, m)、ξ2∈(m,xi)\xi_2 \in (m, x_i) があり、もう一度使って w′′(η)=0w''(\eta) = 0 となる η∈(ξ1,ξ2)\eta \in (\xi_1, \xi_2) があります。Ih(2)uI_h^{(2)}u は ee の上で 2 次式なので w′′′=u′′′w''' = u'''(a.e.)で、w′′(x)=∫ηxu′′′ dtw''(x) = \int_\eta^xu'''\,dt と Cauchy–Schwarz の不等式から ∣w′′(x)∣≤ℓ ∥u′′′∥e|w''(x)| \le \sqrt{\ell}\,\|u'''\|_e、したがって ∥w′′∥e≤ℓ∥u′′′∥e\|w''\|_e \le \ell\|u'''\|_e です。同じく w′(x)=∫ξ1xw′′ dtw'(x) = \int_{\xi_1}^xw''\,dt から ∥w′∥e≤ℓ∥w′′∥e≤ℓ2∥u′′′∥e\|w'\|_e \le \ell\|w''\|_e \le \ell^2\|u'''\|_e です。要素について足せば主張を得ます。

3. 要素 ee の上で、Ih(k)uI_h^{(k)}u は k+1k+1 個の点 y0,…,yk∈ey_0, \ldots, y_k \in e での uu の補間多項式です。定理 4.1[多項式補間] から、x∈ex \in e について u(x)−Ih(k)u(x)=u(k+1)(ξ)(k+1)!∏j(x−yj)u(x) - I_h^{(k)}u(x) = \frac{u^{(k+1)}(\xi)}{(k+1)!}\prod_j(x - y_j) で、∣x−yj∣≤ℓ≤h|x - y_j| \le \ell \le h なので主張を得ます。

1 の証明は、境界で 00 の関数についての Poincaré 不等式を 2 回使っています。1 回目は ww に、2 回目は部分積分を通して w′w' にです。定数 hπ\frac{h}{\pi} は、要素の上の −d2/dx2-d^2/dx^2 の最小固有値 (π/ℓ)2(\pi/\ell)^2 から来ています。証明は要素ごとなので、uu が全体では H2H^2 に属さなくても、各要素の上で C1C^1 級かつ u′u' が絶対連続で u′′∈L2u'' \in L^2 ならば、同じ評価が成り立ちます(係数が不連続な問題で使います)。

次数は「L2L^2 ノルムで hk+1h^{k+1}、導関数で hkh^k」で、導関数を 1 回とるごとに hh の冪が 1 つ減ります。kk を上げれば次数は上がりますが、それに見合う滑らかさ(u(k+1)u^{(k+1)})が要ります。

定理 5.1(エネルギーノルムでの誤差)✓

モデル問題(p∈C1p \in C^1、p≥p0>0p \ge p_0 > 0、q≥0q \ge 0、f∈L2f \in L^2)の弱解を uu、区分 1 次の有限要素解を uh∈Vhu_h \in V_h とすると

∥(u−uh)′∥≤Mp0 hπ∥u′′∥≤Mp0 CRπ h∥f∥,∥u−uh∥H1≤1+π−2 ∥(u−uh)′∥\|(u - u_h)'\| \le \sqrt{\frac{M}{p_0}}\,\frac{h}{\pi}\|u''\| \le \sqrt{\frac{M}{p_0}}\,\frac{C_R}{\pi}\,h\|f\|,\qquad\|u - u_h\|_{H^1} \le \sqrt{1 + \pi^{-2}}\,\|(u - u_h)'\|

が成り立つ。ここで ∥w∥H12=∥w∥2+∥w′∥2\|w\|_{H^1}^2 = \|w\|^2 + \|w'\|^2、M=p1+q1/π2M = p_1 + q_1/\pi^2、CRC_R は 補題 3.3 の定数である。

証明(定理 5.1)

補題 3.3 から u∈H2(0,1)u \in H^2(0,1)、∥u′′∥≤CR∥f∥\|u''\| \le C_R\|f\| です。aa は対称なので、定理 2.5 の後半を vh=Ihu∈Vhv_h = I_hu \in V_h で使い、定理 4.1 の 1 を代入すると

∥(u−uh)′∥=∥u−uh∥V≤Mp0∥u−Ihu∥V=Mp0∥(u−Ihu)′∥≤Mp0 hπ∥u′′∥\|(u - u_h)'\| = \|u - u_h\|_V \le \sqrt{\frac{M}{p_0}}\|u - I_hu\|_V = \sqrt{\frac{M}{p_0}}\|(u - I_hu)'\| \le \sqrt{\frac{M}{p_0}}\,\frac{h}{\pi}\|u''\|

です。H1H^1 ノルムについては、u−uh∈H01(0,1)u - u_h \in H^1_0(0,1) に 補題 3.2 の 2 を使って ∥u−uh∥2≤π−2∥(u−uh)′∥2\|u - u_h\|^2 \le \pi^{-2}\|(u - u_h)'\|^2 です。

この評価だけからは、L2L^2 ノルムの誤差も O(h)O(h) としか言えません。しかし補間誤差は L2L^2 で O(h2)O(h^2) でした。有限要素解でも O(h2)O(h^2) になることを、次の論法で示します。

定理 5.2(L² ノルムでの誤差(Aubin–Nitsche の双対論法))✓

定理 5.1 と同じ設定で

∥u−uh∥≤MCR hπ ∥(u−uh)′∥≤MCRMp0 (hπ)2∥u′′∥\|u - u_h\| \le MC_R\,\frac{h}{\pi}\,\|(u - u_h)'\| \le MC_R\sqrt{\frac{M}{p_0}}\,\Bigl(\frac{h}{\pi}\Bigr)^2\|u''\|

が成り立つ。

誤差そのものを右辺にした双対問題を考える理由は、∥e∥2=⟨e,e⟩\|e\|^2 = \langle e, e\rangle を a(e,z)a(e, z) の形に書き直したいからです。そうすると Galerkin 直交性で zz を z−Ihzz - I_hz に取り替えられ、補間誤差の評価が 1 回余分に使えます。

証明(定理 5.2)

e=u−uh∈Ve = u - u_h \in V とおきます。ee を右辺とする双対問題

z∈V,a(v,z)=⟨e,v⟩(すべての v∈V)z \in V,\qquad a(v, z) = \langle e, v\rangle\quad(\text{すべての}\ v \in V)

を考えます。aa は対称なので、これは右辺 f=e∈L2f = e \in L^2 のモデル問題そのもので、補題 3.3 から z∈H2(0,1)z \in H^2(0,1)、∥z′′∥≤CR∥e∥\|z''\| \le C_R\|e\| です。v=ev = e とし、補題 2.4 の a(e,Ihz)=0a(e, I_hz) = 0(Ihz∈VhI_hz \in V_h)を使うと

∥e∥2=a(e,z)=a(e,z−Ihz)≤M∥e′∥ ∥(z−Ihz)′∥≤M∥e′∥ hπ∥z′′∥≤MCR hπ ∥e′∥ ∥e∥\|e\|^2 = a(e, z) = a(e, z - I_hz) \le M\|e'\|\,\|(z - I_hz)'\| \le M\|e'\|\,\frac{h}{\pi}\|z''\| \le MC_R\,\frac{h}{\pi}\,\|e'\|\,\|e\|

です(定理 4.1 の 1 を zz に使いました)。∥e∥\|e\| で割って最初の不等式、定理 5.1 を代入して 2 番目の不等式を得ます。

双対論法の要点は、誤差 ee の L2L^2 ノルムを a(e,z)a(e, z) と書き、Galerkin 直交性で zz を z−Ihzz - I_hz に置き換えて、hh をもう 1 つ稼ぐことです。そのためには、双対問題の解が H2H^2 に属して ∥z′′∥≤CR∥e∥\|z''\| \le C_R\|e\| となること(H2H^2 正則性)が必要です。1 次元では常に成り立ちますが、多次元の角のある領域では成り立たないことがあり、そのとき L2L^2 の次数は 22 より下がります(次章)。

命題 5.3(−u″ = f では節点で厳密)✓

a(u,v)=∫01u′v′ dxa(u,v) = \int_0^1u'v'\,dx(p=1p = 1、q=0q = 0)とし、F∈V∗F \in V^{*} を任意の右辺、u∈Vu \in V を a(u,v)=F(v)a(u,v) = F(v) の解とする。任意の分割について、区分 1 次の有限要素解は uh=Ihuu_h = I_hu である。特に uh(xi)=u(xi)u_h(x_i) = u(x_i)(すべての節点)。

証明(命題 5.3)

vh∈Vhv_h \in V_h をとると、vh′v_h' は各要素 eie_i の上で定数 sis_i です。u−Ihuu - I_hu は絶対連続で、節点で 00 なので

a(u−Ihu,vh)=∑i=1Nsi∫ei(u−Ihu)′ dx=∑i=1Nsi[(u−Ihu)(xi)−(u−Ihu)(xi−1)]=0a(u - I_hu, v_h) = \sum_{i=1}^{N}s_i\int_{e_i}(u - I_hu)'\,dx = \sum_{i=1}^{N}s_i\bigl[(u - I_hu)(x_i) - (u - I_hu)(x_{i-1})\bigr] = 0

です。したがって a(Ihu,vh)=a(u,vh)=F(vh)a(I_hu, v_h) = a(u, v_h) = F(v_h) で、Ihu∈VhI_hu \in V_h は Galerkin 近似の式を満たします。一意性(定理 2.3)から uh=Ihuu_h = I_hu です。

等間隔の分割の場合が 補題 6.2[無限次元の固有値の包含] で、そこでは固有値の下界に使われています。背景には Green 関数があります。−u′′=f-u'' = f の Green 関数 gy(x)g_y(x)(yy を固定)は、x<yx < y で x(1−y)x(1-y)、x>yx > y で y(1−x)y(1-x) という折れ線で、yy が節点なら VhV_h に属します。u(y)−uh(y)=a(u−uh,gy)u(y) - u_h(y) = a(u - u_h, g_y) なので、Galerkin 直交性から節点で誤差が 00 になります。

例 5.4(係数が定数でないと節点で厳密にならない)

−u′′+u=1-u'' + u = 1、u(0)=u(1)=0u(0) = u(1) = 0 を、N=2N = 2(節点は x1=12x_1 = \frac12 だけ)で解きます。φ1\varphi_1 は 12\frac12 で 11 の帽子関数で、∫φ1′2=4\int\varphi_1'^2 = 4、∫φ12=13\int\varphi_1^2 = \frac13、∫φ1=12\int\varphi_1 = \frac12 です。Galerkin 近似は (4+13)c1=12(4 + \frac13)c_1 = \frac12 から uh(12)=c1=326=0.11538…u_h(\frac12) = c_1 = \frac{3}{26} = 0.11538\ldots です。厳密解は u(x)=1−cosh⁡(x−1/2)cosh⁡(1/2)u(x) = 1 - \frac{\cosh(x - 1/2)}{\cosh(1/2)} で u(12)=1−1cosh⁡(1/2)=0.11318…u(\frac12) = 1 - \frac{1}{\cosh(1/2)} = 0.11318\ldots です。一致しません。−u′′+u-u'' + u の Green 関数は sinh⁡\sinh で書かれ、折れ線ではないからです。

例 5.5(不連続な係数:分割を不連続点に合わせる)

pp が点 c∈(0,1)c \in (0,1) で跳ぶ区分的な定数(x<cx < c で 11、x>cx > c で p2≠1p_2 \ne 1)のとき、弱解 uu では pu′pu' が連続で、u′u' は cc で跳びます。uu は H2(0,1)H^2(0,1) に属しませんが、[0,c][0,c] と [c,1][c,1] のそれぞれの上では滑らかです。強圧性と有界性は成り立つので 定理 2.5 はそのまま使えます。

cc が節点なら、各要素の上で uu は滑らかなので、定理 4.1 の 1 の証明がそのまま通り、∥(u−Ihu)′∥≤hπ(∥u′′∥L2(0,c)2+∥u′′∥L2(c,1)2)1/2\|(u - I_hu)'\| \le \frac{h}{\pi}(\|u''\|_{L^2(0,c)}^2 + \|u''\|_{L^2(c,1)}^2)^{1/2} で、誤差は O(h)O(h) です。cc が要素 ee の内部にあると、ee の上で u′u' は跳び J≠0J \ne 0 をもつのに (Ihu)′(I_hu)' は定数なので、∥(u−Ihu)′∥e2\|(u - I_hu)'\|_e^2 は J2ℓJ^2\ell に比例する大きさになり、全体の誤差は O(h1/2)O(h^{1/2}) に落ちます。有限要素法では、係数や右辺の不連続点、境界の角を節点に合わせるのが原則です。

6. 行列の性質と差分法との関係

Section titled “6. 行列の性質と差分法との関係”

命題 6.1(等間隔の区分 1 次要素と 3 点差分)✓

等間隔の分割(hi=h=1/Nh_i = h = 1/N)で、a(u,v)=∫01u′v′ dxa(u,v) = \int_0^1u'v'\,dx とする。

  1. 剛性行列は K=1h tridiag(−1,2,−1)K = \frac{1}{h}\,\mathrm{tridiag}(-1, 2, -1)、質量行列は G=h6 tridiag(1,4,1)G = \frac{h}{6}\,\mathrm{tridiag}(1, 4, 1) である。KK は 3 点差分の行列 A=1h2tridiag(−1,2,−1)A = \frac{1}{h^2}\mathrm{tridiag}(-1,2,-1) の hh 倍である。
  2. 荷重ベクトルは bi=∫01fφi dx=hfˉib_i = \int_0^1f\varphi_i\,dx = h\bar f_i で、fˉi=1h∫fφi dx\bar f_i = \frac1h\int f\varphi_i\,dx は重み φi/h\varphi_i/h(積分が 11)による ff の平均である。
  3. bib_i を hf(xi)hf(x_i) で置き換える(質量行列の各行の和を対角に集める集中化にあたる)と、Kc=bK\boldsymbol{c} = \boldsymbol{b} は 3 点差分の方程式 −ci+1−2ci+ci−1h2=f(xi)-\frac{c_{i+1} - 2c_i + c_{i-1}}{h^2} = f(x_i) に一致する。
証明(命題 6.1)

1. φi′\varphi_i' は eie_i で 1h\frac1h、ei+1e_{i+1} で −1h-\frac1h、他で 00 です。Kii=∫φi′2=2⋅1h2⋅h=2hK_{ii} = \int\varphi_i'^2 = 2\cdot\frac{1}{h^2}\cdot h = \frac2h、Ki,i+1=∫ei+1φi′φi+1′=−1h2⋅h=−1hK_{i,i+1} = \int_{e_{i+1}}\varphi_i'\varphi_{i+1}' = -\frac{1}{h^2}\cdot h = -\frac1h、∣i−j∣≥2|i - j| \ge 2 では台が重ならないので 00 です。質量行列は、ei+1e_{i+1} の上で s=x−xis = x - x_i として φi=1−sh\varphi_i = 1 - \frac sh、φi+1=sh\varphi_{i+1} = \frac sh なので、Gii=2∫0h(sh)2ds=2h3G_{ii} = 2\int_0^h(\frac sh)^2ds = \frac{2h}{3}、Gi,i+1=∫0hsh(1−sh)ds=h2−h3=h6G_{i,i+1} = \int_0^h\frac sh(1 - \frac sh)ds = \frac h2 - \frac h3 = \frac h6 です。

2. ∫φi dx=h\int\varphi_i\,dx = h(底辺 2h2h、高さ 11 の三角形)なので、φi/h\varphi_i/h は積分が 11 の非負の重みです。

3. GG の第 ii 行の和は h6+4h6+h6=h\frac h6 + \frac{4h}{6} + \frac h6 = h です(両隣がともに内部の節点である行で。要素ごとに見れば、要素質量行列 h6(2112)\frac h6\bigl(\begin{smallmatrix}2&1\\1&2\end{smallmatrix}\bigr) を対角行列 h2I\frac h2I に置き換えることにあたり、これは台形則による積分です)。bi=hf(xi)b_i = hf(x_i) とすると、Kc=bK\boldsymbol{c} = \boldsymbol{b} の第 ii 式は 1h(−ci−1+2ci−ci+1)=hf(xi)\frac1h(-c_{i-1} + 2c_i - c_{i+1}) = hf(x_i) で、hh で割れば 3 点差分の式です。

つまり、等間隔の 1 次元では、有限要素法と差分法の違いは右辺の扱いだけです。有限要素法は ff の点での値の代わりに重みつき平均を使うので、ff が連続でなくても(L2L^2 でさえなくても、V∗V^{*} の元なら)意味をもちます。違いがはっきり出るのは、不等間隔の分割や多次元の不規則な分割で、そこでは差分法の式を書くこと自体が難しくなりますが、有限要素法は同じ手続き(要素ごとの積分と組み立て)で済みます。

連立方程式 Kc=bK\boldsymbol{c} = \boldsymbol{b} を反復法で解く手間は、条件数で決まります。次の補題は、その条件数を評価するための道具です。証明は、1 つの要素の上の 1 次式を両端の値 α\alpha、β\beta で表し、導関数の二乗の積分と関数の二乗の積分を計算して比べるだけです。

補題 6.2(1 次元の逆不等式と質量行列の評価)✓

vh∈Vhv_h \in V_h(区分 1 次)、c=(vh(xi))i=1N−1\boldsymbol{c} = (v_h(x_i))_{i=1}^{N-1}、hmin⁡=min⁡ihih_{\min} = \min_ih_i とする。

  1. 長さ ℓ\ell の要素 ee の上で ∥vh′∥e≤12ℓ∥vh∥e\|v_h'\|_e \le \dfrac{\sqrt{12}}{\ell}\|v_h\|_e。したがって ∥vh′∥≤12hmin⁡∥vh∥\|v_h'\| \le \dfrac{\sqrt{12}}{h_{\min}}\|v_h\|。
  2. hmin⁡3∣c∣2≤∥vh∥2=cTGc≤h∣c∣2\dfrac{h_{\min}}{3}|\boldsymbol{c}|^2 \le \|v_h\|^2 = \boldsymbol{c}^{\mathsf{T}}G\boldsymbol{c} \le h|\boldsymbol{c}|^2(∣⋅∣|\cdot| は Euclid ノルム)。
証明(補題 6.2)

要素 ee の両端での vhv_h の値を α\alpha、β\beta とします。vhv_h は ee の上で 1 次式なので、t∈[0,1]t \in [0,1] を ee の上の座標として vh=α(1−t)+βtv_h = \alpha(1 - t) + \beta t で、

∫evh′2dx=(β−α)2ℓ,∫evh2dx=ℓ∫01(α(1−t)+βt)2dt=ℓ3(α2+αβ+β2)\int_ev_h'^2dx = \frac{(\beta - \alpha)^2}{\ell},\qquad\int_ev_h^2dx = \ell\int_0^1\bigl(\alpha(1-t) + \beta t\bigr)^2dt = \frac{\ell}{3}(\alpha^2 + \alpha\beta + \beta^2)

です。

1. 4(α2+αβ+β2)−(β−α)2=3α2+6αβ+3β2=3(α+β)2≥04(\alpha^2 + \alpha\beta + \beta^2) - (\beta - \alpha)^2 = 3\alpha^2 + 6\alpha\beta + 3\beta^2 = 3(\alpha + \beta)^2 \ge 0 なので、∫evh′2≤4ℓ(α2+αβ+β2)=12ℓ2∫evh2\int_ev_h'^2 \le \frac{4}{\ell}(\alpha^2 + \alpha\beta + \beta^2) = \frac{12}{\ell^2}\int_ev_h^2 です。ℓ≥hmin⁡\ell \ge h_{\min} を使って要素について足せば後半です。

2. ∣αβ∣≤12(α2+β2)|\alpha\beta| \le \frac12(\alpha^2 + \beta^2) から ℓ6(α2+β2)≤∫evh2≤ℓ2(α2+β2)\frac{\ell}{6}(\alpha^2 + \beta^2) \le \int_ev_h^2 \le \frac{\ell}{2}(\alpha^2 + \beta^2) です。要素について足すと、内部の節点 xix_i の値 ci2c_i^2 は両側の 2 つの要素から現れ(境界の値は 00)、

∑i=1N−1hi+hi+16ci2≤∥vh∥2≤∑i=1N−1hi+hi+12ci2\sum_{i=1}^{N-1}\frac{h_i + h_{i+1}}{6}c_i^2 \le \|v_h\|^2 \le \sum_{i=1}^{N-1}\frac{h_i + h_{i+1}}{2}c_i^2

です。2hmin⁡≤hi+hi+1≤2h2h_{\min} \le h_i + h_{i+1} \le 2h から主張を得ます。∥vh∥2=∑i,jcicj⟨φj,φi⟩=cTGc\|v_h\|^2 = \sum_{i,j}c_ic_j\langle\varphi_j,\varphi_i\rangle = \boldsymbol{c}^{\mathsf{T}}G\boldsymbol{c} です。

1 は逆不等式と呼ばれます。一般の関数では導関数のノルムを関数のノルムで抑えることはできませんが、有限要素空間の中では、h−1h^{-1} を払えばできます。折れ線が最も激しく振動するのは節点ごとに符号を変えるとき(α=−β\alpha = -\beta)で、そのとき等号になります。

命題 6.3(剛性行列と質量行列の条件数)✓

モデル問題の双線形形式 aa と、定数 γ\gamma で準一様な分割について、2 ノルムの条件数(定義 6.2[行列のノルムと線形方程式の条件数])は

κ2(G)≤3γ,κ2(K)≤36Mγ3π2p0 h−2\kappa_2(G) \le 3\gamma,\qquad\kappa_2(K) \le \frac{36M\gamma^3}{\pi^2p_0}\,h^{-2}

を満たす。等間隔の分割で p=1p = 1、q=0q = 0 のときは κ2(K)=cot⁡2πh2≥4π2h2−1\kappa_2(K) = \cot^2\frac{\pi h}{2} \ge \frac{4}{\pi^2h^2} - 1 で、h−2h^{-2} の次数は改良できない。共役勾配法で Kc=bK\boldsymbol{c} = \boldsymbol{b} を解くと、エネルギーノルムの誤差を ε\varepsilon 倍にするのに k≥κ2(K)2ln⁡2εk \ge \frac{\sqrt{\kappa_2(K)}}{2}\ln\frac{2}{\varepsilon} 回の反復で足り、これは h−1h^{-1} に比例する。

証明(命題 6.3)

正定値対称行列 SS は、系 4.3[スペクトル定理] の正規直交固有基底で展開すると ∥S∥2=λmax⁡\|S\|_2 = \lambda_{\max}、∥S−1∥2=1/λmin⁡\|S^{-1}\|_2 = 1/\lambda_{\min} なので κ2(S)=λmax⁡/λmin⁡\kappa_2(S) = \lambda_{\max}/\lambda_{\min} です。また、すべての c≠0\boldsymbol{c} \ne \boldsymbol{0} で r∣c∣2≤cTSc≤R∣c∣2r|\boldsymbol{c}|^2 \le \boldsymbol{c}^{\mathsf{T}}S\boldsymbol{c} \le R|\boldsymbol{c}|^2 ならば、固有ベクトルを代入して r≤λmin⁡r \le \lambda_{\min}、λmax⁡≤R\lambda_{\max} \le R で、κ2(S)≤R/r\kappa_2(S) \le R/r です。

質量行列. 補題 6.2 の 2 から r=hmin⁡/3r = h_{\min}/3、R=hR = h で、κ2(G)≤3h/hmin⁡≤3γ\kappa_2(G) \le 3h/h_{\min} \le 3\gamma です。

剛性行列. vh=∑iciφiv_h = \sum_ic_i\varphi_i とすると cTKc=a(vh,vh)\boldsymbol{c}^{\mathsf{T}}K\boldsymbol{c} = a(v_h, v_h) です。下からは、強圧性、補題 3.2 の 2、補題 6.2 の 2 から

a(vh,vh)≥p0∥vh′∥2≥p0π2∥vh∥2≥p0π2hmin⁡3∣c∣2a(v_h, v_h) \ge p_0\|v_h'\|^2 \ge p_0\pi^2\|v_h\|^2 \ge \frac{p_0\pi^2h_{\min}}{3}|\boldsymbol{c}|^2

です。上からは、有界性と 補題 6.2 の 1、2 から

a(vh,vh)≤M∥vh′∥2≤12Mhmin⁡2∥vh∥2≤12Mhhmin⁡2∣c∣2a(v_h, v_h) \le M\|v_h'\|^2 \le \frac{12M}{h_{\min}^2}\|v_h\|^2 \le \frac{12Mh}{h_{\min}^2}|\boldsymbol{c}|^2

です。したがって κ2(K)≤36Mhπ2p0hmin⁡3≤36Mγ3π2p0h−2\kappa_2(K) \le \frac{36Mh}{\pi^2p_0h_{\min}^3} \le \frac{36M\gamma^3}{\pi^2p_0}h^{-2} です。

等間隔の場合. 命題 6.1 から K=hAK = hA で、条件数は定数倍で変わらないので、系 8.4[行列のノルムと線形方程式の条件数] から κ2(K)=cot⁡2πh2≥4π2h2−1\kappa_2(K) = \cot^2\frac{\pi h}{2} \ge \frac{4}{\pi^2h^2} - 1 です。

共役勾配法. 定理 5.4[Krylov 部分空間法 I] から ∥ek∥K≤2ρk∥e0∥K\|\boldsymbol{e}_k\|_K \le 2\rho^k\|\boldsymbol{e}_0\|_K、ρ=κ−1κ+1\rho = \frac{\sqrt\kappa - 1}{\sqrt\kappa + 1} です。s=1/κ∈(0,1]s = 1/\sqrt\kappa \in (0, 1] として ln⁡ρ=ln⁡(1−s)−ln⁡(1+s)≤−2s\ln\rho = \ln(1 - s) - \ln(1 + s) \le -2s(ln⁡(1−s)≤−s\ln(1-s) \le -s、ln⁡(1+s)≤s\ln(1+s) \le s から)なので ρk≤e−2k/κ\rho^k \le e^{-2k/\sqrt\kappa} で、k≥κ2ln⁡2εk \ge \frac{\sqrt\kappa}{2}\ln\frac2\varepsilon なら 2ρk≤ε2\rho^k \le \varepsilon です。κ2(K)=O(h−1)\sqrt{\kappa_2(K)} = O(h^{-1}) です。

条件数 h−2h^{-2} は、差分法のときと同じく、微分作用素が 2 階であることの反映です(最小固有値は hh によらず、最大固有値が逆不等式の定数 h−2h^{-2} で増える)。質量行列の条件数が hh によらないことは、放物型問題の陽的な時間積分や、L2L^2 射影の計算で使われます。

実際の計算では、aa や ff の積分を求積則で近似します。このとき解いているのは元の Galerkin 近似ではないので、Céa の補題に補正が要ります。次の命題の aha_h、fhf_h は、aa と ff の積分を求積則で置き換えたもの(積分を和で近似したもの)で、元の問題との違いはこの置き換えだけです。

命題 6.4(求積による誤差(第 1 Strang の補題))✓

aa は VV 上で有界(定数 MM)、uu は変分問題の解とする。Vh×VhV_h \times V_h 上の双線形形式 aha_h と VhV_h 上の線形汎関数 fhf_h が与えられ、定数 α~>0\tilde\alpha > 0 について ah(vh,vh)≥α~∥vh∥V2a_h(v_h, v_h) \ge \tilde\alpha\|v_h\|_V^2(vh∈Vhv_h \in V_h)とする。このとき ah(u~h,vh)=fh(vh)a_h(\tilde u_h, v_h) = f_h(v_h)(すべての vh∈Vhv_h \in V_h)を満たす u~h∈Vh\tilde u_h \in V_h がただ 1 つあり

∥u−u~h∥V≤inf⁡vh∈Vh[(1+Mα~)∥u−vh∥V+1α~sup⁡wh≠0∣a(vh,wh)−ah(vh,wh)∣∥wh∥V]+1α~sup⁡wh≠0∣f(wh)−fh(wh)∣∥wh∥V\|u - \tilde u_h\|_V \le \inf_{v_h \in V_h}\Bigl[\Bigl(1 + \frac{M}{\tilde\alpha}\Bigr)\|u - v_h\|_V + \frac{1}{\tilde\alpha}\sup_{w_h \ne 0}\frac{|a(v_h, w_h) - a_h(v_h, w_h)|}{\|w_h\|_V}\Bigr] + \frac{1}{\tilde\alpha}\sup_{w_h \ne 0}\frac{|f(w_h) - f_h(w_h)|}{\|w_h\|_V}

が成り立つ(上限は wh∈Vhw_h \in V_h についてとる)。

証明(命題 6.4)

ah(vh,vh)>0a_h(v_h, v_h) > 0(vh≠0v_h \ne 0)なので、定理 2.3 の証明と同じく係数行列は正則で、u~h\tilde u_h はただ 1 つあります。vh∈Vhv_h \in V_h を任意にとり、wh=u~h−vhw_h = \tilde u_h - v_h とおきます。a(u,wh)=f(wh)a(u, w_h) = f(w_h) を使って

α~∥wh∥V2≤ah(u~h,wh)−ah(vh,wh)=fh(wh)−ah(vh,wh)=[fh(wh)−f(wh)]+a(u−vh,wh)+[a(vh,wh)−ah(vh,wh)]\tilde\alpha\|w_h\|_V^2 \le a_h(\tilde u_h, w_h) - a_h(v_h, w_h) = f_h(w_h) - a_h(v_h, w_h) = \bigl[f_h(w_h) - f(w_h)\bigr] + a(u - v_h, w_h) + \bigl[a(v_h, w_h) - a_h(v_h, w_h)\bigr]

です。∣a(u−vh,wh)∣≤M∥u−vh∥V∥wh∥V|a(u - v_h, w_h)| \le M\|u - v_h\|_V\|w_h\|_V を使い、∥wh∥V\|w_h\|_V で割ると(wh=0w_h = 0 なら自明)

∥wh∥V≤1α~(M∥u−vh∥V+sup⁡wh∣a(vh,wh)−ah(vh,wh)∣∥wh∥V+sup⁡wh∣f(wh)−fh(wh)∣∥wh∥V)\|w_h\|_V \le \frac{1}{\tilde\alpha}\Bigl(M\|u - v_h\|_V + \sup_{w_h}\frac{|a(v_h, w_h) - a_h(v_h, w_h)|}{\|w_h\|_V} + \sup_{w_h}\frac{|f(w_h) - f_h(w_h)|}{\|w_h\|_V}\Bigr)

で、∥u−u~h∥V≤∥u−vh∥V+∥wh∥V\|u - \tilde u_h\|_V \le \|u - v_h\|_V + \|w_h\|_V と合わせて、vhv_h について下限をとれば主張を得ます。

誤差は「近似の誤差」と「aa と ff の整合性の誤差」の和になります。差分法の「整合性 + 安定性」の形がここで再び現れます(α~\tilde\alpha が安定性の定数)。求積則の精度が足りないと、整合性の項が近似の誤差より大きくなり、次数が落ちます。

例 6.5(求積の次数が足りない例)

a(u,v)=∫u′v′a(u,v) = \int u'v' は正確に計算し、荷重だけを各要素の左端点の値で近似します。fh(wh)=∑ef(xeL)∫ewh dxf_h(w_h) = \sum_ef(x_e^{\mathrm{L}})\int_ew_h\,dx(xeLx_e^{\mathrm{L}} は ee の左端)です。f∈C1f \in C^1 なら ∣f(x)−f(xeL)∣≤hmax⁡∣f′∣|f(x) - f(x_e^{\mathrm{L}})| \le h\max|f'| なので

∣f(wh)−fh(wh)∣=∣∑e∫e(f−f(xeL))wh dx∣≤hmax⁡∣f′∣∫01∣wh∣ dx≤hπmax⁡∣f′∣ ∥wh∥V|f(w_h) - f_h(w_h)| = \Bigl|\sum_e\int_e\bigl(f - f(x_e^{\mathrm{L}})\bigr)w_h\,dx\Bigr| \le h\max|f'|\int_0^1|w_h|\,dx \le \frac{h}{\pi}\max|f'|\,\|w_h\|_V

で、命題 6.4 の整合性の項は O(h)O(h) です。エネルギーノルムの誤差は O(h)O(h) のままですが、節点での厳密性と L2L^2 の O(h2)O(h^2) は保証されなくなります。§7 の表で、実際に L2L^2 と節点の誤差が O(h)O(h) に落ちることを確かめます。要素の中点で ff を評価する(中点則)か、両端の平均(台形則。命題 6.1 の 3 の集中化)にすれば、求積の誤差は O(h2)O(h^2) になり、次数は保たれます。求積則の誤差の一般的な表し方は Peano 核定理(定理 3.2[数値積分 I])です。

区分 1 次要素の剛性行列・質量行列・荷重ベクトルを要素ごとに組み立て、3 つの問題を解きます。どれも厳密解は u=sin⁡πxu = \sin\pi x です。(i) −u′′=π2sin⁡πx-u'' = \pi^2\sin\pi x、(ii) p=1+xp = 1 + x、q=1q = 1 の問題、(iii) (i) の荷重を左端点の 1 点則で近似したもの。積分は各要素で 5 点の Gauss 則を使います。

# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1(arm64)
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
GX, GW = np.polynomial.legendre.leggauss(5) # [-1,1] の 5 点 Gauss 則
def assemble(x, p, q, f, load="gauss"):
"""節点 x の区分 1 次要素で -(p u')' + q u = f の K, M, b を要素ごとに組み立てる。"""
N = len(x) - 1
rows, cols, kv, mv = [], [], [], []
b = np.zeros(N + 1)
for e in range(N):
a, c = x[e], x[e + 1]
ell = c - a
s = 0.5 * (a + c) + 0.5 * ell * GX # 求積点
w = 0.5 * ell * GW
phi = np.array([(c - s) / ell, (s - a) / ell]) # 2 つの形状関数
dphi = np.array([-1.0 / ell, 1.0 / ell])
for i in range(2):
for j in range(2):
rows.append(e + i); cols.append(e + j)
kv.append(np.sum(w * (p(s) * dphi[i] * dphi[j] + q(s) * phi[i] * phi[j])))
mv.append(np.sum(w * phi[i] * phi[j]))
if load == "gauss":
b[e + i] += np.sum(w * f(s) * phi[i])
else: # 左端点の 1 点則(次数不足)
b[e + i] += 0.5 * ell * f(a)
K = sp.coo_matrix((kv, (rows, cols)), shape=(N + 1, N + 1)).tocsr()
M = sp.coo_matrix((mv, (rows, cols)), shape=(N + 1, N + 1)).tocsr()
return K[1:-1, 1:-1], M[1:-1, 1:-1], b[1:-1] # Dirichlet 条件
def errors(N, p, q, f, u, du, load="gauss"):
x = np.linspace(0.0, 1.0, N + 1)
K, M, b = assemble(x, p, q, f, load)
c = np.concatenate([[0.0], spsolve(K.tocsc(), b), [0.0]])
e0 = e1 = 0.0
for e in range(N):
ell = x[e + 1] - x[e]
s = 0.5 * (x[e] + x[e + 1]) + 0.5 * ell * GX
w = 0.5 * ell * GW
uh = (c[e] * (x[e + 1] - s) + c[e + 1] * (s - x[e])) / ell
e0 += np.sum(w * (u(s) - uh)**2)
e1 += np.sum(w * (du(s) - (c[e + 1] - c[e]) / ell)**2)
return np.sqrt(e1), np.sqrt(e0), np.max(np.abs(c - u(x)))
pi = np.pi
u, du = lambda x: np.sin(pi * x), lambda x: pi * np.cos(pi * x)
one, zero = lambda x: 1.0 + 0 * x, lambda x: 0 * x
probs = [("-u''=f", one, zero, lambda x: pi**2 * np.sin(pi * x), "gauss"),
("p=1+x,q=1", lambda x: 1 + x, one,
lambda x: -pi * np.cos(pi * x) + ((1 + x) * pi**2 + 1) * np.sin(pi * x), "gauss"),
("-u''=f left", one, zero, lambda x: pi**2 * np.sin(pi * x), "left")]
for name, p, q, f, load in probs:
E = np.array([errors(N, p, q, f, u, du, load) for N in (10, 20, 40, 80, 160)])
print(name)
for k, lab in enumerate(("H1", "L2", "node")):
o = np.log2(E[:-1, k] / E[1:, k])
print(f" {lab:4s}" + " ".join(f"{v:9.2e}" for v in E[:, k])
+ " |" + " ".join(f"{v:5.2f}" for v in o))
print("H1 上界 (h/π)‖u''‖, h=1/10:", f"{0.1 / pi * pi**2 / np.sqrt(2):.3e}")
-u''=f
H1 2.01e-01 1.01e-01 5.04e-02 2.52e-02 1.26e-02 | 1.00 1.00 1.00 1.00
L2 6.36e-03 1.59e-03 3.98e-04 9.95e-05 2.49e-05 | 2.00 2.00 2.00 2.00
node 4.44e-16 1.11e-15 3.44e-15 2.66e-15 9.83e-14 |-1.32 -1.63 0.37 -5.20
p=1+x,q=1
H1 2.01e-01 1.01e-01 5.04e-02 2.52e-02 1.26e-02 | 1.00 1.00 1.00 1.00
L2 5.96e-03 1.49e-03 3.73e-04 9.32e-05 2.33e-05 | 2.00 2.00 2.00 2.00
node 8.93e-04 2.24e-04 5.62e-05 1.41e-05 3.52e-06 | 1.99 1.99 2.00 2.00
-u''=f left
H1 2.52e-01 1.26e-01 6.31e-02 3.15e-02 1.58e-02 | 1.00 1.00 1.00 1.00
L2 2.87e-02 1.25e-02 6.01e-03 2.97e-03 1.48e-03 | 1.19 1.06 1.02 1.00
node 4.25e-02 1.91e-02 8.93e-03 4.30e-03 2.11e-03 | 1.15 1.10 1.05 1.03
H1 上界 (h/π)‖u''‖, h=1/10: 2.221e-01

検算. (a) 参照解は製造解 sin⁡πx\sin\pi x です。(ii) の右辺は −(pu′)′+qu=−p′u′−pu′′+qu=−πcos⁡πx+(1+x)π2sin⁡πx+sin⁡πx-(pu')' + qu = -p'u' - pu'' + qu = -\pi\cos\pi x + (1+x)\pi^2\sin\pi x + \sin\pi x です。(b) 観測次数は H1H^1 で 11、L2L^2 で 22 で、定理 5.1 と 定理 5.2 に合います。(i) の H1H^1 誤差 0.2010.201(h=0.1h = 0.1)は、上界 hπ∥u′′∥=0.1π⋅π22=0.222\frac{h}{\pi}\|u''\| = \frac{0.1}{\pi}\cdot\frac{\pi^2}{\sqrt2} = 0.222 の内側です(p=1p = 1、q=0q = 0 では M=p0=1M = p_0 = 1)。(c) (i) の節点での誤差は丸めの大きさ(10−1610^{-16}〜10−1310^{-13}。次数の欄に意味はありません)で、命題 5.3 のとおりです。NN とともに少しずつ増えるのは、条件数 O(h−2)O(h^{-2})(命題 6.3)で丸め誤差が拡大されるためです。(ii) の節点での誤差は O(h2)O(h^2) で、厳密ではありません(例 5.4)。(d) 別の方法との照合として、(i) では uh=Ihuu_h = I_hu なので ∥(u−uh)′∥2=∥u′∥2−∑i(u(xi)−u(xi−1))2/h\|(u - u_h)'\|^2 = \|u'\|^2 - \sum_i(u(x_i) - u(x_{i-1}))^2/h と閉じた形で書けます。h=0.1h = 0.1 では ∥u′∥2=π2/2=4.9348\|u'\|^2 = \pi^2/2 = 4.9348、第 2 項は 1h∑i(sin⁡πxi−sin⁡πxi−1)2=4sin⁡2(πh/2)h⋅12h=4.8944\frac{1}{h}\sum_i(\sin\pi x_i - \sin\pi x_{i-1})^2 = \frac{4\sin^2(\pi h/2)}{h}\cdot\frac{1}{2h} = 4.8944 で、差の平方根は 0.2010.201 となり表と一致します。(e) (iii) では L2L^2 と節点の誤差が O(h)O(h) に落ち、例 6.5 の説明と合います。

次に、等間隔の分割での KK と GG の条件数と、共役勾配法(相対残差 10−810^{-8} で停止)の反復回数です。

# Python 3.9.6 / NumPy 2.0.2 / SciPy 1.13.1(arm64)
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import cg
print(" N cond(K) cot^2(pi h/2) cond(M) CG 反復 sqrt(cond)/2*ln(2e8)")
for N in (16, 32, 64, 128, 256):
h = 1.0 / N
K = sp.diags([-1, 2, -1], [-1, 0, 1], shape=(N - 1, N - 1)) / h
M = sp.diags([1, 4, 1], [-1, 0, 1], shape=(N - 1, N - 1)) * h / 6
ek, em = np.linalg.eigvalsh(K.toarray()), np.linalg.eigvalsh(M.toarray())
b = h * np.ones(N - 1)
b[: N // 3] *= -2.0 # 対称性のない右辺
it = [0]
x, info = cg(K.tocsr(), b, rtol=1e-8, atol=0.0,
callback=lambda xk: it.__setitem__(0, it[0] + 1))
kap = ek[-1] / ek[0]
print(f"{N:4d} {kap:10.1f} {1 / np.tan(np.pi * h / 2)**2:12.1f} "
f"{em[-1] / em[0]:9.3f} {it[0]:7d} {0.5 * np.sqrt(kap) * np.log(2e8):10.1f}")
N cond(K) cot^2(pi h/2) cond(M) CG 反復 sqrt(cond)/2*ln(2e8)
16 103.1 103.1 2.925 15 97.0
32 414.3 414.3 2.981 31 194.5
64 1659.4 1659.4 2.995 63 389.3
128 6639.5 6639.5 2.999 127 778.7
256 26560.1 26560.1 3.000 255 1557.5

検算. (a) KK の条件数は閉じた形 cot⁡2πh2\cot^2\frac{\pi h}{2} と一致し、NN を 2 倍にすると約 4 倍になります(h−2h^{-2})。(b) 質量行列(コードの M、本文の GG)の条件数は 33 に近づき、命題 6.3 の上界 3γ=33\gamma = 3 の内側です。G=h6tridiag(1,4,1)G = \frac h6\mathrm{tridiag}(1,4,1) の固有値は 補題 8.3[行列のノルムと線形方程式の条件数] から h6(4+2cos⁡kπh)\frac h6(4 + 2\cos k\pi h) で、比は 4+2cos⁡πh4−2cos⁡πh→3\frac{4 + 2\cos\pi h}{4 - 2\cos\pi h} \to 3 です。(c) 共役勾配法の反復回数は N−1N - 1(未知数の個数)で、h−1h^{-1} に比例します。上界 κ2ln⁡2ε\frac{\sqrt\kappa}{2}\ln\frac{2}{\varepsilon} よりずっと少ないのは、共役勾配法が高々未知数の個数の反復で厳密解に達するからです。1 次元の三重対角行列では 1 回の反復で情報が隣の節点にしか伝わらないので、NN に比例する回数は避けられません。

例 7.1(H² に属さない解と、特異点に寄せた分割)

u(x)=x0.6−xu(x) = x^{0.6} - x は H01(0,1)H^1_0(0,1) に属します(u′=0.6x−0.4−1u' = 0.6x^{-0.4} - 1、∫01x−0.8dx<∞\int_0^1x^{-0.8}dx < \infty)が、u′′=−0.24x−1.4u'' = -0.24x^{-1.4} は L2L^2 に属さないので H2H^2 には属しません。F(v)=∫01u′v′ dxF(v) = \int_0^1u'v'\,dx は VV 上の有界線形汎関数で、uu は a(u,v)=∫u′v′=F(v)a(u,v) = \int u'v' = F(v) の解です。命題 5.3 から有限要素解は uh=Ihuu_h = I_hu なので、誤差 ∥(u−Ihu)′∥\|(u - I_hu)'\| を直接計算できます。要素ごとに ∫e(u′−se)2=∫eu′2−(u(xi)−u(xi−1))2/hi\int_e(u' - s_e)^2 = \int_eu'^2 - (u(x_i) - u(x_{i-1}))^2/h_i(ses_e は IhuI_hu の傾き)で、∫u′2\int u'^2 は原始関数 0.360.2x0.2−2x0.6+x\frac{0.36}{0.2}x^{0.2} - 2x^{0.6} + x で求まります。分割は xi=(i/N)βx_i = (i/N)^\beta とし、β=1\beta = 1(等間隔)と、00 に節点を寄せた β=5\beta = 5、1212 を比べます。

# Python 3.9.6 / NumPy 2.0.2(arm64)
import numpy as np
al = 0.6 # u = x^al - x は H^2 に属さない
u = lambda x: x**al - x
F = lambda x: al**2 * x**(2 * al - 1) / (2 * al - 1) - 2 * x**al + x # ∫ u'^2 の原始関数
def h1_error(x):
"""節点補間(= 有限要素解)の H^1 半ノルム誤差。要素ごとに ∫u'^2 - (Δu)^2/ℓ。"""
ell, du = np.diff(x), np.diff(u(x))
return np.sqrt(np.sum(np.diff(F(x)) - du**2 / ell))
for beta in (1.0, 5.0, 12.0):
Ns = (16, 32, 64, 128, 256)
e = [h1_error((np.arange(N + 1) / N)**beta) for N in Ns]
o = [np.log2(e[i] / e[i + 1]) for i in range(4)]
print(f"beta={beta:4.1f} " + " ".join(f"{v:9.2e}" for v in e)
+ " |" + " ".join(f"{v:5.2f}" for v in o))
beta= 1.0 6.79e-01 6.33e-01 5.91e-01 5.51e-01 5.14e-01 | 0.10 0.10 0.10 0.10
beta= 5.0 2.81e-01 2.00e-01 1.42e-01 1.01e-01 7.13e-02 | 0.49 0.49 0.50 0.50
beta=12.0 2.11e-01 1.15e-01 6.12e-02 3.19e-02 1.64e-02 | 0.87 0.91 0.94 0.96

等間隔では次数は 0.10.1 で、ほとんど収束しません。最初の要素 [0,h][0, h] だけで ∫0h(u′−se)2dx\int_0^h(u' - s_e)^2dx が h0.2h^{0.2} に比例する大きさになるからです。u′u' の特異な部分 0.6x−0.40.6x^{-0.4} の 2 乗の積分は 1.8h0.21.8h^{0.2}、傾きは se≈h−0.4s_e \approx h^{-0.4} で se2h≈h0.2s_e^2h \approx h^{0.2} なので、差はおよそ 0.8h0.20.8h^{0.2} で、その平方根が h0.1h^{0.1} に比例します(N=16N = 16 で 0.8⋅16−0.2=0.68\sqrt{0.8\cdot16^{-0.2}} = 0.68。表の 0.6790.679 と合います)。β=5\beta = 5 では、最初の要素の長さが N−5N^{-5} なので同じ量が N−5×0.1=N−0.5N^{-5\times0.1} = N^{-0.5} になり、表の次数 0.50.5 と合います。β=12\beta = 12 では次数が 11 に近づき、滑らかな解と同じ速さが(節点の個数 NN について)回復します。解の特異性が分かっていれば、分割をそこに寄せることで次数を取り戻せます。特異性が分からないときに誤差の大きい要素を自動的に見つけて細分するのが、事後誤差評価と適応的な細分です(次章の注意)。

演習 8.1易

長さ ℓ\ell の要素の上の区分 1 次要素について、a(u,v)=∫u′v′a(u,v) = \int u'v' の要素剛性行列と要素質量行列を求めよ。それを使って、N=3N = 3 の等間隔の分割の剛性行列(2×22 \times 2)を組み立てよ。

解答

要素の上の 2 つの形状関数は ψ1=1−t\psi_1 = 1 - t、ψ2=t\psi_2 = t(t=(x−xi−1)/ℓt = (x - x_{i-1})/\ell)で、ψ1′=−1/ℓ\psi_1' = -1/\ell、ψ2′=1/ℓ\psi_2' = 1/\ell です。

Ke=1ℓ(1−1−11),Ge=ℓ6(2112)K_e = \frac{1}{\ell}\begin{pmatrix}1 & -1\\ -1 & 1\end{pmatrix},\qquad G_e = \frac{\ell}{6}\begin{pmatrix}2 & 1\\ 1 & 2\end{pmatrix}

です(∫01t2dt=13\int_0^1t^2dt = \frac13、∫01t(1−t)dt=16\int_0^1t(1-t)dt = \frac16 に ℓ\ell を掛ける)。N=3N = 3、ℓ=13\ell = \frac13 では、節点 x0,…,x3x_0, \ldots, x_3 について 4×44 \times 4 の行列に 3 つの Ke=3(1−1−11)K_e = 3\bigl(\begin{smallmatrix}1&-1\\-1&1\end{smallmatrix}\bigr) を足し込み、境界の節点 x0x_0、x3x_3 の行と列を除くと

K=(3+3−3−33+3)=(6−3−36)=1h(2−1−12)K = \begin{pmatrix}3 + 3 & -3\\ -3 & 3 + 3\end{pmatrix} = \begin{pmatrix}6 & -3\\ -3 & 6\end{pmatrix} = \frac1h\begin{pmatrix}2 & -1\\ -1 & 2\end{pmatrix}

で、命題 6.1 の 1 と一致します。

演習 8.2標準

ε>0\varepsilon > 0 とし、V=H01(0,1)V = H^1_0(0,1) 上の双線形形式 a(u,v)=∫01(εu′v′+u′v) dxa(u,v) = \int_0^1(\varepsilon u'v' + u'v)\,dx(−εu′′+u′=f-\varepsilon u'' + u' = f の弱形式)を考える。aa は対称でないが、有界で強圧的であることを示し、定数 MM、α\alpha と Céa の補題の定数 M/αM/\alpha を求めよ。ε→0\varepsilon \to 0 で何が起きるか。

解答

v∈Vv \in V について、補題 3.2 の 1 から ∫01v′v dx=[v2]01−∫01vv′ dx\int_0^1v'v\,dx = \bigl[v^2\bigr]_0^1 - \int_0^1vv'\,dx で、v(0)=v(1)=0v(0) = v(1) = 0 なので ∫01v′v dx=0\int_0^1v'v\,dx = 0 です。したがって a(v,v)=ε∥v′∥2a(v,v) = \varepsilon\|v'\|^2 で、α=ε\alpha = \varepsilon で強圧的です。有界性は、同じ補題の 2 から ∣a(u,v)∣≤ε∥u′∥∥v′∥+∥u′∥∥v∥≤(ε+1π)∥u′∥∥v′∥|a(u,v)| \le \varepsilon\|u'\|\|v'\| + \|u'\|\|v\| \le (\varepsilon + \frac1\pi)\|u'\|\|v'\| で、M=ε+1πM = \varepsilon + \frac1\pi です。定理 2.5 の前半(対称性を使わない形)から

∥(u−uh)′∥≤(1+1πε)inf⁡vh∥(u−vh)′∥\|(u - u_h)'\| \le \Bigl(1 + \frac{1}{\pi\varepsilon}\Bigr)\inf_{v_h}\|(u - v_h)'\|

です。ε→0\varepsilon \to 0 で定数が 1πε\frac{1}{\pi\varepsilon} のように大きくなり、hh が ε\varepsilon に比べて大きいあいだは、最良近似がよくても有限要素解がよいとは保証されません。実際、この領域では解に振動が現れます(注意 2.7)。

演習 8.3標準

f∈C[0,1]f \in C[0,1] とし、−u′′=f-u'' = f、u(0)=u(1)=0u(0) = u(1) = 0 の区分 1 次の有限要素解を uhu_h とする。max⁡[0,1]∣u−uh∣≤h28max⁡∣f∣\max_{[0,1]}|u - u_h| \le \frac{h^2}{8}\max|f| を示せ。

解答

ff が連続なので、補題 3.3 の証明から u′=c−Gu' = c - G(G(x)=∫0xfG(x) = \int_0^xf)は C1C^1 級で、u∈C2[0,1]u \in C^2[0,1]、u′′=−fu'' = -f です。命題 5.3 から uh=Ihuu_h = I_hu です。長さ ℓ\ell の要素 [xi−1,xi][x_{i-1}, x_i] の上で、IhuI_hu は 2 点 xi−1x_{i-1}、xix_i での uu の 1 次の補間多項式なので、定理 4.1[多項式補間] から(この定理は節点を含む開区間での 2 回微分可能性を仮定するので、uu を x<0x < 0 では u(0)+u′(0)x+12u′′(0)x2u(0) + u'(0)x + \frac12u''(0)x^2、x>1x > 1 でも同様の 2 次式で延長して [0,1][0,1] を含む開区間で C2C^2 級にしておきます。ξ\xi は [xi−1,xi][x_{i-1}, x_i] にとれるので、延長の仕方は結論に影響しません)

u(x)−Ihu(x)=u′′(ξ)2(x−xi−1)(x−xi)u(x) - I_hu(x) = \frac{u''(\xi)}{2}(x - x_{i-1})(x - x_i)

です(要素の内部の点 xx について。端点では両辺 00)。∣(x−xi−1)(x−xi)∣≤ℓ2/4|(x - x_{i-1})(x - x_i)| \le \ell^2/4(中点で最大)なので ∣u−uh∣≤ℓ28max⁡∣u′′∣≤h28max⁡∣f∣|u - u_h| \le \frac{\ell^2}{8}\max|u''| \le \frac{h^2}{8}\max|f| です。1 次元の −u′′=f-u'' = f では、最大値ノルムでも O(h2)O(h^2) が初等的に得られます。

  • S. C. Brenner, L. R. Scott『The Mathematical Theory of Finite Element Methods』3rd ed., Springer, 2008 — 第 0 章(1 次元のモデル問題と双対論法)、第 2 章(変分問題と Céa の補題)。
  • P. G. Ciarlet『The Finite Element Method for Elliptic Problems』North-Holland, 1978(SIAM Classics, 2002)— Galerkin 近似の抽象論、補間誤差、求積の影響。
  • G. Strang, G. J. Fix『An Analysis of the Finite Element Method』Prentice-Hall, 1973 — Ritz 法としての有限要素法、1 次元の誤差解析、条件数、求積と「変分法の罪」。
  • A. Ern, J.-L. Guermond『Theory and Practice of Finite Elements』Springer, 2004 — inf-sup 条件と Petrov–Galerkin 法、移流が支配的な問題の安定化。
  • 菊地文雄『有限要素法の数理』培風館、1994 — 日本語での数学的な基礎。
  • J. Céa, “Approximation variationnelle des problèmes aux limites”, Annales de l’Institut Fourier 14 (1964), 345–444 — Céa の補題の原論文。

このテーマの 100 本ノック(50 問)

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

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

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