重み付き残差法を理解する｜第1回

はじめに

有限要素法を勉強していると、割と最初のほうに重み付き残差法（特に、Galerkin法）という用語が登場します。
それが有限要素法とどう関係しているのか、よく分からなくないですか？
私は分かりませんでした。

そこでいろいろと調べてみると、重み付き残差法というのは微分方程式の数値計算の一種であって、有限要素法という文脈でしか登場しないものというわけではないようです。

このシリーズでは2回に分けて重み付き残差法やGalerkin法について私が理解したことをシェアします。
最終的には有限要素法とのつながりも説明したいですが、まずは微分方程式の解法としての重み付き残差法を説明します。

第1回のこの記事では、重み付き残差法の背景や全体像を説明します。第2回はこちらからどうぞ。

目次

解析的な微分方程式の解法
数値的な微分方程式の解法
数値代入法の問題点
積分残差
重み付き残差法

解析的な微分方程式の解法

はじめにで書いたように、重み付き残差法は微分方程式の数値解法の一つです。
そこで、（主として私のために）まずは解析的な微分方程式の解法をおさらいしておきましょう。
よくご存じの方はすっ飛ばしていただいて大丈夫です。

例題として、次の微分方程式を考えます。（問題はこちらからお借りしました）

\frac{d^2u}{dx^2} - \frac{du}{dx} - 6u = 0 \tag{1}

まずは解として $u=e^{\lambda x}$ を仮定します。
みんな大好き、天下り的発想です。
この解を元の方程式に代入しごちゃごちゃ計算すると、

(\lambda^2 - \lambda - 6) e^{\lambda x} = 0 \tag{2}

が得られます。

これは任意の $x$ で成立する「恒等式」ですが、恒等式の解き方には主に2通りあります。
「係数比較法」と「数値代入法」です。
名前から大体わかると思いますが、不安な方はこういったサイトを参考にしてください。
恒等式（2）の場合、係数比較法で

\lambda^2 - \lambda - 6 = 0 \tag{3}

が得られます。数値代入法でも、$x=0$ を代入することで同じ式が得られます。これを解いて、$\lambda = -2$, $3$ と求めれば、与えられた微分方程式（1）の解が $y=e^{-2x}$, $e^{3x}$ と求められるわけです。

今回は $u=e^{\lambda x}$ という1つの基底関数による解を仮定しましたが、解として最初に仮定する関数は、複数の基底関数の線形和でも構いません。
例えば、$u = e^{\lambda x} + a\sin(x)$ なんていうのでもいいわけです。
すなわち、$n$ 個の基底関数 $\phi_i$ とその係数 $a_i$ を用いて、

u(x) = \sum_{i=1}^n a_i \phi_i \tag{4}

と関数を定めるわけです。
基底関数の選び方は問題によります。
上の例題のように指数関数にしてもいいし、三角関数でもいいでしょう。
とにかく、正しい解になりそうな基底関数を先に選んでおきます。

この場合、未知数は基底関数の係数 $a_i$ で、これは全部で $n$ 個あります。
従って、恒等式から $n$ 個の方程式を見つけ出し、それらを連立する必要があります。
係数比較法の場合はそれぞれの基底関数の係数を比較すればよく、数値代入法の場合は $n$ 個の数値を代入すればよいです。

ここまでの解法を一般化してまとめておきましょう。

1. （任意の）$n$ 個の基底関数 $\phi_i$ を定める。
2. $u(x) = \sum_{i=1}^n a_i \phi_i$ と解を仮定する。
3. $u(x)$ を元の微分方程式に代入し、恒等式を得る。
4. 係数比較法か数値代入法を用いて、恒等式から係数 $a_i$ に関する連立方程式を $n$ 個得る。
5. 得られた連立方程式を解いて、係数 $a_i$ を求める。

以上が解析的な微分方程式の解法です。

数値的な微分方程式の解法

解析的な微分方程式の解法を確認したところで、それと同じことを数値計算でできるか考えてみましょう。
ここで問題になるのは恒等式の解き方です。
係数比較法の場合、項を基底関数ごとに整理して比較する必要がありますが、これはコンピュータさんには難しい作業です。
そこで、コンピュータさんの得意そうな数値代入法を使うことにします。

数値代入法の問題点

しかし、数値代入法には困った問題があります。
例として、$u(x)=1$ が解となる微分方程式を考えます。
その際、次の図のように仮の関数を定めると、どうなるでしょうか？

これは明らかに間違った関数です。
しかし、仮の関数と真の解が交わる点 $x=0$、$x=1.5$ 付近、$x=3$ 付近を用いて数値を代入すると、コンピュータさんはこの仮の関数を正しいものと認識してしまいます。

すなわち、「真の解」⇒「数値代入法の解」は成り立つけれども、「数値代入法の解」⇒「真の解」は成り立たないわけです。
簡単な方程式であれば数値代入法の解を用いた関数が常に微分方程式を満たすことを解析的に確認することができますが、工学で登場するような複雑な関数ではそうもいきません。

積分残差

この問題はどうすれば解決できるのでしょうか？
問題が起こる原因は、特定の数値を代入していたことでした。
あらゆる $x$ を代入して確認することはできますが、計算時間の観点から現実的ではありません。

そこで、元の微分方程式に代入した後の値を全ての地点でをチェックするのではなく、それが全領域で平均的に0になってればいいだろうと近似的に考えます。
すなわち、積分して0になるか確認するわけです。
あくまでも近似的な考え方であることにご注意ください。

次の形で表される微分方程式を考えます。

\frac{d^2u(x)}{dx^2}-f(x) = 0 \tag{5}

先ほどと同様に、仮の解を

u(x) = \sum_{i=1}^n a_i \phi_i \tag{6}

とおきましょう。
繰り返しになりますが、$\phi_i$ は自分で定めた基底関数、$a_i$ は最終的に求める未知数です。
あくまで仮の答えなので、式（6）を式（5）に代入してもその右辺がゼロになるとは限りません。
そこで、”残”ってしまった解析解との”差”という意味で「残差（residual）」$r(\boldsymbol{a}, x)$ を次のように定義します。

r(\boldsymbol{a}, x) = \frac{d^2}{dx^2} \sum_{i=1}^n a_i\phi_i - f(x) \tag{7}

$\boldsymbol{a}$ は $a_1$ から $a_n$ までをまとめて1つのベクトルにしたものです。
表記上の簡略化が目的であって、深い意味はありません。
先ほど書いたように、式（7）を積分してその値を評価することにしましょう。
その積分値を「積分残差」$R(\boldsymbol{a})$ と呼ぶことにします。

R(\boldsymbol{a}) = \int_{-\infty}^{\infty} r(\boldsymbol{a}, x) dx \tag{8}

ややこしいのは、$r$ も $R$ も「残差」と呼ばれることです。
ここでは、$r$ を「残差」、$R$ を「積分残差」として区別します（積分残差は一般的な呼び方ではないと思うのでご注意ください）。

重み付き残差法

$R(\boldsymbol{a})=0$ からを用いれば基底関数の係数 $a_i$ に関する式が1つ得られます。
しかし、$a_i$ は基底関数の数と同じく $n$ 個分あるので、$n$ 個の式を求めて連立する必要がありました。
そこで登場するのが「重み付き残差法（Method of Weighted Residuals: MWR）」です。

まず、任意の関数 $w(x)$ を用意します。
仮の関数が正しければ $r(\boldsymbol{a}, x)=0$ となり、当然 $r(\boldsymbol{a}, x)w(x)=0$ となるはずです。
従って、これを積分してもゼロになります。

R(\boldsymbol{a}) = \int_{-\infty}^{\infty} r(\boldsymbol{a}, x)w(x)dx = 0 \tag{9}

一方で、式（8）と式（9）は積分の中身が異なるので積分後の式は違った形になります。
すなわち、係数 $a_i$ に関する別の式が得られたわけですね。

次に先ほどの関数 $w(x)$ とは異なった関数 $w_1(x)$ を用意します。
こちらも同様に、

R(\boldsymbol{a}) = \int_{-\infty}^{\infty} r(\boldsymbol{a}, x)w_1(x)dx = 0 \tag{10}

となることから係数 $a_i$ に関する別の式を得られます。

これを $n$ 個の式が得られるまで繰り返せば、$a_i$ についての連立方程式が完成し、未知数 $a_i$ を求めることができます！

すなわち、異なった形の $w_i(x)$ を $n$ 個用意し、それぞれの $w_i(x)$ を使って計算した積分残差 $R_i$ がゼロとなる条件から、$a_i$ に関する $n$ 個の方程式を得られるのです。

\begin{align}
R_1(\boldsymbol{a}) &= \int_{-\infty}^{\infty} r(\boldsymbol{a}, x)w_1(x)dx = 0 \\
R_2(\boldsymbol{a}) &= \int_{-\infty}^{\infty} r(\boldsymbol{a}, x)w_2(x)dx = 0 \\
&\vdots \\
R_n(\boldsymbol{a}) &= \int_{-\infty}^{\infty} r(\boldsymbol{a}, x)w_n(x)dx = 0
\end{align}

式（9）で表される積分残差 $R(\boldsymbol{a})$ は残差 $r(\boldsymbol{a}, x)$ に $w(x)$ だけの"重み"を付けて計算した残差なので、「重み付き残差」と呼びます（この記事のルールに則ると「重み付き積分残差」と呼びたいところですが、これについては一般的な慣習にに倣うことにします）。
この時、"重み"として用いる関数 $w(x)$ は「重み関数（weighting function）」と呼ばれます。　

ちなみに、式（8）で定義した積分残差は、$w(x)=1$ の場合の重み付き残差ともいえます。

解析的な解法と同様に、重み付き残差法の解法を一般化してまとめておきましょう。

1. （任意の）$n$ 個の基底関数 $\phi_i$ を定める。
2. $u(x) = \sum_{i=1}^n a_i \phi_i$ と解を仮定する。
3. $u(x)$ を元の微分方程式に代入し、残差 $r(\boldsymbol{a}, x)$ を求める。

       r(\boldsymbol{a}, x) = \frac{d^2}{dx^2} \sum_{i=1}^n a_i\phi_i - f(x) \tag{11}
   

4. $n$ 個の重み関数 $w_j(x)$ を用意し、重み付き残差 $R_j(\boldsymbol{a})$ を求める。

       R_j(\boldsymbol{a}) = \int_{-\infty}^{\infty} r(\boldsymbol{a}, x)w_j(x)dx \tag{12}
   

5. 各 $R_j(\boldsymbol{a})$ がゼロとなる条件から、$a_i$ に関する $n$ 個の連立方程式を得る。
6. 連立方程式を解いて $a_i$ を求める。

以上が重み付き残差法の説明です。
では重み付き残差法で鍵を握る重み関数はどのように選べばいいのでしょうか？
次の記事↓ではこの点について解説します。

参考資料