# 有限要素法 II：三角形要素・補間誤差・L² 評価と放物型問題 > 多角形領域の三角形分割の上の区分 1 次要素について、補間誤差 O(h²)・O(h²/ρ) を定数つきで証明し、Poisson 方程式の誤差評価、L² での双対論法、逆不等式、余接の公式と離散最大値原理、熱方程式の半離散・全離散の誤差評価を導く。 > https://rikai.mugen-giken.com/computer-science/pde-discretization/finite-element-triangles ## 0. この記事の要点 板の上の熱の広がりや、膜のたわみのような「平面の上の現象」を計算機で解くとき、領域を小さな三角形に分けて、各三角形の上では 1 次式（傾いた平面）で近似する方法があります。この記事は、そうして得た近似解が真の解からどれだけ離れるかを、三角形の大きさ $h$ を小さくしたときの速さとして数式で保証します。読み終えると、「三角形をつぶれた形にしてはいけない理由」「角のへこんだ領域で精度が落ちる理由」を、定数まで追って説明できるようになります。 - 多角形領域を三角形に分け、各三角形の上で 1 次式になる連続関数の空間 $V_h$ をとると、$V_h \subset H^1_0(\Omega)$ で、前章の Galerkin 近似の枠組みがそのまま使えます。 - 三角形 $T$（直径 $h_T$、内接円の直径 $\rho_T$）の上の節点補間は $\|u - I_hu\|_{L^2(T)} \le h_T^2|u|_{H^2(T)}$、$|u - I_hu|_{H^1(T)} \le 3\frac{h_T^2}{\rho_T}|u|_{H^2(T)}$ を満たします。$H^1$ の評価には三角形がつぶれていないこと（**形状正則性** $h_T/\rho_T \le \sigma$）が要ります。 - したがって、解が $H^2$ に属せば有限要素解の誤差は $H^1$ で $O(h)$ です。$L^2$ で $O(h^2)$ になるには、双対問題（誤差そのものを右辺にした補助の方程式）の $H^2$ 正則性（解が 2 階微分まで 2 乗可積分になること）が要ります。正方形では成り立ち（正弦級数で証明します）、L 字形の領域では成り立たず、次数が $2/3$ と $4/3$ に落ちます。 - 要素の剛性行列の成分は向かい合う角の余接で書けます（$-\frac12\cot\theta$）。辺に向かい合う 2 つの角の和が $\pi$ 以下なら非対角成分は $0$ 以下で、離散最大値原理（近似解の最大値が境界でとられ、内部に余計な山ができないこと）が成り立ちます。直角二等辺三角形の分割では剛性行列が 5 点差分に一致します。 - 熱方程式では、誤差を「Ritz 射影（有限要素空間の中で、エネルギーの意味でいちばん解に近い関数をとること）の誤差」と「有限要素空間の中の差」に分ける方法（Wheeler）で、半離散解が $L^2$ で $O(h^2)$、陰的 Euler・Crank–Nicolson による全離散解が $O(h^2 + \Delta t)$・$O(h^2 + \Delta t^2)$ であることが分かります。 ## 1. 動機 [前章](/computer-science/pde-discretization/finite-element-galerkin)では、Galerkin 近似の誤差が「有限要素空間で解をどれだけ近似できるか」に帰着すること（Céa の補題）を示し、1 次元でその近似誤差を評価しました。2 次元以上で有限要素法が力を発揮するのは、領域を三角形に分けることで、どんな多角形にも同じ手続きが使えるからです。 2 次元では、1 次元になかった問題が 3 つ現れます。 1. **点での値.** 節点補間には節点での値が要りますが、2 次元の $H^1$ 関数は連続とは限りません。$H^2$ なら連続です。これを三角形の上で、定数を明示して示します。 2. **要素の形.** 区間には長さしかありませんが、三角形には形があります。細長くつぶれた三角形では補間した関数の勾配が悪くなります。 3. **領域の角.** 1 次元では弱解は常に $H^2$ に属しました。2 次元では、領域の角が凹んでいると解が $H^2$ に属さず、収束の次数が落ちます。 この章の誤差評価は、三角形の上の関数を「平均をとった 1 次の Taylor 多項式」と剰余に分ける方法で行います。この方法だと定数がすべて明示でき、三角形の形への依存がどこから来るかが見えます。弱解と Poincaré の不等式は [楕円型方程式の弱解](/mathematics/pde/elliptic-weak-solutions)、正則性の一般論は [楕円型方程式の正則性](/mathematics/pde/elliptic-regularity) にあります。 ## 2. 三角形分割と区分 1 次の有限要素空間 以下 $\Omega \subset \mathbb{R}^2$ は有界な**多角形領域**（境界が有限個の線分からなる連結な開集合）です。$\|\cdot\|$ と $\langle\cdot,\cdot\rangle$ は $L^2(\Omega)$ のノルムと内積、$H^1(\Omega)$、$H^1_0(\Omega)$、$H^2(\Omega)$ は弱微分による Sobolev 空間（）です。集合 $S$ の上で $$ |u|_{H^1(S)} = \|\nabla u\|_{L^2(S)},\qquad|u|_{H^2(S)} = \|D^2u\|_{L^2(S)} = \Bigl(\sum_{i,j=1}^{2}\|\partial_i\partial_ju\|_{L^2(S)}^2\Bigr)^{1/2} $$ と書きます。 閉じた三角形（内部が空でないもの）の有限集合 $\mathcal{T}_h$ が $\bar\Omega$ の（適合する）**三角形分割**であるとは、$\bar\Omega = \bigcup_{T\in\mathcal{T}_h}T$ で、相異なる 2 つの三角形の共通部分が、空集合か、共通の頂点 1 つか、共通の辺 1 つの全体であることをいう。三角形の頂点を**節点**、$\partial\Omega$ の上にない節点を**内部の節点**という。$T \in \mathcal{T}_h$ について $$ h_T = \operatorname{diam}T\ (\text{最長の辺の長さ}),\qquad\rho_T = T\ \text{の内接円の直径},\qquad h = \max_Th_T $$ とおく。分割の族 $\{\mathcal{T}_h\}$ が、$h$ によらない定数 $\sigma$ について $h_T/\rho_T \le \sigma$（すべての $T$）を満たすとき**形状正則**、さらに定数 $\gamma$ について $h_T \ge h/\gamma$ を満たすとき**準一様**であるという。 適合性は、ある三角形の頂点が隣の三角形の辺の途中に来ること（ぶら下がった節点）を禁じます。形状正則性は、最小の角が $0$ に近づかないことと同値です。内接円の中心から頂点までの距離は $\frac{\rho_T/2}{\sin(\theta/2)}$（$\theta$ はその頂点の角）で $h_T$ 以下なので、$\sin\frac\theta2 \ge \frac{1}{2\sigma}$ です。したがって 1 つの節点を共有する三角形の個数は $n_\sigma = \pi/\arcsin\frac{1}{2\sigma}$ 以下です。 頂点 $a_0, a_1, a_2$ の三角形 $T$ について、$\lambda_i(a_j) = \delta_{ij}$ を満たす 1 次式 $\lambda_0, \lambda_1, \lambda_2$ を**重心座標**という。$\lambda_0 + \lambda_1 + \lambda_2 = 1$ で、$T$ の上で $\lambda_i \ge 0$ である。 $$ X_h = \{v \in C(\bar\Omega) : v|_T\ \text{は 1 次式}\ (T \in \mathcal{T}_h)\},\qquad V_h = \{v \in X_h : \partial\Omega\ \text{の上で}\ v = 0\} $$ を**区分 1 次の（連続な）有限要素空間**という。節点 $z_i$ ごとに $\varphi_i(z_j) = \delta_{ij}$ を満たす $\varphi_i \in X_h$ を**節点基底**という。$\bar\Omega$ 上の連続関数 $u$ に対し $I_hu = \sum_iu(z_i)\varphi_i \in X_h$ を**節点補間**という。$T$ の上では $I_hu = \sum_{i=0}^{2}u(a_i)\lambda_i$ である。 1 次式は 3 頂点での値で決まり、辺の上では両端の値で決まります。したがって、節点での値を任意に与えると、三角形ごとの 1 次式が辺の上で一致して連続な関数になり、$X_h$ の元がただ 1 つ定まります。$V_h$ の次元は内部の節点の個数です。$u$ が $\partial\Omega$ の上で $0$ なら $I_hu \in V_h$ です。$T$ の上で $\varphi_i$ は重心座標 $\lambda_i$ に一致するので、$T$ の上で $|I_hu| \le \max_j|u(a_j)|$ です。
適合する分割 適合しない分割（ぶら下がった節点）
次の補題は、三角形の境目で値がつながってさえいれば、区分的に滑らかな関数は折れ目をまたいでも弱微分をもつ、ということを述べています。証明は 2 段で、まず水平な線に沿った 1 変数の部分積分で弱微分を確かめ、次に、境界で $0$ の関数の絶対値の小さい部分を切り落とした関数で近似して $H^1_0(\Omega)$ に入ることを示します。 $v \in C(\bar\Omega)$ が、三角形分割 $\mathcal{T}_h$ の各三角形 $T$ の上で $C^1$ 級（$T$ の内部で $C^1$ 級で、$\nabla v$ が $T$ まで連続に延びる）とする。このとき $v \in H^1(\Omega)$ で、弱微分は三角形ごとの導関数に（ほとんどいたるところ）一致する。さらに $\partial\Omega$ の上で $v = 0$ ならば $v \in H^1_0(\Omega)$ である。特に $X_h \subset H^1(\Omega)$、$V_h \subset H^1_0(\Omega)$ である。 **弱微分.** $g_1$ を、各三角形の内部で $\partial_1v$ に等しい関数とします（辺の上は測度 $0$ なので定めなくてよく、$g_1$ は有界です）。$\varphi \in C_c^\infty(\Omega)$ をとり、$w = v\varphi$ を $\Omega$ の外で $0$ とおくと、$w$ は $\mathbb{R}^2$ 上の連続関数で台はコンパクトです。$x_2$ を、どの節点の第 2 座標とも異なる値に固定します（除かれる $x_2$ は有限個）。水平線 $\mathbb{R}\times\{x_2\}$ は、水平な辺を含まず、各辺と高々 1 点で交わります。したがって $x_1 \mapsto w(x_1, x_2)$ は連続で、有限個の点を除いて $C^1$ 級、導関数は $g_1\varphi + v\,\partial_1\varphi$（$\Omega$ の外では $0$）です。区分的に $C^1$ 級の連続関数には微積分学の基本定理が区間ごとに使え、連続性から端点の値が打ち消し合うので、台がコンパクトであることから $$ \int_{\mathbb{R}}\bigl(g_1\varphi + v\,\partial_1\varphi\bigr)(x_1, x_2)\,dx_1 = 0 $$ です。被積分関数は有界で台がコンパクトなので、 により $x_2$ について積分でき、$\int_\Omega v\,\partial_1\varphi\,dx = -\int_\Omega g_1\varphi\,dx$ です。これは $g_1$ が $v$ の弱微分 $\partial_1v$ であること（）を示しています。$\partial_2$ も同様で、$v$、$g_1$、$g_2$ は有界なので $v \in H^1(\Omega)$ です。 **境界で $0$ の場合.** $\varepsilon > 0$ に対し $v_\varepsilon = (v - \varepsilon)^+ - (-v - \varepsilon)^+$ とおきます。 から $v_\varepsilon \in H^1(\Omega)$ で、$\nabla v_\varepsilon = \mathbf{1}_{\{|v| > \varepsilon\}}\nabla v$（a.e.）です。集合 $\{|v| \ge \varepsilon\}$ は $\bar\Omega$ の閉部分集合で $\partial\Omega$ と交わらないので、$\Omega$ のコンパクト部分集合です。$v_\varepsilon$ の台はそこに含まれるので、 の最後の主張から $v_\varepsilon \in H^1_0(\Omega)$ です。$|v - v_\varepsilon| \le \varepsilon$ なので $\|v - v_\varepsilon\| \to 0$ です。また $\nabla(v - v_\varepsilon) = \mathbf{1}_{\{|v| \le \varepsilon\}}\nabla v$ は、$\varepsilon \to 0$ で各点で $\mathbf{1}_{\{v = 0\}}\nabla v$ に収束し、$|\nabla v|$ で抑えられるので、優収束定理から $\|\nabla(v - v_\varepsilon)\|^2 \to \int_{\{v=0\}}|\nabla v|^2dx$ です。 を $v$ と $-v$ に使うと $\nabla v = \nabla v^+ - \nabla(-v)^+ = \mathbf{1}_{\{v \ne 0\}}\nabla v$（a.e.）なので、この積分は $0$ です。したがって $v_\varepsilon \to v$（$H^1$）で、$H^1_0(\Omega)$ は閉じているので $v \in H^1_0(\Omega)$ です。 **モデル問題.** $f \in L^2(\Omega)$ に対する Poisson 方程式 $-\Delta u = f$、$u|_{\partial\Omega} = 0$ の[弱解](/mathematics/pde/elliptic-weak-solutions#def-weak-solution)は $$ u \in H^1_0(\Omega),\qquad a(u,v) := \langle\nabla u, \nabla v\rangle = \langle f, v\rangle\quad(v \in H^1_0(\Omega)) $$ です（）。$\Omega$ が幅 $d$ の帯に含まれるとき、Poincaré の不等式 $\|v\| \le d\|\nabla v\|$（）から、$\|v\|_V = \|\nabla v\|$ は $V = H^1_0(\Omega)$ のノルムで、$a$ はこのノルムの内積そのものです（$M = \alpha = 1$）。弱解はただ 1 つ存在します（）。 から $V_h \subset V$ なので、有限要素解 $u_h \in V_h$（$a(u_h, v_h) = \langle f, v_h\rangle$）について、前章の適切性（）、Galerkin 直交性（）、Céa の補題（）が成り立ちます。$a$ が内積なので、Céa の補題は $$ \|\nabla(u - u_h)\| = \min_{v_h \in V_h}\|\nabla(u - v_h)\| $$ という等式になります。残る仕事は右辺の評価です。 ## 3. 三角形の上の H² 関数：連続性と 1 次式による近似 節点補間を定義するには、$u$ の点での値が要ります。全空間 $\mathbb{R}^n$ では $s > n/2$ のとき $H^s$ の関数が連続です（）が、ここでは三角形の上で、定数が三角形の大きさにどうよるかまで必要なので、直接証明します。 $S$ を閉じた三角形または $\bar\Omega$ とする。$S$ の近傍で定義された $C^2$ 級関数の $S$ への制限の全体を $C^2(S)$ と書く。$S$ の内部の上の $H^2$ 関数のうち、$C^2(S)$ の関数列の $H^2$ ノルムでの極限になっているものの全体を $\mathcal{H}^2(S)$ と書く。$\partial\Omega$ の上で $0$ である $C^2(\bar\Omega)$ の関数列の $H^2(\Omega)$ ノルムでの極限になっているものの全体を $\mathcal{H}^2_{\mathrm{D}}(\Omega)$ と書く。 三角形や凸多角形のように境界が Lipschitz 連続な領域では、滑らかな関数が $H^2$ で稠密なので $\mathcal{H}^2(S)$ は $H^2$ 全体に一致し、凸多角形では $\mathcal{H}^2_{\mathrm{D}}(\Omega) = H^2(\Omega)\cap H^1_0(\Omega)$ です（Grisvard、参考文献）。この章では、稠密性の一般論を使わずに済むよう、補間誤差の評価を $\mathcal{H}^2$ の元について述べます。正方形の上の Poisson 方程式の弱解が $\mathcal{H}^2_{\mathrm{D}}$ に属することは で直接示します。$\mathcal{H}^2_{\mathrm{D}}(\Omega)$ は $H^2(\Omega)$ の閉部分空間で、 から $H^1_0(\Omega)$ に含まれます。 この節の目標は、三角形の上の $H^2$ 関数 $u$ を「1 次式 + 小さな剰余」に分けることです。1 次式は節点補間でそのまま再現されるので、残るのは剰余の大きさだけで、これを 2 階導関数で抑えます。次の補題の証明は、線分に沿った 1 変数の Taylor の公式を、出発点 $y$ について $T$ の上で平均する、という流れです。変数変換 $z = y + t(x - y)$ は、平均に現れる $y$ を、$x$ に向かって縮めた三角形の上の変数 $z$ に取り替えて、$|D^2u(z)|$ の積分にまとめるためのものです。 $T$ を閉じた三角形（面積 $|T|$、直径 $h_T$）とし、$u \in C^2(T)$ に対し $$ Q_Tu(x) = \frac{1}{|T|}\int_T\bigl(u(y) + \nabla u(y)\cdot(x - y)\bigr)\,dy $$ とおく。$Q_Tu$ は $x$ の 1 次式で、$u$ が 1 次式なら $Q_Tu = u$ である。すべての $x \in T$ について $$ |u(x) - Q_Tu(x)| \le \frac{h_T^2}{2|T|^{1/2}}|u|_{H^2(T)},\qquad|u(x)| \le \frac{1}{|T|^{1/2}}\Bigl(\|u\|_{L^2(T)} + h_T|u|_{H^1(T)} + \frac{h_T^2}{2}|u|_{H^2(T)}\Bigr) $$ が成り立つ。したがって、$u \in \mathcal{H}^2(T)$ は $T$ の上で連続な代表元をもち、$Q_Tu$ を同じ式で定めれば、その代表元について上の 2 つの評価が成り立つ。 $Q_Tu$ が 1 次式であることは定義から、1 次式 $u$ では $u(y) + \nabla u\cdot(x - y) = u(x)$ なので $Q_Tu = u$ です。 **積分表示.** $u \in C^2(T)$、$x, y \in T$ とし、$\phi(t) = u(y + t(x - y))$（$0 \le t \le 1$。$T$ は凸なので線分は $T$ に含まれます）とおきます。部分積分で $\int_0^1(1 - t)\phi''(t)\,dt = -\phi'(0) + \phi(1) - \phi(0)$ なので $$ u(x) = u(y) + \nabla u(y)\cdot(x - y) + \int_0^1(1 - t)\,(x - y)^{\mathsf{T}}D^2u\bigl(y + t(x - y)\bigr)(x - y)\,dt $$ です。$y$ について $T$ の上で平均すると $u(x) = Q_Tu(x) + R(x)$、 $$ |R(x)| \le \frac{1}{|T|}\int_0^1(1 - t)\int_T|x - y|^2\,\bigl|D^2u\bigl(y + t(x - y)\bigr)\bigr|\,dy\,dt $$ です（$|D^2u|$ は成分の 2 乗和の平方根で、$|\xi^{\mathsf{T}}D^2u\,\xi| \le |D^2u||\xi|^2$）。 **変数変換.** $t \in [0, 1)$ を固定し、$z = y + t(x - y) = x + (1 - t)(y - x)$ とおきます。$y \mapsto z$ は $T$ から $T_t = x + (1 - t)(T - x)$ への全単射で、$T$ が凸で $x \in T$ なので $T_t \subset T$、また $dy = (1-t)^{-2}dz$、$x - y = (x - z)/(1 - t)$ です。したがって内側の積分は $(1-t)^{-4}\int_{T_t}|x - z|^2|D^2u(z)|\,dz$ で、 により積分の順序を交換して $$ |R(x)| \le \frac{1}{|T|}\int_T|x - z|^2|D^2u(z)|\Bigl(\int_0^1(1 - t)^{-3}\mathbf{1}_{T_t}(z)\,dt\Bigr)dz $$ です。$z \in T_t$ なら $|z - x| \le (1 - t)h_T$、すなわち $1 - t \ge |z - x|/h_T$ なので、$s = 1 - t$ として括弧の中は $\int_{|z-x|/h_T}^{1}s^{-3}ds \le \frac{h_T^2}{2|z - x|^2}$ 以下です。したがって、Cauchy–Schwarz の不等式から $$ |R(x)| \le \frac{h_T^2}{2|T|}\int_T|D^2u(z)|\,dz \le \frac{h_T^2}{2|T|^{1/2}}|u|_{H^2(T)} $$ です。これが第 1 の評価です。 **第 2 の評価.** $|x - y| \le h_T$ と Cauchy–Schwarz の不等式から $|Q_Tu(x)| \le \frac{1}{|T|}\bigl(|T|^{1/2}\|u\|_{L^2(T)} + h_T|T|^{1/2}|u|_{H^1(T)}\bigr)$ で、$|u(x)| \le |Q_Tu(x)| + |R(x)|$ です。 **$\mathcal{H}^2(T)$ の元.** $u_n \in C^2(T)$、$u_n \to u$（$H^2$）とします。第 2 の評価を $u_n - u_m$ に使うと、$(u_n)$ は $T$ の上で一様に Cauchy 列をなすので、連続関数 $\tilde u$ に一様収束します。$u_n \to u$ は $L^2$ でも成り立つので $\tilde u = u$（a.e.）です。$Q_Tu_n \to Q_Tu$ は一様収束（係数が $u_n$、$\nabla u_n$ の積分で、$L^2$ 収束から収束する）なので、2 つの評価は極限でも成り立ちます。 2 つ目の不等式が、三角形の上の Sobolev の埋め込み $H^2 \subset C$ です。2 次元であることは、$\int s^{-3}ds$ の特異性 $|z - x|^{-2}$ が、分子の $|x - z|^2$ とちょうど打ち消し合うところに現れています（3 次元では $dy = (1-t)^{-3}dz$ となって $|z - x|^{-1}$ が残り、別の評価が要ります）。$Q_Tu$ は、各点 $y$ のまわりの 1 次の Taylor 多項式を $y$ について平均したもので、**平均 Taylor 多項式**と呼ばれます。 勾配の誤差には、平均が $0$ の関数についての Poincaré の不等式を使います。 $T$ を閉じた三角形、$v \in C^1(T)$、$\bar v = \frac{1}{|T|}\int_Tv\,dx$ とすると $$ \|v - \bar v\|_{L^2(T)} \le \sqrt2\,h_T\|\nabla v\|_{L^2(T)}. $$ $x, y \in T$ について $v(x) - v(y) = \int_0^1\nabla v(y + t(x - y))\cdot(x - y)\,dt$ なので、Cauchy–Schwarz の不等式から $|v(x) - v(y)|^2 \le h_T^2\int_0^1|\nabla v(y + t(x - y))|^2dt$ です。$x$、$y$ について積分し、$t$ の範囲を 2 つに分けます。 $t \in [\frac12, 1]$ では、$y$ を固定して $z = y + t(x - y)$ と変換すると、$x \mapsto z$ は $T$ から $y + t(T - y) \subset T$ への全単射で $dx = t^{-2}dz \le 4\,dz$ なので、$\int_T|\nabla v(y + t(x-y))|^2dx \le 4\|\nabla v\|_{L^2(T)}^2$ です。$y$ と $t$ について積分すると $\frac12\cdot|T|\cdot4\|\nabla v\|^2 = 2|T|\|\nabla v\|^2$ 以下です。$t \in [0, \frac12]$ では、$x$ を固定して $z = x + (1-t)(y - x)$ と変換すれば $dy = (1-t)^{-2}dz \le 4\,dz$ で、同じ評価が得られます。合わせて $$ \int_T\int_T|v(x) - v(y)|^2dx\,dy \le 4h_T^2|T|\,\|\nabla v\|_{L^2(T)}^2 $$ です。左辺は、展開すると $2|T|\int_Tv^2 - 2(\int_Tv)^2 = 2|T|\int_T(v - \bar v)^2dx$ なので、$\|v - \bar v\|_{L^2(T)}^2 \le 2h_T^2\|\nabla v\|_{L^2(T)}^2$ です。 $T$ を閉じた三角形、$u \in \mathcal{H}^2(T)$ とする。 の 1 次式 $p = Q_Tu$ は $$ \|u - p\|_{L^\infty(T)} \le \frac{h_T^2}{2|T|^{1/2}}|u|_{H^2(T)},\qquad\|u - p\|_{L^2(T)} \le \frac{h_T^2}{2}|u|_{H^2(T)},\qquad|u - p|_{H^1(T)} \le \sqrt2\,h_T|u|_{H^2(T)} $$ を満たす。 まず $u \in C^2(T)$ とします。第 1 の評価は で、第 2 の評価は $\|w\|_{L^2(T)} \le |T|^{1/2}\|w\|_{L^\infty(T)}$ から従います。$Q_Tu$ の定義を $x$ で微分すると $\nabla(Q_Tu) = \frac{1}{|T|}\int_T\nabla u(y)\,dy$ で、これは $\nabla u$ の平均です。$\partial_1u$、$\partial_2u \in C^1(T)$ のそれぞれに を使うと $$ |u - p|_{H^1(T)}^2 = \sum_{i=1}^{2}\|\partial_iu - \overline{\partial_iu}\|_{L^2(T)}^2 \le 2h_T^2\sum_{i=1}^{2}\|\nabla\partial_iu\|_{L^2(T)}^2 = 2h_T^2|u|_{H^2(T)}^2 $$ です。$u \in \mathcal{H}^2(T)$ のときは、$u_n \in C^2(T)$、$u_n \to u$（$H^2$）について評価を書き、極限をとります（$u_n \to u$ は から一様収束でもあり、$Q_Tu_n \to Q_Tu$ は一様かつ $H^1$ で収束します）。 Bramble–Hilbert の補題は、ふつう「1 次式を法とする商ノルムが $H^2$ の半ノルムと同値である」という形で述べられ、Rellich のコンパクト性定理による背理法で証明されます。その証明は短い代わりに定数の値が分かりません。ここでは 1 次式 $p$ を具体的に作ったので、定数が $h_T$ と $|T|$ だけで書けています。右辺に $\|u\|$ や $|u|_{H^1}$ が現れず、2 階の導関数だけで抑えられていることが要点です。 ## 4. アフィン写像と補間誤差 三角形の上の量を計算するときは、1 つの**参照三角形** $\hat T$（頂点 $(0,0)$、$(1,0)$、$(0,1)$）に移すのが便利です。$\hat T$ の直径は $\hat h = \sqrt2$、内接円の直径は $\hat\rho = 2 - \sqrt2$、面積は $\frac12$ です。頂点 $a_0, a_1, a_2$ の三角形 $T$ に対し、$B = (a_1 - a_0\ \ a_2 - a_0)$（$2 \times 2$ 行列）、$F(\hat x) = B\hat x + a_0$ とすると、$F$ は $\hat T$ を $T$ に写すアフィン写像で、$|T| = \frac12|\det B|$ です。 $F(\hat x) = B\hat x + b$（$B$ は正則）が三角形 $\hat T$ を $T$ に写すとし、$T$ 上の関数 $v$ に $\hat v = v \circ F$ を対応させる。$\|B\|$ を Euclid ノルムから誘導される行列ノルムとする。 1. $v \in C^m(T)$（$m = 0, 1, 2$）について $|\hat v|_{H^m(\hat T)} \le \|B\|^m|\det B|^{-1/2}|v|_{H^m(T)}$、$|v|_{H^m(T)} \le \|B^{-1}\|^m|\det B|^{1/2}|\hat v|_{H^m(\hat T)}$（$|\cdot|_{H^0}$ は $L^2$ ノルム）。 2. $\|B\| \le \dfrac{h_T}{\hat\rho}$、$\|B^{-1}\| \le \dfrac{\hat h}{\rho_T}$。ここで $\hat h$、$\hat\rho$ は $\hat T$ の直径と内接円の直径である。 **1.** 連鎖律から $\hat\nabla\hat v(\hat x) = B^{\mathsf{T}}(\nabla v)(F(\hat x))$、$\hat D^2\hat v(\hat x) = B^{\mathsf{T}}(D^2v)(F(\hat x))B$ です。$|B^{\mathsf{T}}\xi| \le \|B\||\xi|$ で、行列 $N$ の成分の 2 乗和の平方根 $|N|$ について $|B^{\mathsf{T}}NB| \le \|B\|^2|N|$ です（$|B^{\mathsf{T}}N| \le \|B\||N|$ を列ごとに、$|N'B| \le |N'|\|B\|$ を行ごとに使う）。したがって $|\hat D^m\hat v(\hat x)| \le \|B\|^m|D^mv(F(\hat x))|$ で、変数変換 $x = F(\hat x)$、$dx = |\det B|\,d\hat x$ から $$ |\hat v|_{H^m(\hat T)}^2 \le \|B\|^{2m}\int_{\hat T}|D^mv(F(\hat x))|^2d\hat x = \|B\|^{2m}|\det B|^{-1}|v|_{H^m(T)}^2 $$ です。逆向きは、$v = \hat v\circ F^{-1}$、$F^{-1}(x) = B^{-1}x - B^{-1}b$ に同じ議論を使います。 **2.** $\hat T$ の内接円の中心を $\hat c$ とします。$|\xi| = \hat\rho$ である任意のベクトル $\xi$ について、$\hat x = \hat c + \xi/2$ と $\hat y = \hat c - \xi/2$ は内接円の上の点で $\hat T$ に属し、$B\xi = F(\hat x) - F(\hat y)$ は $T$ の 2 点の差なので $|B\xi| \le h_T$ です。したがって $\|B\| = \max_{|\xi| = \hat\rho}|B\xi|/\hat\rho \le h_T/\hat\rho$ です。$F^{-1}$ は $T$ を $\hat T$ に写すので、役割を入れ替えて $\|B^{-1}\| \le \hat h/\rho_T$ です。 1 は、導関数を 1 回とるごとに $\|B\| \approx h_T$ が掛かること（長さのスケール）と、積分で $|\det B| \approx h_T^2$ が掛かること（面積のスケール）を表しています。2 が、三角形の形（$\rho_T$）が評価に入ってくる場所です。逆不等式（§6）ではこの補題を使います。補間誤差は、前節の評価が $T$ の上で直接成り立つので、写像を経由せずに証明できます。 $T$ を閉じた三角形、$u \in \mathcal{H}^2(T)$ とし、$I_Tu = \sum_{i=0}^{2}u(a_i)\lambda_i$（$u$ は連続な代表元）とする。このとき $$ \|u - I_Tu\|_{L^2(T)} \le h_T^2|u|_{H^2(T)},\qquad|u - I_Tu|_{H^1(T)} \le 3\,\frac{h_T^2}{\rho_T}|u|_{H^2(T)} $$ が成り立つ。形状正則な分割（$h_T \le \sigma\rho_T$）では右辺は $3\sigma h_T|u|_{H^2(T)}$ 以下である。 $p = Q_Tu$、$w = u - p$ とおきます。$p$ は 1 次式なので $I_Tp = p$ で、$I_T$ は線形なので $$ u - I_Tu = (u - p) - I_T(u - p) = w - I_Tw $$ です。 から $\|w\|_{L^\infty(T)} \le \frac{h_T^2}{2|T|^{1/2}}|u|_{H^2(T)}$ です。 **$L^2$ ノルム.** $T$ の上で $\lambda_i \ge 0$、$\sum_i\lambda_i = 1$ なので $|I_Tw| \le \|w\|_{L^\infty(T)}$ で、$\|I_Tw\|_{L^2(T)} \le |T|^{1/2}\|w\|_{L^\infty(T)} \le \frac{h_T^2}{2}|u|_{H^2(T)}$ です。 の $\|w\|_{L^2(T)} \le \frac{h_T^2}{2}|u|_{H^2(T)}$ と合わせて $\|u - I_Tu\|_{L^2(T)} \le h_T^2|u|_{H^2(T)}$ です。 **$H^1$ 半ノルム.** $\lambda_i$ は、頂点 $a_i$ で $1$、向かい合う辺の上で $0$ の 1 次式なので、$|\nabla\lambda_i| = 1/H_i$（$H_i$ は $a_i$ から向かい合う辺への高さ）です。内接円は $T$ に含まれ、その直径は高さを超えないので $H_i \ge \rho_T$、したがって $|\nabla\lambda_i| \le 1/\rho_T$ です。$\nabla(I_Tw) = \sum_iw(a_i)\nabla\lambda_i$ から $|\nabla(I_Tw)| \le \frac{3}{\rho_T}\|w\|_{L^\infty(T)}$ で、 $$ |I_Tw|_{H^1(T)} \le |T|^{1/2}\cdot\frac{3}{\rho_T}\cdot\frac{h_T^2}{2|T|^{1/2}}|u|_{H^2(T)} = \frac{3h_T^2}{2\rho_T}|u|_{H^2(T)} $$ です。 の $|w|_{H^1(T)} \le \sqrt2\,h_T|u|_{H^2(T)}$ と、$\rho_T \le h_T$ から $h_T \le h_T^2/\rho_T$ であることを合わせて $$ |u - I_Tu|_{H^1(T)} \le \Bigl(\sqrt2 + \frac32\Bigr)\frac{h_T^2}{\rho_T}|u|_{H^2(T)} \le 3\,\frac{h_T^2}{\rho_T}|u|_{H^2(T)} $$ です。 $L^2$ の評価には三角形の形が現れず、$H^1$ の評価にだけ $1/\rho_T$ が現れます。原因は補間した関数の[勾配](/mathematics/calculus/vector-analysis#def-grad-div-rot) $\nabla\lambda_i$ で、高さの低い三角形では、頂点での値のわずかな違いが大きな勾配になります。 $0 < \varepsilon \le 1$ とし、頂点 $(-1, 0)$、$(1, 0)$、$(0, \varepsilon)$ の三角形 $T_\varepsilon$ と $u(x,y) = x^2$ を考えます。頂点での値は $1, 1, 0$ なので、補間は $I_Tu = 1 - y/\varepsilon$ です。$\partial_yu = 0$、$\partial_y(I_Tu) = -1/\varepsilon$ なので $$ |u - I_Tu|_{H^1(T_\varepsilon)}^2 \ge \int_{T_\varepsilon}\frac{1}{\varepsilon^2}\,dx = \frac{|T_\varepsilon|}{\varepsilon^2} = \frac{1}{\varepsilon},\qquad|u|_{H^2(T_\varepsilon)}^2 = \int_{T_\varepsilon}2^2dx = 4\varepsilon $$ です。$h_T = 2$ は一定なのに、比 $|u - I_Tu|_{H^1}/|u|_{H^2} \ge \frac{1}{2\varepsilon}$ は $\varepsilon \to 0$ で発散します。内接円の直径は $\rho_T \le \varepsilon$（高さ以下）なので、これは の $h_T^2/\rho_T \ge 4/\varepsilon$ と矛盾せず、$1/\rho_T$ の因子が証明の都合でないことを示しています。この三角形は $(0,\varepsilon)$ の角が $\pi$ に近づきます。 つぶれた三角形がすべて悪いわけではありません。頂点 $(0,0)$、$(1,0)$、$(0,\varepsilon)$ の細長い直角三角形では、$H^1$ の補間誤差は $\varepsilon$ によらず $O(h_T)$ のままです。補間誤差が $O(h_T)$ であるために本質的なのは、最小の角が $0$ から離れていること（形状正則性）ではなく、**最大の角が $\pi$ から離れていること**です（Babuška–Aziz、参考文献）。境界層のように一方向にだけ急に変わる解には、この性質を使って細長い要素をわざと使います。この章では扱いやすい形状正則性を仮定します。 ## 5. Poisson 方程式の誤差評価 $\Omega$ を多角形領域、$\mathcal{T}_h$ を $h_T/\rho_T \le \sigma$ を満たす三角形分割とする。任意の $w \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$ は $\bar\Omega$ で連続で $\partial\Omega$ の上で $0$ の代表元をもち、$I_hw \in V_h$ が定まって $\|\nabla(w - I_hw)\| \le 3\sigma h|w|_{H^2(\Omega)}$ である。したがって、$f \in L^2(\Omega)$ に対する Poisson 方程式の弱解 $u$ が $\mathcal{H}^2_{\mathrm{D}}(\Omega)$ に属するならば、区分 1 次の有限要素解 $u_h \in V_h$ は $$ \|\nabla(u - u_h)\| \le \|\nabla(u - I_hu)\| \le 3\sigma h|u|_{H^2(\Omega)} $$ を満たす。 $u \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$ を任意にとり（前半の $w$ を $u$ と書きます）、$u_n \in C^2(\bar\Omega)$、$u_n|_{\partial\Omega} = 0$、$u_n \to u$（$H^2(\Omega)$）とします。各 $T \in \mathcal{T}_h$ の上で $u|_T \in \mathcal{H}^2(T)$ で、 から $u_n$ は $T$ の上で一様収束します。三角形は有限個なので、$u_n$ は $\bar\Omega$ の上で一様収束し、極限 $\tilde u$ は $\bar\Omega$ で連続、$\partial\Omega$ の上で $0$、$\tilde u = u$（a.e.）です。したがって $I_hu$ が定まって $I_hu \in V_h$ で、$T$ の上で $I_hu = I_Tu$ です。 から $I_hu$ の弱微分は三角形ごとの勾配なので、 を三角形ごとに使って足すと $$ \|\nabla(u - I_hu)\|^2 = \sum_{T}|u - I_Tu|_{H^1(T)}^2 \le \sum_T(3\sigma h_T)^2|u|_{H^2(T)}^2 \le (3\sigma h)^2|u|_{H^2(\Omega)}^2 $$ です。これが前半です。$u$ が弱解のときは、$a$ が $V$ の内積なので、Céa の補題（ の対称な場合）から $\|\nabla(u - u_h)\| \le \|\nabla(u - I_hu)\|$ で、後半を得ます。 仮定「$u \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$」がいつ成り立つかが問題です。正方形では、固有関数展開で直接示せます。 $\Omega = (0,1)^2$、$f \in L^2(\Omega)$ とする。$-\Delta u = f$ の弱解 $u \in H^1_0(\Omega)$ は $\mathcal{H}^2_{\mathrm{D}}(\Omega)$ に属し、 $$ |u|_{H^2(\Omega)} = \|D^2u\| = \|f\| $$ を満たす。 $k = (k_1, k_2) \in \mathbb{N}^2$ について $s_k(x) = 2\sin(k_1\pi x_1)\sin(k_2\pi x_2)$、$\mu_k = \pi^2(k_1^2 + k_2^2)$ とおきます。 から $\{s_k\}$ は $L^2(\Omega)$ の正規直交基底で、余弦の積 $e_k(x) = 2\cos(k_1\pi x_1)\cos(k_2\pi x_2)$（$k_1, k_2 \ge 1$）は正規直交系です。$s_k$ は $\bar\Omega$ の近傍で滑らかで $\partial\Omega$ の上で $0$、$-\Delta s_k = \mu_ks_k$ です。 **有限和での等式.** 有限和 $w = \sum_kc_ks_k$ について、$\partial_1^2w = -\sum_kc_k\pi^2k_1^2s_k$、$\partial_2^2w = -\sum_kc_k\pi^2k_2^2s_k$、$\partial_1\partial_2w = \sum_kc_k\pi^2k_1k_2e_k$ なので、正規直交性から $$ |w|_{H^2(\Omega)}^2 = \|\partial_1^2w\|^2 + 2\|\partial_1\partial_2w\|^2 + \|\partial_2^2w\|^2 = \sum_kc_k^2\pi^4(k_1^4 + 2k_1^2k_2^2 + k_2^4) = \sum_kc_k^2\mu_k^2 = \|\Delta w\|^2 $$ です。同様に、 から $\|\nabla w\|^2 = \sum_kc_k^2\mu_k$、$\|w\|^2 = \sum_kc_k^2$ で、$\mu_k \ge 2\pi^2 > 1$ なのでどちらも $\|\Delta w\|^2$ 以下です。 **極限.** $f = \sum_kf_ks_k$（$L^2$ 収束）と展開し、$f_N = \sum_{|k| \le N}f_ks_k$、$u_N = \sum_{|k| \le N}\frac{f_k}{\mu_k}s_k$ とおくと、$-\Delta u_N = f_N$ です。上の等式を $u_N - u_M$ に使うと $\|u_N - u_M\|^2 + \|\nabla(u_N - u_M)\|^2 + |u_N - u_M|_{H^2}^2 \le 3\|f_N - f_M\|^2 \to 0$ で、$(u_N)$ は $H^2(\Omega)$ の Cauchy 列です。$H^2(\Omega)$ は完備（）なので $u_N \to \tilde u$（$H^2$）となる $\tilde u$ があり、$u_N \in C^2(\bar\Omega)$ は $\partial\Omega$ の上で $0$ なので $\tilde u \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$、$|\tilde u|_{H^2} = \lim_N\|f_N\| = \|f\|$ です。 **$\tilde u$ は弱解.** $u_N \in H^1_0(\Omega)$（）で、$v \in C_c^\infty(\Omega)$ について、弱微分の定義から $\langle\nabla u_N, \nabla v\rangle = -\langle\Delta u_N, v\rangle = \langle f_N, v\rangle$ です。$N \to \infty$ として $\langle\nabla\tilde u, \nabla v\rangle = \langle f, v\rangle$ で、$C_c^\infty(\Omega)$ は $H^1_0(\Omega)$ で稠密なのですべての $v \in H^1_0(\Omega)$ で成り立ちます。$H^1_0(\Omega)$ は閉じているので $\tilde u \in H^1_0(\Omega)$ で、弱解の一意性（）から $\tilde u = u$ です。 辺の長さが異なる長方形でも、各方向の正弦の積が[正規直交基底](/mathematics/functional-analysis/hilbert-spaces#def-orthonormal-system)になることを使えば、同じ証明が通ります（下の ）。滑らかな境界をもつ凸領域での $\|D^2u\| \le \|\Delta u\|$ は で証明されていて、一般の凸多角形にも同じ評価が成り立ちます（同じ章の注意で述べられている Grisvard の定理）。そこでこの章では、$L^2$ の評価を次の仮定のもとで述べます。 **仮定 (R).** 定数 $C_R$ があって、すべての $g \in L^2(\Omega)$ について、$-\Delta z = g$ の弱解 $z \in H^1_0(\Omega)$ は $\mathcal{H}^2_{\mathrm{D}}(\Omega)$ に属し、$|z|_{H^2(\Omega)} \le C_R\|g\|$ を満たす。 から、正方形では (R) が $C_R = 1$ で成り立ちます。辺の長さが異なる長方形でも同じです。 $a, b > 0$、$\Omega = (0,a)\times(0,b)$、$f \in L^2(\Omega)$ とする。$-\Delta u = f$ の弱解 $u \in H^1_0(\Omega)$ は $\mathcal{H}^2_{\mathrm{D}}(\Omega)$ に属し、$|u|_{H^2(\Omega)} = \|f\|$ を満たす。したがって、この長方形で仮定 (R) が $C_R = 1$ で成り立つ。 の証明を、方向ごとに周期の違う正弦で書き直します。$\ell > 0$、$m \ge 1$ について $\sigma^{(\ell)}_m(t) = \sqrt{2/\ell}\,\sin(m\pi t/\ell)$、$\gamma^{(\ell)}_m(t) = \sqrt{2/\ell}\,\cos(m\pi t/\ell)$ とおくと、 の前半（区間 $(0,\ell)$ の場合）から、$\{\sigma^{(\ell)}_m\}_{m\ge1}$ と $\{\gamma^{(\ell)}_m\}_{m\ge1}$ はそれぞれ $L^2(0,\ell)$ の正規直交系です。$k = (k_1,k_2) \in \mathbb{N}^2$ について $$ s_k = \sigma^{(a)}_{k_1}(x_1)\sigma^{(b)}_{k_2}(x_2),\quad e_k = \gamma^{(a)}_{k_1}(x_1)\gamma^{(b)}_{k_2}(x_2),\quad p_k = \gamma^{(a)}_{k_1}(x_1)\sigma^{(b)}_{k_2}(x_2),\quad q_k = \sigma^{(a)}_{k_1}(x_1)\gamma^{(b)}_{k_2}(x_2),\qquad\mu_k = \pi^2\Bigl(\frac{k_1^2}{a^2} + \frac{k_2^2}{b^2}\Bigr) $$ とおきます。 **基底.** 積の形 $g_1(x_1)g_2(x_2)$ の連続関数どうしの $L^2(\Omega)$ の内積は、 により $L^2(0,a)$ の内積と $L^2(0,b)$ の内積の積になるので、$\{s_k\}$、$\{e_k\}$、$\{p_k\}$、$\{q_k\}$ はそれぞれ $L^2(\Omega)$ の正規直交系です。$\{s_k\}$ が基底であることを示します。$g \in L^2((0,1)^2)$ に対し $(Ug)(x) = (ab)^{-1/2}g(x_1/a, x_2/b)$ とおくと、変数変換 $y = (x_1/a, x_2/b)$、$dx = ab\,dy$ から $\|Ug\|_{L^2(\Omega)} = \|g\|_{L^2((0,1)^2)}$ で、$U$ は逆写像 $(U^{-1}w)(y) = (ab)^{1/2}w(ay_1, by_2)$ をもつので、$L^2((0,1)^2)$ から $L^2(\Omega)$ の上へのユニタリ写像です。単位正方形の関数 $s^{\square}_k(y) = 2\sin(k_1\pi y_1)\sin(k_2\pi y_2)$ について $Us^{\square}_k = s_k$ です。$w \in L^2(\Omega)$ がすべての $s_k$ と直交すれば、$U^{-1}w$ はすべての $s^{\square}_k$ と直交し、$\{s^{\square}_k\}$ は $L^2((0,1)^2)$ の正規直交基底（）なので $U^{-1}w = 0$、すなわち $w = 0$ です。したがって $\{s_k\}$ は $L^2(\Omega)$ の正規直交基底です。 **有限和での等式.** 有限和 $w = \sum_k\beta_ks_k$ について $$ \partial_1w = \sum_k\beta_k\frac{\pi k_1}{a}p_k,\quad\partial_2w = \sum_k\beta_k\frac{\pi k_2}{b}q_k,\quad\partial_1^2w = -\sum_k\beta_k\frac{\pi^2k_1^2}{a^2}s_k,\quad\partial_2^2w = -\sum_k\beta_k\frac{\pi^2k_2^2}{b^2}s_k,\quad\partial_1\partial_2w = \sum_k\beta_k\frac{\pi^2k_1k_2}{ab}e_k $$ なので、正規直交性から $\|w\|^2 = \sum_k\beta_k^2$、$\|\nabla w\|^2 = \sum_k\beta_k^2\mu_k$、 $$ |w|_{H^2(\Omega)}^2 = \|\partial_1^2w\|^2 + 2\|\partial_1\partial_2w\|^2 + \|\partial_2^2w\|^2 = \sum_k\beta_k^2\pi^4\Bigl(\frac{k_1^4}{a^4} + \frac{2k_1^2k_2^2}{a^2b^2} + \frac{k_2^4}{b^4}\Bigr) = \sum_k\beta_k^2\mu_k^2 = \|\Delta w\|^2 $$ です。$\mu_{*} = \pi^2(a^{-2} + b^{-2})$ とおくと $\mu_k \ge \mu_{*}$ なので、$\|w\|^2 \le \mu_{*}^{-2}\|\Delta w\|^2$、$\|\nabla w\|^2 \le \mu_{*}^{-1}\|\Delta w\|^2$ です。 **極限.** $f = \sum_kf_ks_k$（$L^2$ 収束）と展開し、$f_N = \sum_{|k|\le N}f_ks_k$、$u_N = \sum_{|k|\le N}\frac{f_k}{\mu_k}s_k$ とおくと $-\Delta u_N = f_N$ です。上の評価を $u_N - u_M$ に使うと $$ \|u_N - u_M\|^2 + \|\nabla(u_N - u_M)\|^2 + |u_N - u_M|_{H^2}^2 \le \bigl(\mu_{*}^{-2} + \mu_{*}^{-1} + 1\bigr)\|f_N - f_M\|^2 \to 0 $$ で、$(u_N)$ は $H^2(\Omega)$ の Cauchy 列です。 から $u_N \to \tilde u$（$H^2$）となる $\tilde u$ があります。$u_N$ は $\bar\Omega$ の近傍で滑らかで $\partial\Omega$ の上で $0$ なので、$\tilde u \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$ で、$|\tilde u|_{H^2} = \lim_N|u_N|_{H^2} = \lim_N\|f_N\| = \|f\|$ です。 **$\tilde u$ は弱解.** $\bar\Omega$ を対角線で 2 つの三角形に分けた分割に を使うと $u_N \in H^1_0(\Omega)$ です。$v \in C_c^\infty(\Omega)$ について、$\partial_iu_N$ の弱微分が $\partial_i^2u_N$ であること（滑らかな関数なので部分積分から従います）から $\langle\nabla u_N, \nabla v\rangle = -\langle\Delta u_N, v\rangle = \langle f_N, v\rangle$ です。$N \to \infty$ として $\langle\nabla\tilde u, \nabla v\rangle = \langle f, v\rangle$ で、$C_c^\infty(\Omega)$ は $H^1_0(\Omega)$ で稠密（ の $H^1_0$ の定義）なので、すべての $v \in H^1_0(\Omega)$ で成り立ちます。$H^1_0(\Omega)$ は $H^1(\Omega)$ の閉部分空間なので $\tilde u \in H^1_0(\Omega)$ で、弱解の一意性（）から $\tilde u = u$ です。 平行移動しても弱解・$H^2$ 半ノルム・$L^2$ ノルムは変わらないので、辺が座標軸に平行な任意の長方形 $x_0 + (0,a)\times(0,b)$ でも同じ結論が成り立ちます。 ここまでの $H^1$ の誤差は $O(h)$ でした。$L^2$ の誤差 $\|e\|$ では、さらに $h$ が 1 つ分よくなることを示します。そのために $\|e\|^2 = \langle e, e\rangle$ に注目し、$e$ 自身を右辺にもつ補助の方程式 $-\Delta z = e$（双対問題）を立てます。すると $\|e\|^2$ が「$z$ の補間誤差」と「$\|\nabla e\|$」の積に書き直せ、$z$ の補間誤差のほうにもう 1 つ $h$ が出ます。 $\Omega$ を多角形領域、$\mathcal{T}_h$ を $h_T/\rho_T \le \sigma$ を満たす三角形分割とし、仮定 (R) が成り立つとする。$w \in H^1_0(\Omega)$ とその Ritz 射影 $R_hw \in V_h$（$\langle\nabla(R_hw - w), \nabla v_h\rangle = 0$、$v_h \in V_h$）について $$ \|w - R_hw\| \le 3\sigma C_Rh\,\|\nabla(w - R_hw)\| $$ が成り立つ。さらに $w \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$ ならば $\|w - R_hw\| \le C_Ah^2|w|_{H^2(\Omega)}$、$C_A = 9\sigma^2C_R$ である。特に、Poisson 方程式の弱解 $u$（(R) から $\mathcal{H}^2_{\mathrm{D}}(\Omega)$ に属する）とその有限要素解 $u_h = R_hu$ について $$ \|u - u_h\| \le 9\sigma^2C_Rh^2|u|_{H^2(\Omega)} \le 9\sigma^2C_R^2h^2\|f\|. $$ $e = w - R_hw \in H^1_0(\Omega)$ とし、$-\Delta z = e$ の弱解を $z$ とします。(R) から $z \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$、$|z|_{H^2} \le C_R\|e\|$ です。弱解の式で試験関数を $e$ とし、Ritz 射影の定義 $\langle\nabla e, \nabla I_hz\rangle = 0$（$I_hz \in V_h$。 の証明のとおり定まる）を使うと $$ \|e\|^2 = \langle\nabla z, \nabla e\rangle = \langle\nabla(z - I_hz), \nabla e\rangle \le \|\nabla(z - I_hz)\|\,\|\nabla e\| \le 3\sigma h\,C_R\|e\|\,\|\nabla e\| $$ です（ の前半を $z$ に使いました）。$\|e\|$ で割って前半です。$w \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$ なら、Ritz 射影は $\|\nabla\cdot\|$ についての最良近似なので $\|\nabla e\| \le \|\nabla(w - I_hw)\| \le 3\sigma h|w|_{H^2}$ で、後半を得ます。$u_h = R_hu$ は Galerkin 直交性（）です。 内角が $\pi$ を超える角（凹角）をもつ多角形では、(R) は成り立ちません。内角 $\omega > \pi$ の角を原点とする極座標で、$u = r^{\pi/\omega}\sin(\pi\theta/\omega)$ は調和で、角をはさむ 2 辺の上で $0$ ですが、2 階導関数は $r^{\pi/\omega - 2}$ の大きさで、$\pi/\omega < 1$ なので 2 乗可積分ではありません。L 字形の領域（$\omega = 3\pi/2$）では $u = r^{2/3}\sin(2\theta/3)$ です。このとき、準一様な分割での有限要素解の誤差は $H^1$ で $O(h^{2/3})$、$L^2$ で $O(h^{4/3})$ になります。証明には重みつきの Sobolev 空間か、指数が整数でない Sobolev 空間が要るので、ここでは §8 で数値的に確かめるにとどめます。Grisvard（参考文献）が角のある領域の正則性を体系的に扱っています。角に向かって分割を細かくすれば、前章の 1 次元の例と同じく次数を回復できます。 ## 6. 剛性行列：余接の公式、離散最大値原理、逆不等式 すべての節点（境界の上のものも含む）について $\tilde K_{ij} = \langle\nabla\varphi_j, \nabla\varphi_i\rangle$ とおき、内部の節点に対応する部分を $K$（剛性行列）、$G_{ij} = \langle\varphi_j,\varphi_i\rangle$（内部の節点）を質量行列とします。 次の命題の 3 の証明では、境界の最大値 $m$ を超える内部の節点があると仮定して、超えた分だけを集めたベクトル $\boldsymbol{p}$ をつくり、二次形式 $\boldsymbol{p}^{\mathsf{T}}K\boldsymbol{p}$ を調べます。2 で示す「非対角成分が $0$ 以下」という符号のおかげで、ほかの項がすべて $0$ 以上になり、$K$ が正定値であることから $\boldsymbol{p} = \boldsymbol{0}$ が従います。 1. 頂点 $a_i, a_j, a_k$ の三角形 $T$ で、$a_k$ の角を $\theta_k$ とすると $\displaystyle\int_T\nabla\lambda_i\cdot\nabla\lambda_j\,dx = -\frac12\cot\theta_k$。したがって、辺で結ばれた相異なる節点 $i$、$j$ について $\tilde K_{ij} = -\frac12\sum\cot\theta$（和はその辺を含む三角形についてとり、$\theta$ はその三角形で辺に向かい合う角）で、辺で結ばれていなければ $\tilde K_{ij} = 0$。 2. 内部の節点を端点とするすべての辺について、その辺に向かい合う 2 つの角の和が $\pi$ 以下ならば、内部の節点 $i$ と任意の節点 $j \ne i$ について $\tilde K_{ij} \le 0$ である。このとき $K$ は非対角成分が $0$ 以下の正定値対称行列（Stieltjes 行列。）である。 3. 2 の仮定のもとで、$u_h \in X_h$ がすべての内部の節点 $i$ について $\langle\nabla u_h, \nabla\varphi_i\rangle \le 0$ を満たすならば、$\max_{\bar\Omega}u_h = \max_{\partial\Omega}u_h$ である（**離散最大値原理**）。 **1.** $a_i$ に向かい合う辺を $E_i$、その長さを $|E_i|$、$a_i$ からの高さを $H_i$ とします。$\nabla\lambda_i$ は $E_i$ に垂直で三角形の内側（$a_i$ の側）を向き、大きさは $1/H_i = |E_i|/(2|T|)$ です。$E_i$ と $E_j$ は頂点 $a_k$ で角 $\theta_k$ をなして交わり、それぞれの内向きの法線のなす角は $\pi - \theta_k$ です。したがって $$ \nabla\lambda_i\cdot\nabla\lambda_j = \frac{|E_i||E_j|}{4|T|^2}\cos(\pi - \theta_k) = -\frac{|E_i||E_j|\cos\theta_k}{4|T|^2} $$ で、$|T| = \frac12|E_i||E_j|\sin\theta_k$ を使うと $\int_T\nabla\lambda_i\cdot\nabla\lambda_j\,dx = -\frac{|E_i||E_j|\cos\theta_k}{4|T|} = -\frac12\cot\theta_k$ です。$\varphi_i$ と $\varphi_j$ の台が重なるのは、両方の節点を頂点にもつ三角形の上だけで、そこでは $\varphi_i = \lambda_i$、$\varphi_j = \lambda_j$ なので、後半が従います。 **2.** 内部の節点 $i$ を端点とする辺は $\partial\Omega$ の上にないので、ちょうど 2 つの三角形に含まれます。向かい合う角を $\theta$、$\theta'$（$0 < \theta, \theta' < \pi$）とすると $$ \tilde K_{ij} = -\frac12(\cot\theta + \cot\theta') = -\frac{\sin(\theta + \theta')}{2\sin\theta\sin\theta'} $$ で、$\theta + \theta' \le \pi$ なら $\tilde K_{ij} \le 0$ です。$K$ は対称で正定値（）です。 **3.** $\sum_j\varphi_j = 1$（どの節点でも値が $1$ の $X_h$ の元）なので、$\sum_j\tilde K_{ij} = \langle\nabla1, \nabla\varphi_i\rangle = 0$ です。$m = \max_{\partial\Omega}u_h$ とします。$u_h$ は境界の辺の上で 1 次式なので、$m$ は境界の節点での値の最大です。$d_j = u_h(z_j) - m$ とおくと、境界の節点では $d_j \le 0$ で、行の和が $0$ であることから、内部の節点 $i$ について $$ \sum_j\tilde K_{ij}d_j = \sum_j\tilde K_{ij}u_h(z_j) = \langle\nabla u_h, \nabla\varphi_i\rangle \le 0 $$ です。内部の節点について $p_i = \max(d_i, 0)$ とおき、上の式に $p_i \ge 0$ を掛けて足すと $$ 0 \ge \sum_{i,j\,:\,\text{内部}}p_iK_{ij}p_j + \sum_{i,j\,:\,\text{内部}}p_iK_{ij}(d_j - p_j) + \sum_{i\,:\,\text{内部}}\ \sum_{j\,:\,\text{境界}}p_i\tilde K_{ij}d_j $$ です。第 2 の和では、$d_j - p_j = \min(d_j, 0) \le 0$ で、$i = j$ の項は $p_i\min(d_i, 0) = 0$、$i \ne j$ の項は $p_i \ge 0$、$K_{ij} \le 0$、$d_j - p_j \le 0$ の積なので $0$ 以上です。第 3 の和も $p_i \ge 0$、$\tilde K_{ij} \le 0$、$d_j \le 0$ の積の和で $0$ 以上です。したがって $\boldsymbol{p}^{\mathsf{T}}K\boldsymbol{p} \le 0$ で、$K$ が正定値なので $\boldsymbol{p} = \boldsymbol{0}$、すなわち内部の節点でも $u_h(z_i) \le m$ です。$u_h$ は各三角形の上で 1 次式なので、最大値は節点でとられ、$\max_{\bar\Omega}u_h = m$ です。 3 の仮定 $\langle\nabla u_h, \nabla\varphi_i\rangle \le 0$ は、弱い意味での $-\Delta u \le 0$（ の仮定）を試験関数 $\varphi_i \ge 0$ に限ったものです。前々章の差分法の離散最大値原理（）の有限要素版にあたります。2 の条件「向かい合う 2 角の和が $\pi$ 以下」は、Delaunay 三角形分割（どの三角形の外接円の内部にも他の節点がない分割）の特徴づけになっています。すべての角が $\pi/2$ 以下（鈍角がない）なら当然満たされます。鈍角があると破れうることを §8 で見ます。 単位正方形を一辺 $h = 1/n$ の正方形に分け、各正方形を同じ向きの対角線で 2 つの直角二等辺三角形に分けます。内部の節点 $z_i$ から出る辺は 6 本です。座標軸に平行な 4 本は、両側の三角形の直角をはさむ辺で、向かい合う角はどちらも $\pi/4$ なので $\tilde K_{ij} = -\frac12(1 + 1) = -1$ です。対角線の 2 本は両側の三角形の斜辺で、向かい合う角はどちらも $\pi/2$ なので $\tilde K_{ij} = 0$ です。行の和が $0$ なので $\tilde K_{ii} = 4$ です。したがって $K\boldsymbol{c} = \boldsymbol{b}$ の第 $i$ 式は $$ 4c_i - c_{\text{東}} - c_{\text{西}} - c_{\text{南}} - c_{\text{北}} = b_i = \int_\Omega f\varphi_i\,dx $$ で、$\int_\Omega\varphi_i\,dx = 6\cdot\frac13\cdot\frac{h^2}{2} = h^2$（$z_i$ を頂点とする 6 つの三角形の上で、1 次式 $\lambda_i$ の積分は面積の $\frac13$）なので、右辺は $h^2$ と $f$ の重みつき平均の積です。$h^2$ で割れば 5 点差分 $-\Delta_hc = \bar f$ になります。1 次元の場合（）と同じく、規則的な分割では有限要素法と差分法は右辺の扱いだけが違います。 1. 三角形 $T$ の上の 1 次式 $v$（頂点での値 $\alpha_0, \alpha_1, \alpha_2$）について、$\|v\|_{L^2(T)}^2 = \dfrac{|T|}{12}\Bigl(\sum_i\alpha_i^2 + \bigl(\sum_i\alpha_i\bigr)^2\Bigr)$ で、$\dfrac{|T|}{12}\sum_i\alpha_i^2 \le \|v\|_{L^2(T)}^2 \le \dfrac{|T|}{3}\sum_i\alpha_i^2$、$\|\nabla v\|_{L^2(T)} \le \dfrac{\sqrt{96}}{\rho_T}\|v\|_{L^2(T)}$。 2. 形状正則（定数 $\sigma$）で準一様（定数 $\gamma$）な分割では、すべての $v_h \in X_h$ について $\|\nabla v_h\| \le \dfrac{C_{\mathrm{inv}}}{h}\|v_h\|$、$C_{\mathrm{inv}} = \sqrt{96}\,\gamma\sigma$。 3. 2 の仮定のもとで、$\Omega$ が幅 $d$ の帯に含まれるとすると、$\kappa_2(G) \le \dfrac{8n_\sigma\gamma^2\sigma^2}{\pi}$、$\kappa_2(K) \le \dfrac{8n_\sigma\gamma^2\sigma^2}{\pi}\,C_{\mathrm{inv}}^2d^2\,h^{-2}$。ここで $n_\sigma$ は 1 つの節点を共有する三角形の個数の上界である。 4. 2 の仮定のもとで、陽的 Euler 法 $\langle U^{n+1} - U^n, v_h\rangle + \Delta t\langle\nabla U^n, \nabla v_h\rangle = 0$（$v_h \in V_h$）の解は、$\Delta t \le 2h^2/C_{\mathrm{inv}}^2$ ならば $\|U^{n+1}\| \le \|U^n\|$ を満たす。 **1.** まず $\int_T\lambda_i\lambda_j\,dx = \frac{|T|}{12}(1 + \delta_{ij})$ を示します。$T$ を参照三角形 $\hat T$ に写すアフィン写像で、$\lambda_1$、$\lambda_2$ は座標 $\hat x_1$、$\hat x_2$ になり、積分は $|T|/|\hat T| = 2|T|$ 倍されます。$\int_{\hat T}\hat x_1^2d\hat x = \int_0^1s^2(1 - s)\,ds = \frac{1}{12}$、$\int_{\hat T}\hat x_1\hat x_2\,d\hat x = \int_0^1s\,\frac{(1-s)^2}{2}ds = \frac{1}{24}$ なので、$\int_T\lambda_1^2 = \frac{|T|}{6}$、$\int_T\lambda_1\lambda_2 = \frac{|T|}{12}$ です。頂点の番号を付け替えれば他の組も同じです。$v = \sum_i\alpha_i\lambda_i$ から $$ \|v\|_{L^2(T)}^2 = \frac{|T|}{12}\Bigl(2\sum_i\alpha_i^2 + 2\sum_{i から $$ |v|_{H^1(T)} \le \|B^{-1}\||\det B|^{1/2}|\hat v|_{H^1(\hat T)} \le \|B^{-1}\||\det B|^{1/2}\sqrt{48}\,\|\hat v\|_{L^2(\hat T)} = \sqrt{48}\,\|B^{-1}\|\,\|v\|_{L^2(T)} \le \frac{\sqrt{48}\cdot\sqrt2}{\rho_T}\|v\|_{L^2(T)} $$ です（$\|\hat v\|_{L^2(\hat T)} = |\det B|^{-1/2}\|v\|_{L^2(T)}$、$\|B^{-1}\| \le \hat h/\rho_T = \sqrt2/\rho_T$）。 **2.** $\rho_T \ge h_T/\sigma \ge h/(\gamma\sigma)$ なので、1 から $\|\nabla v_h\|_{L^2(T)} \le \frac{\sqrt{96}\gamma\sigma}{h}\|v_h\|_{L^2(T)}$ で、2 乗して三角形について足します（ から弱微分は三角形ごとの勾配です）。 **3.** $v_h = \sum_ic_i\varphi_i \in V_h$ とします。1 を三角形について足すと、内部の節点 $i$ の値 $c_i^2$ は、それを頂点にもつ三角形（1 個以上、$n_\sigma$ 個以下）から現れるので $$ \frac{A_{\min}}{12}|\boldsymbol{c}|^2 \le \|v_h\|^2 = \boldsymbol{c}^{\mathsf{T}}G\boldsymbol{c} \le \frac{n_\sigma A_{\max}}{3}|\boldsymbol{c}|^2 $$ です（$A_{\min}$、$A_{\max}$ は三角形の面積の最小と最大）。三角形は直径 $\rho_T \ge h/(\gamma\sigma)$ の円を含むので $A_{\min} \ge \frac{\pi}{4}\bigl(\frac{h}{\gamma\sigma}\bigr)^2$、底辺も高さも $h_T \le h$ 以下なので $A_{\max} \le \frac{h^2}{2}$ です。前章の の証明の冒頭の議論（Rayleigh 商の上下の評価から条件数を抑える）から $\kappa_2(G) \le \frac{4n_\sigma A_{\max}}{A_{\min}} \le \frac{8n_\sigma\gamma^2\sigma^2}{\pi}$ です。剛性行列は、Poincaré の不等式（）と 2 から $$ \frac{1}{d^2}\|v_h\|^2 \le \boldsymbol{c}^{\mathsf{T}}K\boldsymbol{c} = \|\nabla v_h\|^2 \le \frac{C_{\mathrm{inv}}^2}{h^2}\|v_h\|^2 $$ なので、質量行列の評価と合わせて $\kappa_2(K) \le \frac{C_{\mathrm{inv}}^2d^2}{h^2}\cdot\frac{4n_\sigma A_{\max}}{A_{\min}}$ です。 **4.** $A_h : V_h \to V_h$ を $\langle A_hw, v_h\rangle = \langle\nabla w, \nabla v_h\rangle$（$v_h \in V_h$）で定めると（有限次元の内積空間なので一意に定まります）、スキームは $U^{n+1} = U^n - \Delta tA_hU^n$ です。$v_h = A_hw$ とすると、2 から $\|A_hw\|^2 = \langle\nabla w, \nabla A_hw\rangle \le \|\nabla w\|\cdot\frac{C_{\mathrm{inv}}}{h}\|A_hw\|$、すなわち $\|A_hw\| \le \frac{C_{\mathrm{inv}}}{h}\|\nabla w\|$ です。したがって $$ \|U^{n+1}\|^2 = \|U^n\|^2 - 2\Delta t\|\nabla U^n\|^2 + \Delta t^2\|A_hU^n\|^2 \le \|U^n\|^2 - \Delta t\Bigl(2 - \frac{C_{\mathrm{inv}}^2\Delta t}{h^2}\Bigr)\|\nabla U^n\|^2 $$ で、$\Delta t \le 2h^2/C_{\mathrm{inv}}^2$ なら括弧の中は $0$ 以上です。 剛性行列の[条件数](/mathematics/numerical-linear-algebra/norms-and-conditioning#def-condition-number)が $h^{-2}$ に比例することは、1 次元（前章）でも 2 次元でも同じです。未知数の個数は $h^{-2}$ に増えますが、条件数は次元によりません。4 は、陽的な時間積分の安定条件が、差分法の $\Delta t \le h^2/4$（2 次元）と同じく $\Delta t = O(h^2)$ であることを示しています。しかも有限要素法の陽的 Euler 法は、毎ステップ質量行列 $G$ の連立方程式を解く必要があります（集中化して対角行列にすれば不要になります）。 ## 7. 放物型問題 熱方程式 $u_t - \Delta u = f$（$\Omega\times(0, T_{\mathrm{f}}]$）、$u = 0$（$\partial\Omega$）、$u(0) = u_0$ を、空間だけ有限要素法で離散化します。 $u_h : [0, T_{\mathrm{f}}] \to V_h$ で $$ \langle\partial_tu_h(t), v_h\rangle + \langle\nabla u_h(t), \nabla v_h\rangle = \langle f(t), v_h\rangle\quad(v_h \in V_h),\qquad u_h(0) = u_{0h} $$ を満たすものを**半離散**の有限要素解という。$u_h = \sum_jc_j(t)\varphi_j$ と書くと、これは常微分方程式系 $G\dot{\boldsymbol{c}} + K\boldsymbol{c} = \boldsymbol{b}(t)$ である。時間刻み $\Delta t$、$t_n = n\Delta t$、$0 \le \theta \le 1$ に対し、$U^n \in V_h$ を $$ \Bigl\langle\frac{U^{n+1} - U^n}{\Delta t}, v_h\Bigr\rangle + \bigl\langle\nabla(\theta U^{n+1} + (1-\theta)U^n), \nabla v_h\bigr\rangle = \bigl\langle\theta f(t_{n+1}) + (1-\theta)f(t_n), v_h\bigr\rangle\quad(v_h \in V_h) $$ で定める方法を**全離散**の方法といい、$\theta = 1$ を陰的 Euler 法、$\theta = \frac12$ を Crank–Nicolson 法という。 $G$ は正定値なので、$f \in C([0,T_{\mathrm{f}}]; L^2(\Omega))$ なら半離散の解はただ 1 つ存在します（線形の常微分方程式系 $\dot{\boldsymbol{c}} = -G^{-1}K\boldsymbol{c} + G^{-1}\boldsymbol{b}(t)$。）。全離散の方法は、行列で書けば $(G + \theta\Delta tK)\boldsymbol{c}^{n+1} = (G - (1-\theta)\Delta tK)\boldsymbol{c}^n + \Delta t\,\boldsymbol{b}$ で、前々章の [θ 法](/computer-science/pde-discretization/elliptic-parabolic-difference#def-theta-method)の単位行列が質量行列に置き換わったものです。 次の定理の証明では、誤差 $u_h - u$ を「有限要素空間の中の差 $\theta = u_h - R_hu$」と「Ritz 射影の誤差 $\rho = R_hu - u$」に分けます。$\rho$ は §5 の結果でそのまま抑えられ、$\theta$ は有限要素空間の中の方程式を満たすので、$\rho_t$ を右辺とする簡単な評価に帰着します。 $\Omega$、$\mathcal{T}_h$ は の仮定（形状正則性と (R)）を満たすとし、$C_A = 9\sigma^2C_R$ とする。$u : [0, T_{\mathrm{f}}] \to \mathcal{H}^2_{\mathrm{D}}(\Omega)$ は $H^2(\Omega)$ のノルムで連続微分可能で、$f \in C([0,T_{\mathrm{f}}]; L^2(\Omega))$、すべての $t$ と $v \in H^1_0(\Omega)$ について $\langle u_t, v\rangle + \langle\nabla u, \nabla v\rangle = \langle f, v\rangle$ を満たすとする。 1. （半離散）$\displaystyle\|u_h(t) - u(t)\| \le \|u_{0h} - R_hu_0\| + C_Ah^2\Bigl(|u(t)|_{H^2} + \int_0^t|u_t(s)|_{H^2}\,ds\Bigr)$。特に $u_{0h} = R_hu_0$ なら、右辺は $C_Ah^2\bigl(|u_0|_{H^2} + 2\int_0^t|u_t|_{H^2}ds\bigr)$ 以下である。 2. （陰的 Euler）$u$ が $L^2(\Omega)$ のノルムで 2 回連続微分可能ならば、$t_n \le T_{\mathrm{f}}$ について $\|U^n - u(t_n)\| \le \|U^0 - R_hu_0\| + C_Ah^2\bigl(|u(t_n)|_{H^2} + t_n\sup_s|u_t(s)|_{H^2}\bigr) + \dfrac{\Delta t}{2}\,t_n\sup_s\|u_{tt}(s)\|$。 3. （Crank–Nicolson）$u$ が $L^2(\Omega)$ のノルムで 3 回連続微分可能ならば、 $\|U^n - u(t_n)\| \le \|U^0 - R_hu_0\| + C_Ah^2\bigl(|u(t_n)|_{H^2} + t_n\sup_s|u_t(s)|_{H^2}\bigr) + \dfrac{\Delta t^2}{12}\,t_n\sup_s\|u_{ttt}(s)\|$。 誤差を $\theta = u_h - R_hu \in V_h$ と $\rho = R_hu - u$ に分けます（Wheeler の方法。記号 $\theta$ は、この証明では θ 法のパラメータではなくこの差を表します）。 **$\rho$ の評価.** $u(t) \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$ なので から $\|\rho(t)\| \le C_Ah^2|u(t)|_{H^2}$ です。$R_h$ は $H^1_0(\Omega)$ 上の有界線形作用素（$\|\nabla R_hw\| \le \|\nabla w\|$）で、$u$ は $H^1_0(\Omega)$ 値の関数としても連続微分可能なので、$R_hu$ も連続微分可能で $\partial_t(R_hu) = R_hu_t$ です。$u_t(t)$ は $\mathcal{H}^2_{\mathrm{D}}(\Omega)$ の元の差分商の $H^2$ での極限なので $\mathcal{H}^2_{\mathrm{D}}(\Omega)$ に属し、$\|\rho_t(t)\| = \|R_hu_t - u_t\| \le C_Ah^2|u_t(t)|_{H^2}$ です。 **1. $\theta$ の方程式.** $v_h \in V_h$ について、半離散の式、Ritz 射影の定義 $\langle\nabla R_hu, \nabla v_h\rangle = \langle\nabla u, \nabla v_h\rangle$、$u$ の式を順に使うと $$ \langle\theta_t, v_h\rangle + \langle\nabla\theta, \nabla v_h\rangle = \langle f, v_h\rangle - \langle R_hu_t, v_h\rangle - \langle\nabla u, \nabla v_h\rangle = \langle u_t - R_hu_t, v_h\rangle = -\langle\rho_t, v_h\rangle $$ です。$v_h = \theta(t)$ とすると $\frac12\frac{d}{dt}\|\theta\|^2 + \|\nabla\theta\|^2 = -\langle\rho_t, \theta\rangle \le \|\rho_t\|\|\theta\|$ です。$\varepsilon > 0$ に対し $g_\varepsilon(t) = (\|\theta(t)\|^2 + \varepsilon^2)^{1/2}$ とおくと、$g_\varepsilon' = \frac{1}{2g_\varepsilon}\frac{d}{dt}\|\theta\|^2 \le \frac{\|\rho_t\|\|\theta\|}{g_\varepsilon} \le \|\rho_t\|$ なので、積分して $\varepsilon \to 0$ とすれば $$ \|\theta(t)\| \le \|\theta(0)\| + \int_0^t\|\rho_t(s)\|\,ds \le \|u_{0h} - R_hu_0\| + C_Ah^2\int_0^t|u_t(s)|_{H^2}\,ds $$ です。$\|u_h - u\| \le \|\theta\| + \|\rho\|$ に代入して前半です。後半は、$|u(t)|_{H^2} \le |u_0|_{H^2} + \int_0^t|u_t|_{H^2}ds$（$u(t) - u_0$ の 2 階導関数に、下の「平均値の不等式」と同じ議論を使う）から従います。 **平均値の不等式.** $g : [a,b] \to L^2(\Omega)$ が連続微分可能なら、任意の $w \in L^2(\Omega)$ について、スカラー関数 $s \mapsto \langle g(s), w\rangle$ に微積分学の基本定理を使って $|\langle g(b) - g(a), w\rangle| = |\int_a^b\langle g'(s), w\rangle ds| \le \int_a^b\|g'(s)\|ds\,\|w\|$ なので、$\|g(b) - g(a)\| \le \int_a^b\|g'(s)\|\,ds \le (b - a)\sup_s\|g'(s)\|$ です。 **2. 陰的 Euler.** $\theta^n = U^n - R_hu(t_n)$、$u^n = u(t_n)$ とおきます。全離散の式（$\theta = 1$）と Ritz 射影の定義、$t = t_{n+1}$ での $u$ の式から、$v_h \in V_h$ について $$ \Bigl\langle\frac{\theta^{n+1} - \theta^n}{\Delta t}, v_h\Bigr\rangle + \langle\nabla\theta^{n+1}, \nabla v_h\rangle = \langle\omega^{n+1}, v_h\rangle,\qquad\omega^{n+1} = \Bigl(u_t(t_{n+1}) - \frac{u^{n+1} - u^n}{\Delta t}\Bigr) + \frac{(I - R_h)(u^{n+1} - u^n)}{\Delta t} $$ です。第 1 項を $\omega_1$、第 2 項を $\omega_2$ とします。$w \in L^2(\Omega)$ について、スカラー関数 $\phi(s) = \langle u(s), w\rangle$ の Taylor の定理（）から $\phi'(t_{n+1}) - \frac{\phi(t_{n+1}) - \phi(t_n)}{\Delta t} = \frac{\Delta t}{2}\phi''(c)$ なので、$|\langle\omega_1, w\rangle| \le \frac{\Delta t}{2}\sup_s\|u_{tt}(s)\|\|w\|$、すなわち $\|\omega_1\| \le \frac{\Delta t}{2}\sup\|u_{tt}\|$ です。$u^{n+1} - u^n \in \mathcal{H}^2_{\mathrm{D}}(\Omega)$ に を使い、2 階導関数に平均値の不等式を使うと $\|\omega_2\| \le \frac{C_Ah^2}{\Delta t}|u^{n+1} - u^n|_{H^2} \le C_Ah^2\sup_s|u_t(s)|_{H^2}$ です。$v_h = \theta^{n+1}$ とすると $$ \|\theta^{n+1}\|^2 - \langle\theta^n, \theta^{n+1}\rangle + \Delta t\|\nabla\theta^{n+1}\|^2 = \Delta t\langle\omega^{n+1}, \theta^{n+1}\rangle $$ なので $\|\theta^{n+1}\|^2 \le (\|\theta^n\| + \Delta t\|\omega^{n+1}\|)\|\theta^{n+1}\|$、すなわち $\|\theta^{n+1}\| \le \|\theta^n\| + \Delta t\|\omega^{n+1}\|$ です。くり返して $\|\theta^n\| \le \|\theta^0\| + t_n\max_k\|\omega^k\|$ で、$\|U^n - u^n\| \le \|\theta^n\| + \|\rho(t_n)\|$ に代入して主張を得ます。 **3. Crank–Nicolson.** 同じ計算で、$t_n$ と $t_{n+1}$ での $u$ の式の平均を使うと $$ \Bigl\langle\frac{\theta^{n+1} - \theta^n}{\Delta t}, v_h\Bigr\rangle + \Bigl\langle\nabla\frac{\theta^{n+1} + \theta^n}{2}, \nabla v_h\Bigr\rangle = \langle\omega^{n+1}, v_h\rangle,\qquad\omega^{n+1} = \Bigl(\frac{u_t(t_{n+1}) + u_t(t_n)}{2} - \frac{u^{n+1} - u^n}{\Delta t}\Bigr) + \omega_2 $$ です。$\psi = \phi'$（$\phi(s) = \langle u(s), w\rangle$）について、部分積分を 2 回行うと $$ \int_{t_n}^{t_{n+1}}(s - t_n)(t_{n+1} - s)\psi''(s)\,ds = \Delta t\bigl(\psi(t_n) + \psi(t_{n+1})\bigr) - 2\int_{t_n}^{t_{n+1}}\psi(s)\,ds $$ で、$\int_{t_n}^{t_{n+1}}(s - t_n)(t_{n+1} - s)\,ds = \frac{\Delta t^3}{6}$ なので、$\bigl|\frac{\psi(t_n) + \psi(t_{n+1})}{2} - \frac{\phi(t_{n+1}) - \phi(t_n)}{\Delta t}\bigr| \le \frac{\Delta t^2}{12}\sup|\psi''|$ です（台形則の誤差）。したがって $\|\omega_1\| \le \frac{\Delta t^2}{12}\sup\|u_{ttt}\|$ です。$v_h = \frac{\theta^{n+1} + \theta^n}{2}$ とすると $$ \frac{\|\theta^{n+1}\|^2 - \|\theta^n\|^2}{2\Delta t} + \Bigl\|\nabla\frac{\theta^{n+1} + \theta^n}{2}\Bigr\|^2 \le \|\omega^{n+1}\|\,\frac{\|\theta^{n+1}\| + \|\theta^n\|}{2} $$ で、$\|\theta^{n+1}\| + \|\theta^n\| > 0$ なら割って $\|\theta^{n+1}\| - \|\theta^n\| \le \Delta t\|\omega^{n+1}\|$（$0$ なら自明）です。あとは 2 と同じです。 証明の構造は、前々章の差分法と同じ「局所誤差 $\omega$ の積み上げ」で、安定性は試験関数に解そのものをとるエネルギー法から出ています。違いは、空間の誤差を Ritz 射影 $R_hu$ との比較にまとめたことです。$u$ と直接比べると $\|\nabla(u - u_h)\|$ の $O(h)$ が現れてしまいますが、$R_hu$ と比べると、Ritz 射影の定義によって勾配の項が消え、$L^2$ の $O(h^2)$ だけが残ります。初期値 $u_{0h}$ は、$R_hu_0$ の代わりに節点補間 $I_hu_0$ をとっても、$\|I_hu_0 - R_hu_0\| \le \|I_hu_0 - u_0\| + \|u_0 - R_hu_0\| = O(h^2)$（ と ）なので、次数は変わりません。連続の問題の最大値原理（）の離散版は、質量行列を集中化し、 の角度の条件を仮定すると、前々章の と同じ議論で得られます。 区分 2 次以上の要素では、解が滑らかなら次数が上がります（$k$ 次要素で $H^1$ の誤差 $O(h^k)$）。四角形要素、曲がった境界を近似するアイソパラメトリック要素、計算した解から誤差を見積もる事後誤差評価と、それに基づく適応的な細分、最大値ノルムでの誤差評価（$O(h^2|\log h|)$）は、どれも証明に新しい道具が要るので、この章では扱いません。Brenner–Scott と Ciarlet（参考文献）に証明があります。有限要素法による固有値の計算と、その誤差定数を使った固有値の下界は、[固有値の包含の章](/mathematics/computer-assisted-proofs/eigenvalue-enclosures) で扱っています。そこで必要になる「射影の誤差定数」は、 の前半の形の不等式 $\|w - R_hw\| \le C_h\|\nabla(w - R_hw)\|$ の定数です。 ## 8. 数値実験 まず、以下の実験で共通に使う部品をまとめます。格子から三角形分割を作る関数、剛性行列・質量行列・荷重ベクトルを組み立てる関数、Dirichlet 問題を解く関数、誤差を計算する関数です。組み立てでは、要素ごとの $3 \times 3$ の行列の成分を、行と列の節点番号とともに 1 本の配列に並べ、`coo_matrix` に渡します。同じ位置に複数の値があると足し合わされるので、これがそのまま組み立てになります。このファイルを `fem2d.py` として保存し、後のコードから読み込みます。 ```python # 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 def grid_mesh(xs, ys, keep=lambda xc, yc: True): """格子の各セル（中心が keep を満たすもの）を対角線で 2 つの三角形に分ける。""" nx, ny = len(xs) - 1, len(ys) - 1 idx = lambda i, j: i * (ny + 1) + j X, Y = np.meshgrid(xs, ys, indexing="ij") tri = [] for i in range(nx): for j in range(ny): if keep(0.5 * (xs[i] + xs[i + 1]), 0.5 * (ys[j] + ys[j + 1])): tri.append([idx(i, j), idx(i + 1, j), idx(i + 1, j + 1)]) tri.append([idx(i, j), idx(i + 1, j + 1), idx(i, j + 1)]) tri = np.array(tri) used = np.unique(tri) # 使われる節点だけに番号を振り直す new = -np.ones((nx + 1) * (ny + 1), dtype=int) new[used] = np.arange(len(used)) return np.column_stack([X.ravel(), Y.ravel()])[used], new[tri] def geometry(p, t): """各三角形の面積と、重心座標の勾配（3 頂点ぶん）。""" a, b, c = p[t[:, 0]], p[t[:, 1]], p[t[:, 2]] det = (b[:, 0] - a[:, 0]) * (c[:, 1] - a[:, 1]) - (b[:, 1] - a[:, 1]) * (c[:, 0] - a[:, 0]) g = np.stack([np.column_stack([b[:, 1] - c[:, 1], c[:, 0] - b[:, 0]]), np.column_stack([c[:, 1] - a[:, 1], a[:, 0] - c[:, 0]]), np.column_stack([a[:, 1] - b[:, 1], b[:, 0] - a[:, 0]])], axis=1) return 0.5 * np.abs(det), g / det[:, None, None] QP = np.array([[2/3, 1/6, 1/6], [1/6, 2/3, 1/6], [1/6, 1/6, 2/3]]) # 2 次まで正確な 3 点則 def assemble(p, t, f=None): """剛性行列 K、質量行列 G、荷重ベクトル b（coo_matrix が重複した添字を足し合わせる）。""" area, g = geometry(p, t) n = len(p) I = np.repeat(t, 3, axis=1).ravel() J = np.tile(t, (1, 3)).ravel() Ke = area[:, None, None] * np.einsum("eid,ejd->eij", g, g) Ge = area[:, None, None] / 12.0 * (np.ones((3, 3)) + np.eye(3)) K = sp.coo_matrix((Ke.ravel(), (I, J)), shape=(n, n)).tocsr() G = sp.coo_matrix((Ge.ravel(), (I, J)), shape=(n, n)).tocsr() b = np.zeros(n) if f is not None: for lam in QP: # 求積点 x = Σ λ_i a_i xq = np.einsum("i,eid->ed", lam, p[t]) np.add.at(b, t, area[:, None] / 3.0 * f(xq[:, 0], xq[:, 1])[:, None] * lam[None, :]) return K, G, b def boundary_nodes(t, n): """1 つの三角形にしか属さない辺の端点。""" e = np.sort(np.vstack([t[:, [0, 1]], t[:, [1, 2]], t[:, [2, 0]]]), axis=1) ue, cnt = np.unique(e, axis=0, return_counts=True) return np.unique(ue[cnt == 1]) def solve_dirichlet(p, t, f, gD): K, G, b = assemble(p, t, f) bd = boundary_nodes(t, len(p)) it = np.setdiff1d(np.arange(len(p)), bd) u = np.zeros(len(p)) u[bd] = gD(p[bd, 0], p[bd, 1]) u[it] = spsolve(K[it][:, it].tocsc(), b[it] - K[it][:, bd].dot(u[bd])) return u, K, G def errors(p, t, u, uex, grad): """L2 誤差と H1 半ノルム誤差（3 点則）。""" area, g = geometry(p, t) gu = np.einsum("ei,eid->ed", u[t], g) e0 = e1 = 0.0 for lam in QP: xq = np.einsum("i,eid->ed", lam, p[t]) e0 += np.sum(area / 3.0 * (uex(xq[:, 0], xq[:, 1]) - np.einsum("ei,i->e", u[t], lam))**2) gx, gy = grad(xq[:, 0], xq[:, 1]) e1 += np.sum(area / 3.0 * ((gx - gu[:, 0])**2 + (gy - gu[:, 1])**2)) return np.sqrt(e0), np.sqrt(e1) ``` 要素剛性行列は $|T|\,\nabla\lambda_i\cdot\nabla\lambda_j$、要素質量行列は $\frac{|T|}{12}(1 + \delta_{ij})$（ の証明）です。重心座標の勾配は、頂点 $a, b, c$ について $\nabla\lambda_a = (b_2 - c_2,\ c_1 - b_1)/\det$ などで、$\det$ は $B$ の行列式です。 **正方形での製造解.** $\Omega = (0,1)^2$、$u = \sin\pi x\sin\pi y$、$f = 2\pi^2u$ を、一辺 $1/n$ の直角二等辺三角形の分割で解きます。 ```python from fem2d import * pi = np.pi uex = lambda x, y: np.sin(pi * x) * np.sin(pi * y) grad = lambda x, y: (pi * np.cos(pi * x) * np.sin(pi * y), pi * np.sin(pi * x) * np.cos(pi * y)) f = lambda x, y: 2 * pi**2 * uex(x, y) E0, E1, EP = [], [], [] for n in (8, 16, 32, 64, 128): p, t = grid_mesh(np.linspace(0, 1, n + 1), np.linspace(0, 1, n + 1)) u, K, G = solve_dirichlet(p, t, f, lambda x, y: 0 * x) e0, e1 = errors(p, t, u, uex, grad) E0.append(e0); E1.append(e1) EP.append(np.sqrt(pi**2 / 2 - u.dot(K.dot(u)))) # ‖∇u‖² − ‖∇u_h‖² の平方根 if n == 8: # 5 点差分との照合 i = np.argmin(np.abs(p[:, 0] - 0.5) + np.abs(p[:, 1] - 0.5)) row = K[i].toarray().ravel() print("中央の節点の行の非零成分:", np.round(np.sort(row[np.abs(row) > 1e-12]), 12)) for lab, E in (("H1", E1), ("L2", E0), ("Pyth", EP)): o = [np.log2(E[i] / E[i + 1]) for i in range(4)] print(f"{lab:4s} " + " ".join(f"{v:9.2e}" for v in E) + " |" + " ".join(f"{v:5.2f}" for v in o)) ``` ```text 中央の節点の行の非零成分: [-1. -1. -1. -1. 4.] H1 4.32e-01 2.18e-01 1.09e-01 5.45e-02 2.73e-02 | 0.99 1.00 1.00 1.00 L2 2.04e-02 5.20e-03 1.31e-03 3.27e-04 8.18e-05 | 1.97 1.99 2.00 2.00 Pyth 4.31e-01 2.17e-01 1.09e-01 5.45e-02 2.73e-02 | 0.99 1.00 1.00 1.00 ``` **検算.** (a) 参照解は製造解です。(b) 観測次数は $H^1$ で $1$、$L^2$ で $2$ で、 と （正方形では (R) が $C_R = 1$ で成り立つ）に合います。直角二等辺三角形では $\sigma = \sqrt2/(2 - \sqrt2) = 2.41$、$|u|_{H^2} = \pi^2$（$\|u_{xx}\|^2 + 2\|u_{xy}\|^2 + \|u_{yy}\|^2 = \pi^4(\frac14 + \frac24 + \frac14)$）で、$n = 8$（$h = \sqrt2/8$）での上界は $3\sigma h|u|_{H^2} = 12.6$ と $9\sigma^2h^2|u|_{H^2} = 16.2$ です。観測値 $0.432$、$0.0204$ はその内側ですが、上界はかなり粗い値です。この章の定数は、どんな三角形にも通用するように作ったためです。(c) 別の方法による照合として、Galerkin 直交性から $\|\nabla(u - u_h)\|^2 = \|\nabla u\|^2 - \|\nabla u_h\|^2$（前章の の 3 と同じ内容）で、$\|\nabla u\|^2 = \pi^2/2$、$\|\nabla u_h\|^2 = \boldsymbol{c}^{\mathsf{T}}\tilde K\boldsymbol{c}$ です。`Pyth` の行がこれで、求積で直接計算した `H1` の行と 3 桁まで一致します（荷重と誤差の計算に使った求積則の誤差の範囲）。(d) 中央の節点の行は $4$ と $-1$ が 4 つで、 のとおり 5 点差分です。 **L 字形の領域.** $(-1,1)^2$ から第 4 象限 $[0,1)\times(-1,0]$ を除いた領域で、調和関数 $u = r^{2/3}\sin(2\theta/3)$（$0 \le \theta \le 3\pi/2$）を境界値として与え、$f = 0$ で解きます。 ```python from fem2d import * def uex(x, y): # r^(2/3) sin(2θ/3)、θ ∈ [0, 3π/2] th = np.mod(np.arctan2(y, x), 2 * np.pi) return (x**2 + y**2)**(1 / 3) * np.sin(2 * th / 3) def grad(x, y): r, th = np.sqrt(x**2 + y**2), np.mod(np.arctan2(y, x), 2 * np.pi) ur, ut = 2 / 3 * r**(-1 / 3) * np.sin(2 * th / 3), 2 / 3 * r**(-1 / 3) * np.cos(2 * th / 3) return ur * np.cos(th) - ut * np.sin(th), ur * np.sin(th) + ut * np.cos(th) keep = lambda xc, yc: not (xc > 0 and yc < 0) # L 字形: (-1,1)^2 から第 4 象限を除く E0, E1 = [], [] for n in (8, 16, 32, 64, 128): xs = np.linspace(-1, 1, 2 * n + 1) p, t = grid_mesh(xs, xs, keep) u, K, G = solve_dirichlet(p, t, None, uex) e0, e1 = errors(p, t, u, uex, grad) E0.append(e0); E1.append(e1) for lab, E in (("H1", E1), ("L2", E0)): o = [np.log2(E[i] / E[i + 1]) for i in range(4)] print(f"{lab} " + " ".join(f"{v:9.2e}" for v in E) + " |" + " ".join(f"{v:5.2f}" for v in o)) ``` ```text H1 1.15e-01 7.34e-02 4.67e-02 2.96e-02 1.87e-02 | 0.64 0.65 0.66 0.66 L2 6.50e-03 2.68e-03 1.09e-03 4.43e-04 1.78e-04 | 1.28 1.29 1.30 1.31 ``` 観測次数は $H^1$ で $0.66$、$L^2$ で $1.31$ に近づき、 の $2/3$ と $4/3$ に合います。解が $H^2$ に属さないので の仮定が成り立たず、次数 $1$ と $2$ は得られません。境界値が $0$ でない問題は、境界の節点の値を固定し、その寄与を右辺に移して解いています（`solve_dirichlet`）。**検算**として、$u$ が調和であることは、極座標で $\Delta(r^\alpha\sin\alpha\theta) = (\alpha(\alpha - 1) + \alpha - \alpha^2)r^{\alpha - 2}\sin\alpha\theta = 0$ から、凹角をはさむ 2 辺（$\theta = 0$ と $\theta = 3\pi/2$）で $0$ であることは $\sin0 = \sin\pi = 0$ から確かめられます。誤差の積分は特異点の近くで求積の誤差を含みますが、次数の読み取りには影響しない大きさです（$n$ を変えても次数が安定しています）。 単位正方形を $n_x \times n_y$ 個の長方形（横 $h_x$、縦 $h_y$）に分け、各長方形を中心の点と 4 頂点で 4 つの三角形に分けます。水平な辺は、上下の長方形の中心を頂点とする 2 つの三角形に共有され、向かい合う角 $\theta$ は $\tan\frac\theta2 = h_x/h_y$ を満たします。$h_x > h_y$ なら $\theta > \pi/2$（鈍角）で、和 $2\theta$ は $\pi$ を超えます。 の 1 から、この辺に対応する成分は $-\cot\theta = \frac{r^2 - 1}{2r}$（$r = h_x/h_y$）で、正です。 ```python from fem2d import * def crisscross(nx, ny): """各セルを中心の節点と 4 頂点で 4 つの三角形に分ける。""" xs, ys = np.linspace(0, 1, nx + 1), np.linspace(0, 1, ny + 1) X, Y = np.meshgrid(xs, ys, indexing="ij") Xc, Yc = np.meshgrid(0.5 * (xs[:-1] + xs[1:]), 0.5 * (ys[:-1] + ys[1:]), indexing="ij") p = np.vstack([np.column_stack([X.ravel(), Y.ravel()]), np.column_stack([Xc.ravel(), Yc.ravel()])]) v = lambda i, j: i * (ny + 1) + j c = lambda i, j: (nx + 1) * (ny + 1) + i * ny + j t = [[a, b, c(i, j)] for i in range(nx) for j in range(ny) for a, b in ((v(i, j), v(i + 1, j)), (v(i + 1, j), v(i + 1, j + 1)), (v(i + 1, j + 1), v(i, j + 1)), (v(i, j + 1), v(i, j)))] return p, np.array(t) for nx, ny in ((8, 8), (4, 12), (4, 24)): p, t = crisscross(nx, ny) K, G, _ = assemble(p, t) bd = boundary_nodes(t, len(p)) it = np.setdiff1d(np.arange(len(p)), bd) Kii = K[it][:, it].toarray() off = Kii - np.diag(np.diag(Kii)) Kinv = np.linalg.inv(Kii) print(f"セル {nx}x{ny} (hx/hy={ny / nx:3.1f}) 非対角成分の最大={off.max():7.4f}" f" K^-1 の最小成分={Kinv.min():9.2e}") ``` ```text セル 8x8 (hx/hy=1.0) 非対角成分の最大= 0.0000 K^-1 の最小成分= 5.86e-05 セル 4x12 (hx/hy=3.0) 非対角成分の最大= 1.3333 K^-1 の最小成分= 2.87e-05 セル 4x24 (hx/hy=6.0) 非対角成分の最大= 2.9167 K^-1 の最小成分=-1.40e-02 ``` 非対角成分の最大は公式の値 $\frac{r^2 - 1}{2r}$（$r = 3$ で $\frac{8}{6} = 1.3333$、$r = 6$ で $\frac{35}{12} = 2.9167$）と一致します。$r = 6$ では $K^{-1}$ に負の成分が現れます。これは、ある 1 つの節点に正の荷重をかけた（$\boldsymbol{b} \ge \boldsymbol{0}$）のに、離散解 $K^{-1}\boldsymbol{b}$ がどこかで負になることを意味し、最大値原理（の最小についての形）が破れています。連続の問題では $f \ge 0$ なら $u \ge 0$（）なので、これは離散化が作った現象です。$r = 3$ では非対角成分に正の値があっても $K^{-1}$ は非負のままで、 の角度の条件が十分条件であって必要条件ではないことが分かります。$r = 1$ は直角二等辺三角形の分割で、条件を満たします。 **熱方程式の全離散.** $\Omega = (0,1)^2$、$u_0 = \sin\pi x\sin\pi y$（厳密解は $e^{-2\pi^2t}u_0$）を、$\Delta t = h/4$（$h = 1/n$ は格子の一辺）として $T_{\mathrm{f}} = 0.125$ まで解きます。初期値は節点補間です。 ```python from fem2d import * from scipy.sparse.linalg import splu pi, T = np.pi, 0.125 u0 = lambda x, y: np.sin(pi * x) * np.sin(pi * y) def heat_theta(n, theta): p, t = grid_mesh(np.linspace(0, 1, n + 1), np.linspace(0, 1, n + 1)) K, G, _ = assemble(p, t) it = np.setdiff1d(np.arange(len(p)), boundary_nodes(t, len(p))) K, G = K[it][:, it].tocsc(), G[it][:, it].tocsc() steps = n // 2 # dt = T / steps = h / 4 dt = T / steps lu = splu((G + theta * dt * K).tocsc()) U = u0(p[it, 0], p[it, 1]) for _ in range(steps): U = lu.solve((G - (1 - theta) * dt * K).dot(U)) u = np.zeros(len(p)); u[it] = U ex = lambda x, y: np.exp(-2 * pi**2 * T) * u0(x, y) return errors(p, t, u, ex, lambda x, y: (0 * x, 0 * x))[0] for lab, th in (("CN dt=h/4", 0.5), ("BE dt=h/4", 1.0)): E = [heat_theta(n, th) for n in (8, 16, 32, 64)] o = [np.log2(E[i] / E[i + 1]) for i in range(3)] print(f"{lab} " + " ".join(f"{v:9.2e}" for v in E) + " |" + " ".join(f"{v:5.2f}" for v in o)) ``` ```text CN dt=h/4 8.23e-03 2.10e-03 5.29e-04 1.32e-04 | 1.97 1.99 2.00 BE dt=h/4 2.49e-02 1.44e-02 7.65e-03 3.93e-03 | 0.79 0.91 0.96 ``` **検算.** (a) 参照解は固有関数展開の閉じた形で、$T_{\mathrm{f}} = 0.125$ で $e^{-2\pi^2T_{\mathrm{f}}} = 0.0848$ です。(b) $\Delta t$ を $h$ に比例させているので、 の評価は Crank–Nicolson で $O(h^2 + \Delta t^2) = O(h^2)$、陰的 Euler で $O(h^2 + \Delta t) = O(h)$ で、観測次数 $2$ と $1$（に近づく値）に合います。(c) 行列 $G + \theta\Delta tK$ は $n$ によらず 1 回だけ LU 分解し、毎ステップは代入だけです。この行列は帯幅が $n$ 程度の疎行列で、前々章の 1 次元の三重対角行列にあたります。 ## 9. 演習 参照三角形 $\hat T$（頂点 $a_0 = (0,0)$、$a_1 = (1,0)$、$a_2 = (0,1)$）の重心座標とその勾配を求め、要素剛性行列 $\bigl(\int_{\hat T}\nabla\lambda_j\cdot\nabla\lambda_i\bigr)_{i,j}$ を計算せよ。結果を余接の公式（ の 1）で確かめよ。 $\lambda_1 = x_1$、$\lambda_2 = x_2$、$\lambda_0 = 1 - x_1 - x_2$ で、勾配は $\nabla\lambda_1 = (1,0)$、$\nabla\lambda_2 = (0,1)$、$\nabla\lambda_0 = (-1,-1)$ です。$|\hat T| = \frac12$ を掛けて $$ K_{\hat T} = \frac12\begin{pmatrix}2 & -1 & -1\\ -1 & 1 & 0\\ -1 & 0 & 1\end{pmatrix} $$ です（行と列は $a_0, a_1, a_2$ の順）。角は $a_0$ で $\frac\pi2$、$a_1$ と $a_2$ で $\frac\pi4$ です。余接の公式では、$(a_1, a_2)$ の成分は向かい合う $a_0$ の角から $-\frac12\cot\frac\pi2 = 0$、$(a_0, a_1)$ の成分は向かい合う $a_2$ の角から $-\frac12\cot\frac\pi4 = -\frac12$ で、一致します。対角成分は行の和が $0$ であることから決まります。 $f$ が定数のとき、荷重ベクトルの成分は $b_i = \frac{f}{3}\sum_{T \ni z_i}|T|$ であることを示せ。直角二等辺三角形の規則的な分割（）では $b_i = fh^2$ になることを確かめよ。 $b_i = f\int_\Omega\varphi_i\,dx = f\sum_{T\ni z_i}\int_T\lambda_i\,dx$ です。$\int_T\lambda_i\,dx$ は 3 つの重心座標について等しく（ の証明と同じく、参照三角形に写して頂点の番号を付け替えればよい）、和は $\int_T1\,dx = |T|$ なので、$\int_T\lambda_i\,dx = |T|/3$ です。規則的な分割では、内部の節点は面積 $h^2/2$ の三角形 6 個に共有されるので $b_i = \frac f3\cdot6\cdot\frac{h^2}{2} = fh^2$ です。 半離散の有限要素解（）が $\|u_h(t)\| \le \|u_{0h}\| + \int_0^t\|f(s)\|\,ds$ を満たすことを示せ。$f = 0$ のとき、$\Omega$ が幅 $d$ の帯に含まれるなら $\|u_h(t)\| \le e^{-t/d^2}\|u_{0h}\|$ であることも示せ。 半離散の式で $v_h = u_h(t)$ とすると $\frac12\frac{d}{dt}\|u_h\|^2 + \|\nabla u_h\|^2 = \langle f, u_h\rangle \le \|f\|\|u_h\|$ です。 の証明の 1 と同じく、$g_\varepsilon = (\|u_h\|^2 + \varepsilon^2)^{1/2}$ について $g_\varepsilon' \le \|f\|$ なので、積分して $\varepsilon \to 0$ とすれば前半です。 $f = 0$ のとき、Poincaré の不等式（。$u_h(t) \in V_h \subset H^1_0(\Omega)$）から $\|\nabla u_h\|^2 \ge d^{-2}\|u_h\|^2$ なので、$y(t) = \|u_h(t)\|^2$ は $y' \le -\frac{2}{d^2}y$ を満たします。$\frac{d}{dt}\bigl(e^{2t/d^2}y\bigr) \le 0$ から $y(t) \le e^{-2t/d^2}y(0)$ で、平方根をとって後半です。この評価は $h$ によらず、分割の形にもよりません。半離散の方法の安定性は、楕円型の場合と同じく、連続の問題のエネルギー評価からそのまま受け継がれます。 一辺 $h$ の正三角形で平面を埋めた分割で、内部の節点 $z_i$ の行の成分 $\tilde K_{ii}$、$\tilde K_{ij}$ を求めよ。滑らかな $u$ の節点での値にこの行を作用させたものを $\int\varphi_i = \frac{\sqrt3}{2}h^2$ で割ると、$-\Delta u(z_i)$ を $O(h^2)$ の誤差で近似することを示せ。 $z_i$ から 6 本の辺が出ていて、各辺に向かい合う角は両側とも $\frac\pi3$ です。余接の公式から $\tilde K_{ij} = -\frac12\cdot2\cot\frac\pi3 = -\frac{1}{\sqrt3}$ で、行の和が $0$ なので $\tilde K_{ii} = \frac{6}{\sqrt3} = 2\sqrt3$ です。$z_i$ を頂点とする三角形は 6 個で、面積は $\frac{\sqrt3}{4}h^2$ なので $\int\varphi_i = \frac13\cdot6\cdot\frac{\sqrt3}{4}h^2 = \frac{\sqrt3}{2}h^2$ です。 隣の節点は $z_i + h\boldsymbol{n}_k$、$\boldsymbol{n}_k = (\cos\frac{k\pi}{3}, \sin\frac{k\pi}{3})$（$k = 0, \ldots, 5$）です。$u \in C^4$ として Taylor 展開すると、$\boldsymbol{n}_{k+3} = -\boldsymbol{n}_k$ なので奇数次の項は打ち消し合い、 $$ \sum_{k=0}^{5}\bigl(u(z_i + h\boldsymbol{n}_k) - u(z_i)\bigr) = \frac{h^2}{2}\sum_{k=0}^{5}\boldsymbol{n}_k^{\mathsf{T}}D^2u(z_i)\,\boldsymbol{n}_k + O(h^4) $$ です。$\sum_k\boldsymbol{n}_k\boldsymbol{n}_k^{\mathsf{T}}$ は、$\sum_k\cos^2\frac{k\pi}{3} = \sum_k\sin^2\frac{k\pi}{3} = 3$、$\sum_k\cos\frac{k\pi}{3}\sin\frac{k\pi}{3} = 0$ から $3I$ なので、右辺は $\frac{3h^2}{2}\Delta u(z_i) + O(h^4)$ です。したがって $$ \frac{1}{\frac{\sqrt3}{2}h^2}\sum_j\tilde K_{ij}u(z_j) = \frac{2}{\sqrt3h^2}\cdot\Bigl(-\frac{1}{\sqrt3}\Bigr)\sum_k\bigl(u(z_i + h\boldsymbol{n}_k) - u(z_i)\bigr) = -\Delta u(z_i) + O(h^2) $$ です。正三角形の分割の有限要素法は、6 つの隣接点を使う 7 点差分になります。 ## 参考文献 - S. C. Brenner, L. R. Scott『The Mathematical Theory of Finite Element Methods』3rd ed., Springer, 2008 — 第 4 章（平均 Taylor 多項式と Bramble–Hilbert の補題、補間誤差）、第 5 章（多角形領域での誤差評価と双対論法）。 - P. G. Ciarlet『The Finite Element Method for Elliptic Problems』North-Holland, 1978（SIAM Classics, 2002）— アフィン同値な要素族、補間誤差、逆不等式、離散最大値原理。 - V. Thomée『Galerkin Finite Element Methods for Parabolic Problems』2nd ed., Springer, 2006 — 第 1 章（半離散と全離散の標準的な誤差評価）。 - P. Grisvard『Elliptic Problems in Nonsmooth Domains』Pitman, 1985（SIAM Classics, 2011）— 凸領域の $H^2$ 正則性、角のある領域での解の特異性。 - I. Babuška, A. K. Aziz, "On the angle condition in the finite element method", *SIAM Journal on Numerical Analysis* 13 (1976), 214–226 — 最大角の条件。 - M. F. Wheeler, "A priori $L_2$ error estimates for Galerkin approximations to parabolic partial differential equations", *SIAM Journal on Numerical Analysis* 10 (1973), 723–759 — Ritz 射影との比較による放物型問題の誤差評価。 - T. Dupont, R. Scott, "Polynomial approximation of functions in Sobolev spaces", *Mathematics of Computation* 34 (1980), 441–463 — 平均 Taylor 多項式による近似。