オイラー法とシンプレクティック法を用いた１次元調和振動子の解法

オイラー法とシンプレクティック法を用いた１次元調和振動子の解法

 

１次元調和振動子のハミルトニアン（力学的なエネルギー）は

　　

である。

ここで、

　　

であり、qはバネの自然長からの変位、pは運動量、（１）式のp²/2は運動エネルギー、q²/2はバネの位置エネルギー（ポテンシャルエネルギー）に相当する。

 

（２）式は、

　　

と書き換えることが可能で、（２）式とバネの単振動に関するニュートンの運動方程式と同一のものである。

時刻t=0における初期条件を(q,p)=(1,0)とすれば、（２）、（３）式の解は

　　

である。

 

１次精度のオイラー法を用いて解くには、

　　

つまり、

　　

と近似し、これを逐次繰り返して計算すればよい。

 

たとえば、t=0のとき、(q,p)=(1,0)であれば、Δt=0.1とすれば、

　　

この計算結果をもとに、

　　

といった具合に計算できる。

 

この計算は非常にシンプルな繰り返し計算だから、プログラムを作らずとも、表計算ソフトで簡単に計算が可能で、以下にt=0における初期条件を(q,p)=(1,0)としたものの計算結果を示す。

 

このように非常に粗い近似に基づく計算でありながら、比較的によく計算できていることがわかるだろう。時間の経過とともにqとpの厳密解と１次精度のオイラー法を用いた数値解との食い違いが大きくなるのは、近似計算を繰り返すたびに誤差が伝播し、蓄積するため。このため、オイラー法の誤差は雪だるま式に増加し、計算はやがて発散する。

次に、オイラー法を用い、Δt=0.1とし、t=10まで数値計算結果を(q,p)のグラフとして表したものを以下に示す。

 

解いた問題は、力学的エネルギー、つまり、ハミルトニアンHが一定であり、

　　

(q,p)は単位円上に存在しなければならないのだが、時間の経過とともに軌道q²+p²=1からの逸脱が大きくなってゆく。

つまり、オイラー法を用いた近似計算法はエネルギー保存則を満たす解法ではなないということだ。

 

オイラー法のこの欠点を改良したものにシンプレクティック（Sympretic）法と呼ばれるものがある。

オイラー法で用いる計算式（４）を次のように改良したもの。

 

　　

 

この僅かな改良によってエネルギー保存則からの逸脱がかなり改善され、オイラー法の宿命的な欠点である時間の経過に伴う誤差の伝播、蓄積が解消される。

 

以下に、表計算ソフトを用い、シンプレティック法を用いて解いた計算結果を示す。

 

 

比較参照のために、シンプレクティク法とオイラー法を用いて解いたqと厳密解をのグラフを以下に示す。

 

さらに、シンプレティック法による数値解を(q,p)のグラフで表したものを示す。

 

オイラー法とは異なり、シンプレクティック法を用いた数値解はかなり良くq²+p²=1の円周上に存在していることがわかるだろう。

シンプレクティク法とオイラー法を用いて解いたq,pを元に計算したハミルトニアンHと厳密な値1/2との相対誤差

　　

の時刻による推移をグラフに示す。

 

このグラフが示すように、１次精度のシンプレクティク法は完全にエネルギー保存則を満たすわけでないが、振動をしつつも常にある誤差の範囲内に収まり、オイラー法のように指数関数的に誤差が増大することはない。

参考までに２次精度をのシンプレクティク方による計算結果も示してある。

 

数値計算ではよくあることなのだけれど、（４）式を少し変え（５）式にし計算するだけで計算結果が劇的に変化するのだから、数値計算は奥が深いよね。

そう思いませんか？

そして、少しでもこうした数値計算に興味を持ってもらえれば幸せです。

 

なお、これらの議論をより正確にするためには、ニュートン力学の進化形である解析力学などの知識が必要になるので、感覚的に理解できれば十分ではないでしょうか。

 

参考までに２次精度のシンプレクティック法で解いた計算結果は以下に示す。

コメント 12

　ddt^3です。またご迷惑をおかけします。

［変分法-1］
　自分は有限要素法（FEM）の前に境界要素法（BEM）をやっていたという、変な奴です(^^;)。
　ちなみにFEMは、Finite Element Methodと書いて、ファイナイト・エレメント・メソッドと読みます（フィニットと読んじゃだめ）。BEMは、Boundory Element Methodと書いて、バウンダリー・エレメント・メソッドと読みます。

　BEMにもFEMにも、色んな定式化がありますが、やっぱり変分原理の成り立つケースが、一番綺麗に行きます。変分原理とは何かと言うと、普通は最小エネルギー原理と等価になります。では最小エネルギー原理とは何かと言うと、ポテンシャルエネルギーの定義から明らかだったりするんです。

　バネ定数kのバネが天井に鉛直に固定され、先端に質量mの質点がぶら下がってる状態を想像します。重力をmgとして、バネと質点のポテンシャルエネルギーは、次のように計算しますよね？。釣り合い方程式は、
　　kx＝mg
なので（xはバネの伸び）、上式を伸びxで積分して、
　　∫kx dx－∫mg dx＝E0
　ここでE0は積分定数です。上記積分は簡単にできて、
　　1/2×kx^2－mgx＝E0
です。E0は例えば、x＝0の地点を基準点に選べば0にでき、質点の位置エネルギーの減少分がバネの弾性エネルギーとして貯まった、と読めるようになります。
　次に、
　　E(x)＝1/2×kx^2－mgx
という関数を考え、E(x)の最小点を探してみます。とうぜんE(x)をxで微分してdE/dx＝0としますよね？。
　その結果は、
　　kx－mg＝0
となり、釣り合い方程式に戻れます。バネはx＝mg/k伸びて釣り合う、です。
　でもこれって、当たり前じゃないですか？。だってE(x)の関数形は、kx－mg＝0を積分して作ったんですから。本質的には変分原理とはこれだけです。

　　釣り合い方程式（系の支配方程式） ⇔　エネルギーの最小化（変分原理）

　一般の場合は今のような考えを、3次元の弾性体や運動系にも適用できるように拡張しただけのものです。運動系の場合、支配方程式は運動方程式です。厳めしく変分原理などと言うと難しそうに聞こえますが、エネルギーの定義まで遡って考えれば、「当然じゃん」という事が案外多いです。
　※ 運動系のラグラジアンもよぉ～く考えると、仮想静止系のエネルギーだったりします(^^)。

　そういう訳で、系を解析する方法は2つある事になります。支配方程式（釣り合い方程式や運動方程式）を直接扱うか、エネルギーの最小化（変分原理）という形で陰的に解を求めるか。もちろん変分原理が不便だったら誰も使いません。変分原理による数値解法の方が、精度が良かったりコンピュータとの相性が良かったりするケースがけっこうあります。もちろん支配方程式を直接扱う方が、わかりやすいですが。
by ddtddtddt (2017-08-04 11:01) 

［変分法-2］
　まず具体的な例を紹介します。非圧縮性渦なし完全流体のケースです。このケースでは、速度ポテンシャルφ(x，y)を計算できれば良い事になります（面倒なので2次元にします）。
　流体の任意点の速度ベクトルは(∂φ/∂x，∂φ/∂y)で算定でき、およそ流体力学で知りたい事の全ては、流体各点の速度だからです。このとき最小化するべきエネルギーとして用いるのは、もちろん流体の全運動エネルギーです。
　　J＝∫1/2×ρ((∂φ/∂x)^2＋(∂φ/∂y)^2) dxdy　　　(1)
　ρは流体密度で、上記はちゃんと運動エネルギーの形をしてますよね？。質点の場合との違いは質量mが、密度ρ（単位体積当たりの質量）になってるところです。積分領域R内の微小面積dxdyがちっちゃいちっちゃい質点ρdxdyです。重さを密度で表したので、∫dxdyで積分しますよ～というだけです。積分領域Rは、解析したい流体の全体積（面積）になります。

　ここで前回のxに相当するものはφです。φが系の挙動を決めるからです。そうするとエネルギーJの最小化とは、Jを最小にするような2変数関数φ(x，y)を求めよ、という話になります。この一見解くのが不可能そうな問題は、ゴールドスタインの古典力学（吉岡書店）の方法に従うと割と見通し良く解決できます。もちろん数学的には厳密ではありませんが(^^;)、読んだ中では一番わかりやすかったです。

　まずベタに考えます。可能な全てのφ(x，y)に総当たりして(1)を計算し、Jが最小になった時のφを選べば良い訳ですよね？。しかし可能な全てのφ(x，y)とは、偏微分可能な任意の2変数関数全部の事です。それらは、ほぼ全ての連続関数くらいたくさん（無数に）あり総当りはとても無理ですけれど、可能な全てのφ(x，y)に番号付けする事は、数学的理想化としてはありです。例えばφ(x，y，0)，φ(x，y，1)，φ(x，y，2)，・・・みたいに(^^)。
　でも全ての連続関数くらいたくさんあるので、自然数n＝0，1，2，・・・では番号の数が不足かも知れません。そこで実数αで番号付けする事を考えます。φ(x，y)をφ(x，y，α)とみなすという話ですが、これをテスト関数と呼びます。そうすると(1)は、
　　J(α)＝1/2×ρ∫((∂φ(x，y，α)/∂x)^2＋(∂φ(x，y，α)/∂y)^2)dxdy　　　(2)
です。
　ここで(2)の入出力関係を整理します。流体密度ρは一定（定数）です。∫dxdyの積分領域Rは固定（定数みたいなもの）です。よって(2)右辺の変数は未知関数φ(x，y，α)だけですが、φの変化は番号付けパラメータαが担います。従って積分値Jは、実数αのみの関数J(α)です（← OKですか？(^^;)）。
　次にパラメータαは、人間の方で勝手に設定するものです。その設定はかなり勝手にできます。例えば(2)のJ(α)の最小は、α＝0で起きると決められます。つまり解φ(x，y)を、φ(x，y)＝φ(x，y，0)で表すという事です（任意に決められるから、本当はα≠0でも良い）。
　重要なのは次の想定です。

　せっかくαで番号付けしたんだから、|α|が小さいほどテスト関数φ(x，y，α)は解φ(x，y)＝φ(x，y，0)に近い、として良いだろう。

　・・・そして(^^;)、

　α→0の時のφ(x，y，α)のφ(x，y)への近づき方は、J(α)がα＝0で微分可能になるくらい滑らかに設定できるはずだ。

　・・・と(^^;)。

　以上を認めれば、やる事は一つです。「J(α)の最小を知りたいなら、J(α)をαで微分して0とおけ」。
　　dJ/dα×dα＝0（dJ/dα＝0）　　　(3)

　という訳で、「あれぇ～？。不可能そうな話が、普通の微積になっちゃった(^^)」です。

　(3)はいわばJ(α)の全微分を取った形です。全体が0になるのはdJ/dα＝0だからです。そうすると×dαは余計に見えますが、これを残しておかないと以後の話が続きません。そこは普通の微積とちょっとは違います(^^;)。それに(3)では、まだ絵に描いた餅です。

　(3)の形に対応する(2)の右辺を想像します。(3)はαによるJ(0)からのJ(α)の値のずれを表し、dαは0からのαの微小なずれです。よって微小なdαに対応する、(2)右辺のずれが具体的にわかれば良い訳です。
　　φ(x，y，α)－φ(x，y)＝δφ(x，y，α)　　　(4)
でδφを定義し、変分と呼びます。つまりδφはαのずれに対応するφのずれで、δφは立派な一つの2変数関数です。
　いまα＝dαは微小なので、先の想定からδφも微小です。従って(2)右辺の積分のずれは、δφ(x，y，α)を各点(x，y)ごとで考えて、
　　dJ/dα×dα＝1/2×ρ∫(∂φ/∂x×δ(∂φ/∂x)＋∂φ/∂y×δ(∂φ/∂x))dxdy＝0　　　(5)
　　※ δφが微小という条件しか使っていません(^^)。

になります。
　けっきょく最後にはδφを各点(x，y)ごとに考えるんだから、別にδφをdφと書いても本当は良かったんですよ。しかしそれではφの全微分と混同しそうなので「関数のずれを表してるんだよぉ～」と、慣習的にδφが使われます。それにδでないと、(5)の中辺が妙な感じになるし(^^;)。
　意味から（定義(4)に従えば）、
　　δ(∂φ/∂x)＝∂(δφ)/∂x，δ(∂φ/∂y)＝∂(δφ)/∂y
は明らかなので、
　　∫(∂φ/∂x×∂(δφ)/∂x＋∂φ/∂y×∂(δφ)/∂y)dxdy＝0　　　(6)
を得ます。∫dxdyの積分領域Rは固定で、(6)は「定積分」です。

　ここで具体的計算には役に立たないけれど注意すべきは、(5)の積分が0になる根拠こそ、(3)だという事実です。この事実は定式化上は重要なんですが、具体的計算には役立たずなので、ここで(3)はお役御免です(^^;)。
by ddtddtddt (2017-08-04 11:09) 

［変分法-3］
　　∫(∂φ/∂x×∂(δφ)/∂x＋∂φ/∂y×∂(δφ)/∂y)dxdy＝0　　　(1)
　∫dxdyの積分領域Rは固定で、(1)は「定積分」でした。

　変分法の常套手段として、(1)を部分積分し変分δφの微分を消しますが、その目的は後述します。(1)の場合は、次のグリーンの定理を使います。グリーンの定理とは、ガウスの発散定理の部分積分版です。
　　∫(∇f)・(∇g) dxdy＝∫g (∇f)・ds－∫g Δf dxdy　　　(2)
　∇は勾配で、∫dxdyは固定領域Rで行います。dsはRの境界Cの外法線線素ベクトルで、・は内積です。∫・dsはC上の線積分、Δはラプラシアンです。

　(2)でf＝φ，g＝δφとおけば、確かに、
　　∫(∇f)・(∇g) dxdy
　＝∫(∂φ/∂x，∂φ/∂y)・(∂(δφ)/∂x，∂(δφ)/∂y)dxdy
　＝∫(∂φ/∂x×∂(δφ)/∂x＋∂φ/∂y×∂(δφ)/∂y)dxdy
なので、(1)は、
　　∫δφ (∇φ)・ds－∫δφΔφ dxdy＝0　　　(3)
となります。
　δφは、(3)を満たす解関数φ(x，y)からのテスト関数のずれを表す、微小だけれど立派な2変数関数でした。可能なテスト関数とは、任意の偏微分可能な関数です。よってδφは、微小な任意の偏微分可能な関数を表してる事になります。ここで面倒なので、偏微分可能なら連続だろうからδφの範囲を、任意の微小な連続関数とします(^^;)。

　最初に(3)の1項目です。このまま続行すると最後には、非圧縮性渦なし完全流体の支配方程式を表す偏微分方程式にたどり着くのですが、偏微分方程式の解は一般に、解析領域Rの境界C上の境界条件を与えて初めて、解が一意に決まります。
　どんな境界条件かは知りませんが、今はとにかくどんな境界条件にも適用可能なφの条件が知りたい訳ですから、φはとにかく境界条件を満たすとするのが妥当です。すなわちRの境界C上で、ずれδφ＝0です。そして(3)の1項目はC上の線積分でした。δφ＝0から、それは0とおくのが妥当です。
　　∫δφΔφ dxdy＝0　　　(4)

　ここで(4)が「不定積分」だったら明らかに、Δφ＝0を導けます。しかし(4)は「定積分」でした。
　そこでもう一回、δφの範囲を思い出します。δφは任意の微小な連続関数でした。任意の微小な連続関数って、どうやって作ったら良いでしょう？。
　最も乱暴には、「任意の連続関数fを勝手に一つ持ってきて、それに任意に勝手で微小な定数εをかければOKよね？」という方法があります（← 間違ってはいない(^^;)）。
　f(x，y)を任意の2変数連続関数として、δφ＝εf(x，y)とします。
　　∫εf Δφ dxdy＝0

　でもεは定数で積分とは無関係です。従って、任意の連続関数fに対して、
　　∫f Δφ dxdy＝0　　　(5)
でなければならない事になります。

　ここで［変分学の基本定理］をご紹介します。2次元でも3次元でも同じなので、1次元で紹介します。
［変分学の基本定理］
　f(x)とg(x)を連続関数とする。任意のf(x)について、
　　∫f(x)g(x)dx＝0
である条件は、g(x)＝0。

　この定理は約半世紀前には本当に、［変分学の基本定理］と言われてたようです。しかし［基本定理］という仰々しい名前に対して証明が余りにこっ恥ずかしいからなのか、いつの間にか明確に語られなくなりました（← 自分は非常に問題な気がします）。実際うるさい事さえ言わなければ上記定理は、馬鹿みたいに簡単に証明できます。

　fは任意だから、f＝gの場合がある。この時は、
　　∫f(x)g(x)dx＝∫(g(x))^2 dx＝0

　非負の連続関数の積分が（0以上の値を持つグラフの面積が）0になるなら、g(x)＝0は明らか。逆にg(x)＝0なら、与式は明らか。

　・・・という訳で、(5)よりfが任意の連続関数なので（Δφもきっと連続だろうから）、
　　Δφ＝0
となり、非圧縮性渦なし完全流体の支配方程式が得られます。

　もうおわかりだと思いますが、「変分法の常套手段として、(1)を部分積分し変分δφの微分を消す」目的は、［変分学の基本定理］を使いたいがためなんですよ。

　［変分学の基本定理］の証明を厳密にやろうと思ったら、「非負の連続関数の積分が0になるなら、g(x)＝0」をちゃんとやる必要があり、正式には測度論まで遡らないと駄目です。それが無理な注文である事はわかりますが、だったら「こんなの明らかよね？」って言えば良いじゃないですか。そうせずに部分積分の目的も明確に言わないくらいなら。いかに「こっ恥ずかし」かろうと・・・。
by ddtddtddt (2017-08-04 11:18) 

［変分法-4］
　・・・こうして数々の数学的いい加減さを闇に葬って来ました(^^;)。中でも目玉級に問題なのは、テスト関数の番号付けパラメータ、αの存在でしょう。

　じつはこれ、関数空間上の距離に相当します。関数空間とは乱暴に言うと、2変数の連続関数全体の集合なんかですが、その要素である関数の間に距離を定義できる事がわかっています。距離の定義も色々ありますけれど、どのような距離を取ろうと本質的には同等です。特に2つの関数f，gについての距離をα(f，g)として、α(f，g)＝0　⇔　f＝g　はどんな距離においても成り立ちます。
　関数空間に距離αを入れて距離空間とし、その上でJ(α)について通常の微積をやるというのが、変分法の正式なやり方です。これは実数上の微積分だって、じつは同じです。xとyを実数とすれば実数全体の集合には最初から、距離α＝|x－y|が備え付けになってるので普段は意識しないだけです。

　普通は関数空間までやる余裕はなく、そこは省略されます。それでほとんど変分の説明もないままに、ほとんど突然に「こんなの明らかだよね？」と、
　　∫(∂φ/∂x×∂(δφ)/∂x＋∂φ/∂y×∂(δφ)/∂y)dxdy＝0　　　(1)
などが、ドド～ン！と現れます(^^;)。

　また［変分学の基本定理］が明確に紹介される事もなく、従って(1)を、
　　∫δφ (∇φ)・ds－∫δφΔφ dxdy＝0　　　(2)
の形に変形する目的も暗中模索のまま、いつのまにか訳もわからずに、
　　Δφ＝0　　　(3)
へ達する事になります。普通の神経では耐えられませんな(^^;)。

　ともあれ変分法でΔφ＝0を導出できたので（←ホントか？(^^;)）、このまま境界要素法へ突入して良いでしょうか？（ネコ様へ）。
by ddtddtddt (2017-08-04 11:24) 

記事の提供、感謝いたしますm(__)m

☆このまま境界要素法へ突入して良いでしょうか？（ネコ様へ）。
◇はい、結構です。
ですが、私、８月６日からお盆明けくらいまで雲隠れしますので、この点はよろしくお願いします。


by nemurineko (2017-08-04 18:53) 

　8月11日から私もサマーゲレンデに行ってきました（山籠り？(^^)）。

［境界要素法-1］
　線形問題に限って言えば、境界要素法は恐らく最強です。精度に優れ計算量も少なく、任意点の結果を取りやすい。計算プログラムの構成も明解でシンプル。
　難点は特異積分（広義積分）を扱う必要のある点ですが、今では解析的積分公式やさまざまな計算テクニックも開発されてるので、やりゃ～何とかなります。

　境界要素法（BEM）の原理は有限要素法という計算スタイルより先行し、グリーン関数法です。
　　http://nekodamashi-math.blog.so-net.ne.jp/archive/c2305704962-1　　　(a)
　　の［コーシーの積分公式の周辺-1］
　グリーン関数法は条件さえ合えば確かに便利な手段でしたが、簡単な境界条件でしか便利にならないという計算法としては致命的な所がありました。それでいかなる境界条件にも理論上は簡単に適用可能な、変分原理に基づいた有限要素法（FEM）が先に発展します。そしてFEMがそろそろ常識化しかけた頃、という事は「解析領域を要素に区切って数値計算」という計算スタイルに大抵の人が慣れてきた頃、一部のFEM人がラプラス方程式などをそういう目で見直し、それらも「要素に区切ってやっつけられる」事に気づきます（恐らくブレビアあたり）。1970～80年代の出来事です。で事後検証してみると、「要素に区切ってやっつける定式化」は、グリーン関数法の現代風アレンジになっていたぁ～、という次第です(^^)。

　変分原理から導けるのはラプラス方程式だけではありません。ポアソン方程式も導けます。
　　J＝∫(1/2×ρ((∂φ/∂x)^2＋(∂φ/∂y)^2)＋ρφ×g(x，y))dxdy　　　(1)
としてJを最小化すれば（やり方はこの前といっしょ）、ポアソン方程式、
　　Δφ＝g(x，y)　　　(2)
が得られます。
　(2)は湧き出し密度分布g(x，y)を持つ、非圧縮性渦なし完全流体の支配方程式で、g(x，y)は既知関数です。Jの意味は、流体の運動エネルギーと湧き出しのポテンシャルエネルギーも考慮した、流体全体のエネルギーを表します。
　φが静電ポテンシャルなら、g(x，y)は与えられた電荷密度分布。Jは電荷のエネルギーも考慮した静電場エネルギーです。
　グリーンの定理に戻ります。グリーンの定理はグリーンさんが特にグリーン関数法のために、ストークスさん供に導いたものです。
　　∫(∇φ)・(∇δφ) dxdy＝∫δφ (∇φ)・ds－∫δφΔφ dxdy　　　(3)
　(3)を＝0とおけばラプラス方程式（やポアソン方程式）が得られるのでした。従って(3)とともにラプラスやポアソン方程式を与えれば、変分を取ったのと同等になります（本当は要証明ですが(^^;)）。

　ところで変分δφは実質、連続関数であれば何でもOKでした。それでδφとして、φと同じタイプの方程式を満たすものを考えます。問題としては、偏微分方程式(2)を解きたいとします。
　変分δφは、(2)のタイプとして特に次の条件を満たすとします。
　　Δφ'＝δ(ξ，η)　　　(4)
　δ(ξ，η)は点(ξ，η)に特異点を持つデルタ関数で、φ'(x，y，ξ，η)をラプラス方程式の基本解と呼ぶ事が多いです。ここでφ'が(x，y)だけでなく(ξ，η)の関数にもなってるのは、(4)からφ'はδ(ξ，η)で決定されるので、φ'はきっとその特異点(ξ，η)もパラメータとして含むよね、という意味です。

　形式的にφ'を(3)に代入します。
　　∫(∇φ)・(∇φ') dxdy＝∫φ'(∇φ)・ds－∫φ'Δφ dxdy　　　(5)
　で、φ'はポアソン方程式を満たすのでした。という事は今度はφ'に対してφを、その変分だと考える事も可能です。次式は直接計算しても出せます。
　　∫(∇φ')・(∇φ) dxdy＝∫φ(∇φ')・ds－∫φΔφ' dxdy　　　(6)

　∫(∇φ)・(∇φ') dxdy＝∫(∇φ')・(∇φ) dxdyは明らかです（・は内積）。(5)と(6)の辺々引くと、
　　∫φ'(∇φ)・ds－∫φ(∇φ')・ds－∫φ'Δφ dxdy＋∫φΔφ' dxdy＝0
を得ます。
　ここでφは(2)を満たし、φ'は(4)を満たします。
　　∫φ'(∇φ)・ds－∫φ(∇φ')・ds－∫φ'g(x，y) dxdy＋∫φδ(ξ，η) dxdy＝0
　いま(ξ，η)が積分領域Rの内点なら、デルタ関数の性質から、
　　∫φδ(ξ，η) dxdy＝φ(ξ，η)
です。従って先の関係式から、
　　φ(ξ，η)＝∫φ(∇φ')・ds－∫φ'(∇φ)・ds＋∫φ'g(x，y) dxdy
という事になります。
　ここで、もう少し表現を簡単にします。(∇φ)・dsは、Rの境界Cの線素をdcとした時、方向微分の公式から(∇φ)・ds＝∂φ/∂n×dcになるのがわかります。∂φ/∂nは、C上でのφの外法線方向微分値です。これをqで表しますが、(∇φ')・dsについても同様です。q，q'は、境界上の流速値とか外法線微分とか呼ばれる事が多いです。

　　φ(ξ，η)＝∫φq'dc－∫φ'q dc＋∫φ'g(x，y) dxdy　　　(7)

　(ξ，η)は解析領域Rの内点であれば、どこでもOKでした。よってφ(ξ，η)は、Rの任意の内点(ξ，η)における解関数φ(x，y)の値です。それで(7)を、境界要素法の内点方程式と呼びます。
　(7)の関係式自体は、グリーン関数法の開発過程で相反定理の名のもとに知られていました。当時はデルタ関数の方法がなかったので、こんなに簡単には導けませんでしたが。
　じっさい基本解φ'がφと同じ境界条件を満たすとすれば、これはRの境界C上でφ＝φ'，q＝q'という事なので、(7)のC上の境界積分項は打ち消しあって0となり、参考URL(a)にあげたオリジナルなグリーン関数法の計算式と同等なものになります（グリーン関数の対称性）。
　しかし(4)を満たし、その上φと同じ境界条件を満たすφ'をみつける計算の方が、(2)を解くよりよっぽど難しいという事態はよく起こります。そこがグリーン関数法の実用上の問題でした。

　(7)は、変分関係式(3)に支配方程式(2)を持ちこんだものです。もし変分δφが任意であれば、それは変分を取ったのと同等であり、(7)から陰的に解関数φを計算できる可能性が開けます。しかし(7)では、(4)を満たすという特定の変分φ'しか用いていません。ところが条件(4)があればφ'は十分に任意化されていると、FEMになれた一部の計算技術者や理論家達は気づいたのでした。
by ddtddtddt (2017-08-18 08:56) 

［境界要素法-2］
　境界要素法の内点方程式。
　　φ(ξ，η)＝∫φq'dc－∫φ'q dc＋∫φ'g(x，y) dxdy　　　(1)
　∫dxdyは解析領域R上の積分，∫dcはRの境界C上の線積分，φはΔφ＝g(x，y)を満たす未知関数です。qとq'はφとφ'の外法線微分値。
　φ'はδをデルタ関数として、
　　Δφ'＝δ(ξ，η)　　　(2)
を満たすものなら、何でもOKでした。そこで(2)を満たすφ'(x，y，ξ，η)として出来るだけ簡単なものを考えます。これは実際に計算可能になる事が多いです。例えばφ'として(ξ，η)を中心に等方的なものを選び、境界条件は考えないとすれば、φ'(x，y，ξ，η)＝1/(2π)×log(r)，r＝√((x－ξ)^2＋(y－η)^2)　が容易に得られます。
　　http://nekodamashi-math.blog.so-net.ne.jp/archive/c2305704962-1　　　(a)
　　の［コーシーの積分公式の周辺-1］
　(1)の入出力関係を整理します。φ'は具体的に決められたので、右辺で ' つきのものは全て既知関数です。領域積分項∫φ'g(x，y) dxdyはg(x，y)も既知関数なので、具体的に計算可能です。未知なのは、境界積分項に現れるφとqだけです。そうすると境界C上のφとqさえわかれば、(1)から未知関数φ(x，y)を計算できる事になります。

　解析領域Rの境界Cを折れ線近似し、折れ線の角に節点jを配置した絵を想像して下さい。jは節点番号で、左回りにCを一周します。節点jとj＋1を端点として持つ線分を境界要素と呼び、それに要素番号kを与えます（左回りにCを一周）。節点jでのφとqの値は、φjとqjで表します。
　境界要素kの要素長Lkが十分に小さければ、k上の関数値：φ(c)，q(c)，0≦c≦Lk　は、(φjとφj+1)および(qjとqj+1)で線形近似でもすれば、十分なはずです。t＝c/Lk，0≦t≦1として、
　　要素k上：　φ(c)＝(－t＋1)φj＋tφj+1，t＝c/Lk
　　要素k上：　q(c)＝(－t＋1)qj＋tqj+1，t＝c/Lk

　一方、特異点(ξ，η)＝(ξi，ηi)によるφ'とq'を、φ'i，q'iで表します。iは特異点の場所を識別する番号で、i＝1，2，・・・で十分です。
　これらを(1)の境界積分項に代入し、要素k上で考えてやれば、

　　∫[k]φq'dc
　　＝φj×∫[k] q'i×(－t＋1) dc ＋φj+1×∫[k] q'i×t dc

　　∫[k]φ'qdc
　　＝qj×∫[k] φ'i×(－t＋1) dc ＋qj+1×∫[k] φ'i×t dc

となります。
　上記右辺の積分において、q'iとφ'iは既知関数なので、要素k上で具体的にq'i(t)，φ'i(t)の形に書けます。tは積分パラメータ（0≦t≦1)です。
　従って右辺の積分は、具体的な数値Hij（Hij+1)とBij（Bij+1)になります。q'iに関するφjについての積分をHij，φ'iに関するqjについての積分をBijで表しました。
　境界積分全体の値は、これらの集計です。節点番号jについて和をとる事になります。
　　φi＝Σ[j] Hij×φj－Σ[j] Bij×qj＋Wi
ただしφiは、
　　φi＝φ(ξi，ηi)
Wiは、
　　Wi＝∫φ'(x，y，ξi，ηi) g(x，y) dxdy
です。

　ここで行列記法を思い出して下さい。i＝1，2，・・・だから、
　　Σ[j] Hij×φjかつi＝1，2，・・・は、(Hij)(φj)と書ける。
　　Σ[j] Bij×qjかつi＝1，2，・・・は、(Bij)(qj)と書ける。
・・・ですよね（← OKですか？(^^;)）。

　　(φi)＝(Hij)(φj)－(Bij)(qj)＋(Wi)　　　(3)

となります。
　(3)で(Wi)は、具体的に計算可能な既知量でした。よって(3)は、未知量(φj)と(qj)に関する連立一次方程式の形に、ほぼなっています。係数行列は(Hij)と(Bij)です。もし(3)を解いて(φj)と(qj)を定められれば、線形近似の形で境界上のφとqが得られた事になり、内点方程式(1)から未知関数φ(x，y)の近似解を求められます。左辺の(φi)さえ何とかなれば・・・。
　φiは、
　　φi＝φ(ξi，ηi)＝∫φδ(ξi，ηi)dxdy
なのでした。でもこれは、「(ξi，ηi)が積分領域Rの内点の場合には」なのでした。つまり(ξi，ηi)は、Rの外にあっても（外点でも）OKです。その時は、デルタ関数の性質から、
　　∫φδ(ξi，ηi)dxdy＝0
です(^^)。そういう訳で、

　　－(Hij)(φj)＋(Bij)(qj)＝(Wi)　　　(4)

によって、境界未知量(φj)，(qj)を計算できる事になり、(1)を併用します。
　実際には境界条件として与えられた既知量を(4)に代入し、式の数を調整します。典型的には、境界上のφの値全部かqの値全部が与えられます。理屈の上では特異点(ξi，ηi)は、必要な数だけR外部の任意の位置に、自由に設定すれば良い事になります。

　(4)には境界未知量しか出てきません。そこで(4)を境界方程式と呼びます。(4)の理論的形はグリーン関数法のいわば特殊解法として、昔から知られていたものでもありましたが、FEMになれた人達が気づいたのは、次の点でした。

　φi，i＝1，2，・・・　において(ξi，ηi)≠(ξk，ηk)とすると、(φ'i)は（(q'i)も）線形独立な関数系を構成します。従って、i＝1，2，・・・によって与えられた(4)は、独立な線形条件を成すはずだと。
　またいかに連続関数であろうと、離散化して考えれば節点自由度の有限の自由度しか持たないと。だとすれば、変分δφ＝φ'は完全に任意である必要はなく、i＝1，2，・・・で十分に任意化されている、と。

　このような発想が手軽に出来るようになったのは、やはり関数空間以後ですが、それだけでは理論止まりだったと思われます。「要素に区切ってやっつける」計算スタイルがFEMにより具体的に普及した事が、最後の決定打になった気がします。

　以上のやり方は理論的にも具体的な計算方法としても、最も明解なものです。ところが境界要素法の中で(4)は間接法と言われ、じつはマイナーなんです。というのは特異点の外部配置が余りにも自由過ぎて、数値計算の専門家達にも正確な誤差評価が無理だったからです。
　自分は、実用的な数値解が得られりゃOKよ程度のユーザーなので、(4)も時々使いますが、大概は専門家達に不評な方法なら、趨勢は右にならえです。

　という訳で、境界要素法の直接法という話になります。
by ddtddtddt (2017-08-18 09:12) 

［境界要素法-3］
　境界要素法の内点方程式。
　　φ(ξ，η)＝∫φq'dc－∫φ'q dc＋∫φ'g(x，y) dxdy　　　(1)
　∫dxdyは解析領域R上の積分，∫dcはRの境界C上の線積分，φはΔφ＝g(x，y)を満たす未知関数です。qとq'はφとφ'の外法線微分値。φ'はδをデルタ関数として、
　　Δφ'＝δ(ξ，η)　　　(2)
を満たし、φ'(r)＝1/(2π)×log(r)，r＝√((x－ξ)^2＋(y－η)^2)　を取れます。
　以上が出発点です。

　φ(ξ，η)は、
　　∫φδ(ξ，η)dxdy
の積分結果なのでした。(ξ，η)を積分領域Rの内部に置いた場合（内点とした場合）、(1)が得られますが、解けない形でした。(ξ，η)を積分領域Rの外部に置いた場合（外点とした場合）は、
　　∫φq'dc－∫φ'q dc＋∫φ'g(x，y) dxdy＝0　　　(3)
となり離散化すれば解けますが、不評なのでした。とすればもう後できる事は、(ξ，η)をRの境界Cに置く（境界点とする）しかないじゃないですか(^^;)。
　デルタ関数δ(ξ，η)を境界上で考えRで積分する場合、∫δ(ξ，η)dxdyの再評価が必要になります。評価結果は簡単ですが、ここでたいてい怪しげな計算が登場し（コーシーの主値積分）、一部の理論家からは酷評される事になります。デルタ関数の数学的に厳密な定義にまで遡って再評価するなんて、やってられないですもんね(^^;)。そこでここでは、どうせいい加減になるならばという事で、∫δ(ξ，η)dxdyの再評価にデルタ関数の等方性を使います。

　デルタ関数δ(ξ，η)って、特異点(ξ，η)を中心に等方的ですよね？（本当は関数じゃないけど）。だって(x，y)≠(ξ，η)ではδ＝0で、(x，y)＝(ξ，η)ではδ＝∞なんですから、これを等方的と言わずして何と言う？です(^^)。
　次に(x，y)≠(ξ，η)ではδ＝0ですから、(ξ，η)を内点として含む限りどんな積分領域を取っても、積分結果は同じです。なので特に積分領域として、(ξ，η)を中心とした半径εの円を取れます。(ξ，η)を中心とした円は、(ξ，η)を中心に等方的です。
　等方的な関数を等方的な領域で積分したら、どうなりますか？。例えばr＝√((x－ξ)^2＋(y－η)^2)を、(ξ，η)を中心とした円で積分したら。もし∫r dxdy＝1になったとしたら、半円での積分値は1/2ですよね？。1/4の扇型に積分領域を制限したら、明らかに積分値は1/4ですよね？。

　δ(ξ，η)を境界上に置いた場合、半径εの円が十分小さければ(ξ，η)の近傍で、積分領域Rの境界Cはふつう直線とみなせます。半径εの円は、その直線によって半分だけRに引っかかります。よって、
　　∫δ(ξ，η)dxdy＝1/2
です。εがいくら小さくても直線とみなせないケースもあります。(ξ，η)がCの角点になるケースです。この時は角の内角をkとすれば明らかに、
　　∫δ(ξ，η)dxdy＝k/(2π)
です。従って、∫φδ(ξ，η)dxdy＝k/(2π)×φ(ξ，η)です。デルタ関数はこういう事が、普通の関数と同じく実用的に出来るように、非常に注意深く造られたものです(^^)。

　前回やったように(1)を離散化します。結果は、
　　(φi)＝(Hij)(φj)－(Bij)(qj)＋(Wi)　　　(4)
でした。ここで(Hij)，(Bij)，(Wi)は具体的な数値で与えられます。iは特異点番号，jは節点番号です。
　いま特異点は境界上にあるので、iはどれかのjと一致します。一致するj＝iの内角をkiとすると、(4)左辺のφiは、ki/(2π)×φiにおきかわる事になります。i＝1，2，・・・を考慮して行列記法で書けば、
　　(ki/(2π)×dij－Hij)(φj)＋(Bij)(qj)＝(Wi)　　　(5)
です。iが角でない場合ki＝πなのでki/(2π)＝1/2。dijはクロネッカーのデルタ。
　(5)の未知数も境界未知量のみなので、これも境界方程式と呼ばれます。(5)を解いて境界上のφ，qを定め(1)を併用するのが、境界要素法の直接法です。

　(5)から明らかなように、直接法では自動的に未知数の数に等しい条件数が得られます。また(5)の理論的誤差は、境界Cを折れ線近似した境界要素の長さと近似関数で決まるので、誤差評価も容易になります。
　こうしてユーザーにはよりわかりやすく（← ホントか？(^^;)）、理論家達も満足する定式化が得られました。

　しかし代償もあります。基本解φ'(r)＝1/(2π)×log(r)，r＝√((x－ξ)^2＋(y－η)^2)　はr＝0に特異性を持ちます。間接法の場合、φ'の特異点(ξ，η)は積分領域Rの外部にあるので、(1)右辺の積分は普通に行えます。直接法の場合(ξ，η)は境界上にあるので、r＝0を含む特異積分を処理する必要に迫られます。

　φ'の特異性はlog(r)のオーダーなので可積です。q'の特異性はアークタンジェントの計算に帰着できます。ただしそのような計算はコンピューターには出来ないので、特異点を持つ境界要素の積分では、人間の指示で場合分けし、解析的積分公式をプログラムに書いてやるか、特殊要素を用いる事になります。

　でもやれば、何とかなりますよ(^^)。
by ddtddtddt (2017-08-18 09:24) 

ddt^³さん、こんにちは。

境界要素法の記事をブログの記事にアップしましたが、式の書き間違いの指摘や、ここはこうして欲しいというご要望がありましたら、申し付けください。できる範囲で、応えたいと思います。

なお、投稿記事中のφ’、ｑ’などの記号は、読者が微分と勘違いするおそれがありますので、φ^＊やq^＊に変更しました。

これは私からの要望ですが、BEMを用いたラプラス方程式（ポアソン方程式）の解法の簡単なサンプルプログラムをお示しくだされば、読者の理解はより深まると思います。
汎用的な、高級なものは理解を妨げるだけなので、２次元平面上に１つ、２つ点電荷を置いた静電場のポテンシャルを求めるとか、もっと簡単な極座標で表された定常熱伝導方程式
ΔT=0
を解く、これらの問題だけを解くことに特化したプログラムとか・・・。
使用言語は、FortranやC、または（十進）BASIC。
行列の連立方程式を解くのは、シンプルなガウス消去法で十分です。帯行列の連立方程式を解くとか、誤差を小さくするために、連立方程式を並び替えるとか、そんな高級な処置を施す必要はありませ。最もシンプルなガウス消去法や掃き出し法などで十分です(^^)


by nemurineko (2017-08-19 10:45) 

　ddt^3です。ネコ様ありがとうございます。
　要望というか、ほんとんど元ネタの間違いです。

　「境界要素法入門 part 1」の［変分法-2］の式(5)（2つある）の両方で、右辺の1/2を消して下さい。式のコピペでやっちまった不注意です。

　最初の式(5)の後の「・・・(5)の中辺が妙な感じに・・・」の「中辺」を「右辺」に。
　後の式(5)の式番号を(6)にして、右辺の後に＝0を足して欲しいです。

　「境界要素法入門 part ２」の［変分法-3］の冒頭の式(1)に、右辺として＝0を加えて下さい。

　以上、勝手な事を書きますがお願いします。


＞・・・BEMを用いたラプラス方程式（ポアソン方程式）の解法の簡単なサンプルプログラム・・・

となりますと、図や数式をちゃんと作らないと、と思うのですが、現在のテキストベースではちょっと難しい気がします。何か方法はあるでしょうか？。
　それから言語ですが、本当はFortranプログラマーなんですが（←絶滅危惧種(^^;)）、最近Visual Basicばっかり使ってまして、VBでもOKでしょうか？。

by ddtddtddt (2017-08-22 10:04) 

ddt³さん、こんにちは。

☆＞・・・BEMを用いたラプラス方程式（ポアソン方程式）の解法の簡単なサンプルプログラム・・・

となりますと、図や数式をちゃんと作らないと、と思うのですが、現在のテキストベースではちょっと難しい気がします。何か方法はあるでしょうか？。

◇わたしを信用していただけるのならば、この記事のコメントの設定を非公開にし、設定終了後に投稿されたコメントをブログで非公開、わたし以外のものは誰も見れないようにします。
それで、
わたしを信用していただけるのならば、
ddt³さんのメールアドレス――捨てメールアドレスで結構！！――をコメントに記し、そのコメントを送信してください。
そのコメントをいただければ、わたくし、ネムネコがそのメールアドレスへGmailでお返事のお便りをします。
わたしのGmailアドレス宛にワープロの文章、pdfファイルなどを添付したメールをいただければ、私がそれをブログ用の記事に編集し、ブログへ公開します。

この方法で良ければ、確認のコメントをください。

このとき、念の為に、メールアドレスは記さないでください。
送信したコメントが本当に非公開になっているか、そのコメントを見ることができない状態になっていることを確認してから、よろしければ、次のコメントにメールアドレスを記したコメントを送信してください。

メールアドレスは、使い捨てのもので結構です。
私がお便りするGmaiが届いた後、そのメールアドレスを削除しようが、変更しようが、構いません。


☆　それから言語ですが、本当はFortranプログラマーなんですが（←絶滅危惧種(^^;)）、最近Visual Basicばっかり使ってまして、VBでもOKでしょうか？。
◇私もForan（Fortran77）から入っています。
私はFortranで構いませんが、個人的にFortranのコンパイラを持っているヒトは少ないですからね〜。
Windowsユーザーですと、Cygwinなどをパソコンにインストールしなくてはならない。
ですが、
ddt³さんに、使用言語はお任せします。

ただ、私としたら、
できたら、誰でも簡単に使える十進BASICあたりでプログラムを書いて欲しいな〜
と思っています。

十進BASICのホームページ
https://goo.gl/jYDIFB

十進BASICは、確か１MBくらいの容量ですから、これをダウンロードしたとしてもPCを圧迫することはありません。
また、エミュレータですからパソコンの設定なども必要ありません。レジストリーをいじりませんから、こんなものは使えないと思ったら、ダウンロードしたものを削除するだけで事足ります。

ですが、
覚えることはほとんどありませんし、しかも、強力（？）で簡単なグラフィック機能も備えていますから、数値計算入門にはピッタリの言語ではないですかね。

ここを見ると、かなり本格的なシミュレーションも十進BASICでしているようですよ。
https://goo.gl/FTYYe1


by nemurineko (2017-08-22 15:01) 

ddt³さんへ。

記事の修正はこれ↓でよろしいですか？
https://goo.gl/5DBJzN


by nemurineko (2017-08-22 15:18)