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

重み付き残差法について

はじめに

ここでは重み付き残差法について具体例で学んで行きたいとおもいます。ネットで重み付き残差法について調べていると下記のウェブサイトを発見しました。
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通りのグラフが得られます。

【初心者向け北海道スキー場】リフト券とレンタル料金の比較

おはようございます、ヨーサンです。断然インドア派な私ですが、奥さんの誘いもありウィンタースポーツをやろうと計画しています。どうせやるなら「パウダースノーを体験したい!」ということで北海道に行くことに決めました。

しかし北海道には初心者向けとしても数多くのスキー場があり選びきれません。スキー場の質とかはよく分からないので、まずはリフト券やウェアレンタルの値段で比較してみることにしました。

リフト券とレンタル料金の比較結果

スキー場名 リフト券
(一日)
レンタル
(スキー+ウェア)
レンタル内訳
富良野スキー場 5500 9800
(スキー&ウェアセット)
マウントレースイスキー場 5100 9100
(フルセット)
札幌国際スキー場 4600 9600 5300(スキーセット)+4300(ウェア上下・ゴーグル)
キロロスキー場 4500
(ビギナー一日券)
14400 8700(ウェア用具セット)+5700(用具セット)
札幌テイネスキー場 5400 11000 5600(スキーセット)+4400(ウェアセット)
ニセコアンヌプリ国際スキー場 5100 17000
(フルセット)
ニセコビレッジ 6200 10700 5900(スタンダード)+4800(ウェア上下)
ルスツリゾート 6200 13700
(スタンダードスキー)
サホロリゾート 5800 11900
(スキーフルセット)
星野リゾートトマムスキー場 5900 要問合せ

以上です。 一日リフト券に関してはあまり差はありませんが、レンタル料金には結構な差がありますね。

以上スキー場毎でのリフト券とレンタル料金での比較でした。

以下のウェブサイトを参考にさせていただきました。

www.ski-ichiba.jp

四面体節点要素の補間関数

ヨーサンです。テレワーク中は周囲の目が無いのを良いことに勉強に専念しています。研究者ですからね、勉強も仕事です。
最近は有限要素法による電磁界解析について勉強しており、自分の備忘録的にこの記事を書いています。今回は四面体節点要素の補間関数についてです。

 

突然ですが、「四面体節点における何らかの値」だけが分かっていて「その要素内の任意の位置における何らかの値」を知りたい時ってよくありますよね。そうですよねよくありますよね。図で書くとこんな状態です。

f:id:yo444:20201015212833p:plain

日常でよくあるこの疑問

このような時、どうしますか?私はと言うと、私は大雑把な性格なのでもうAをx,y,zの一次式として勝手に仮定しちゃいます。こんな感じです。
A=w0+w1*x+w2*y+w3*z
こうおくと任意の位置でのAを求めることができます。
ただ今度はw0からw3まである4つの謎の値がでてきてしまいました。これは困りました・・・ 
落ち着きましょう。慌ててはいけません。節点におけるAの値はわかっています、四面体の節点の数は4つです。上式の謎の値の数も4つです。お!なにかいけそうな気がしてきますね。この4つ節点における座標と値の関係を書いてみましょう。

A1=w0+w1*x1+w2*y1+w3*z1
A2=w0+w1*x2+w2*y2+w3*z2
A3=w0+w1*x3+w2*y3+w3*z3 
A4=w0+w1*x4+w2*y4+w3*z4

おお、こうすると4つの未知数に対して4つの方程式を作ることができました。この方程式をそれぞれのwについて解くことで任意の位置におけるAの値を求めることができます。めでたしめでたし。

ところでwという変数を、4つの座標とAの値で表すとどうなるのでしょうか。僕は気になりませんが、世の中には気になる人もいるようです。昔の人はよくこれを手計算でやろうと思いましたね。僕なら絶対に間違えます。
それはさておき、これを計算すると任意の位置でのAの値というのは以下のようなスッキリした形で表現できるようです。
A(x,y,z)=N1*A1+N2*A2+N3*A3+N4*A4
よく分からないNという関数を持ち出してきて「スッキリした!」というのはインチキくさいですが勘弁してください。ところでこのNを補間関数と呼びます。Nは(x,y,z)の関数で、書くのがためらわれるくらい複雑な形です(個人の感想です)。この補間関数を導入することでいろいろ良いことがあるようですが、それはまた別のお話です(未だ理解していません)。

という訳で四面体接点要素における補間関数の考え方を示したつもりです。あくまで考え方だけですね、実際の関数の形も示していませんので。しかし補間関数は節点要素だけでなく、辺要素、面要素にも存在するようです。おまけに要素は四面体要素だけでなく、六面体要素もあります。これらはもっともっと複雑です。頑張って学びたいと思います。

続く

 

 

 

 



 

201015ブログを始めます。

はじめまして。ヨーサンと申します。都内某所で研究者として働いており、趣味は筋トレです。

筋トレのこととか、パソコンのことを自分の備忘録的に書いていけたらなと思っています。

よろしくお願いします。