有限要素法を学ぶ 「重み付き残差法」

重み付き残差法について

はじめに

ここでは重み付き残差法について具体例で学んで行きたいとおもいます。ネットで重み付き残差法について調べていると下記のウェブサイトを発見しました。
http://nkl.cc.u-tokyo.ac.jp/12s/intro/
今回はこちらの「有限要素法入門Ⅱ」を参考にして、書かれている問題をpythonにて解きたいと思います。

問題設定

今回は下記のような問題を考えます。

支配方程式 $$\frac{d^2u(x)}{dx^2}+u(x)+x=0 $$ 境界条件 $$u(0)=u(1)=0$$ ちなみに厳密解は下記の通りです。 $$u(x)=\frac{\sin(x)}{\sin1}-x$$

いきなりですが$u(x)$ はどのような関数が良いでしょうか。まず微分方程式があるので二回以上微分可能であれば良さそうです。 また境界条件も考慮する必要があります。以上を踏まえて$u(x)$ は下記のような多項式の関数であると仮定してしまいましょう。 $$u_{app}(x)=x(1-x)(a_1+a_2x)$$ このようにおけば条件を満たせそうです。あとはなるべく良い感じに$a_1$と$a_2$を選べば近似的に解を求めることができるでしょう。
微分方程式に$u_{app}$ を代入すると下記の通りとなり、これを残差$R$とします。 $$R=-a_2x^3+(a_2-a_1)x^2+(a_1-2a_2+1)x+2(a_2-a_1)$$ この残差$R$は$a_1$と$a_2$により変動し、$R$が小さければ小さいほど良い解ということになります。これら$a_1$、$a_2$の求め方として下記のようなものがあります。

• 選点法
• 最小二乗法
• ガラーキン法

以下で順番に見ていきます。

重み付き残差法について

以下では$0<x<1$の範囲で$u(x)$ を考えます。重み付き残差法では$R$に重み関数を乗じたものの積分 を考え、そこから連立方程式を作って解き、$a_1$と$a_2$を算出します。

選点法(Collection Method)

選点法では任意の点($x_i$)で残差を計算し、その点での残差が0となるように係数を算出します。これは重み関数としてデルタ関数を選択していることになります。 $$\int_0^1R(a_1,a_2,x)\delta(x_i)dx=0 $$

最小二乗法(Least Square Method)

最小二乗法では誤差の積分値を考え、その積分値が最小となるように$a_1$と$a_2$の値を算出します。 これは重み関数として$R$の$a_i$に関する偏微分を選択していることになります。 $$2\int_0^1R(a_1,a_2,x)\frac{\partial R}{\partial a_i}dx=0 $$

ガラーキン法(Galerkin Method)

ここで$u_{app}$ は微分方程式の解を以下の形で仮定しています。 $$u_{app}=x(1-x)a_1+x^2(1-x)a_2=\psi_1a_1+\psi_2a_2$$ ガラーキン法では $\psi_i$を重み関数として選択し、$a_1$と$a_2$の値を算出します。 $$\int_0^1R(a_1,a_2,x)\psi_i(x) dx=0 $$

Pythonで解いてみる

実際に手を動かして解いてみます。手を動かすといっても紙と鉛筆ではなく、キーボードです。

import sympy as sym

x,u=sym.symbols('x u')

#Exact Solution 厳密解
u_es=sym.sin(x)/sym.sin(1)-x

#Approximate Solution 近似解
a1,a2=sym.symbols('a1 a2')
u=x*(1-x)*(a1+a2*x)
#Error -> R
R=sym.diff(u,x,x)+u+x

#Collocation Method 選点法
#適当な2点でR=0
eq1=R.subs([(x,1/4)])
eq2=R.subs([(x,1/2)])
#a1,a2について解く
s_cm=sym.solve([eq1,eq2],[a1,a2])
#近似解の算出
u_cm=u.subs([(a1,s_cm[a1]),(a2,s_cm[a2])])

#LeastSquareMethod 最小二乗法
#二乗誤差
Rs=R**2
#二乗誤差の積分値
Err=sym.integrate(Rs,(x,0,1))
#二乗誤差のa1,a2に関する微分
eq3=sym.diff(Err,a1)
eq4=sym.diff(Err,a2)
#近似解の算出
s_lsm=sym.solve([eq3,eq4],[a1,a2])
u_lsm=u.subs([(a1,s_lsm[a1]),(a2,s_lsm[a2])])

#Galerkin Method ガラーキン法
#重み関数の定義
w1=x*(1-x)
w2=x**2*(1-x)
#重み付き積分
eq5=sym.integrate(R*w1,(x,0,1))
eq6=sym.integrate(R*w2,(x,0,1))
#近似解の算出
s_gm=sym.solve([eq5,eq6],[a1,a2])
u_gm=u.subs([(a1,s_gm[a1]),(a2,s_gm[a2])])

#結果のプロット
p=sym.plot(u_es,u_gm,u_lsm,u_gm,(x,0,1),legend=1,show=0)
for idx in range(4):
    p[idx].line_color='C'+str(idx)
    p[idx].label=['es','cm','lsm','gm'][idx]
p.show()

以上の計算を実行すると4通りのグラフが得られます。