［放物型偏微分方程式にシンプレクティック法は使えないか？（使えない！）］

［放物型偏微分方程式にシンプレクティック法は使えないか？（使えない！）］

　ddt^3です。具体的にシンプレクティック法プログラムで計算してみました。例題はネコ先生と同じです。そして結論は、使えない！となりました。以下はドツボにハマり悪戦苦闘したあげく、大敗した経緯です(^^;)。

　まずシンプレクティック法を強引に使うために導いた時間に関する2階微分方程式、

　　　　　　　　　　　　　　　　　　　

 

 

 

は、熱伝導方程式、

　　　　　　　　　　　　　

 

 

 

を時間tで微分して作ったものなので、(1)は(2)に対する必要条件です。不用意な条件でやるとあらぬ解が出そうですが、(1)の時間積分手順を考えてみると、原理的には初期条件さえ決めればシンプレクティック法で以後の全ての挙動は決まるので、(1)の解が(2)の解と一致するためには、初期条件の設定に敏感そうだと想像できます。

　初期条件は0≦x≦4でT(x，0)＝－x²＋4x，境界条件はT(0，t)＝T(0，4)＝0です。またλ＝0。式(1)は、重み付き残差法の弱形式を通じて有限要素法に持ち込み、節点未知数として図-1の節点温度Tと温度勾配wを取ったので（(3)は形状関数，ξ＝x/4）初期条件として、初期温度と温度勾配分布および初期温度速度と勾配速度分布を与える必要があります。それには熱伝導方程式(2)を用いるのが妥当です。一要素の節点は両側の隣接要素と共有されるので、式(3)はC¹級の3次近似になります。

 

 

 

 

 

 

 

 

 

 

　初期温度は明らかにT(x，0)＝－x²＋4xです。初期温度勾配は次のように決めました。解析領域の両端x＝0，4では温度一定なので(2)より、T(x，t)のxでの2階微分は0になる必要があります。この条件が成立するようにx＝0，4の傾きw₁(0)とw₂(4)を決めます。中間節点位置x＝x_jでは要素間でxでの2階微分が連続になるようにw(x_j)を定めます。

　解析領域をn分割するとすると（n要素）、節点はn+1個あり、一節点2自由度なので未知数は2(n+1)個。そのうち初期温度分布からn+1個の温度値が決まり、両端での2階微分の値と中間での2階微分の連続性から、勾配値n+1個に対する条件も得られので、初期温度分布さえ与えれば温度分布と温度勾配分布を計算できるようになります。ただしこれは、2(n+1)個の節点未知数に関する連立一次方程式を解く事になります。結果はC²級の3次近似です（← すでにドツボの予感(^^;)）。

　温度速度分布については(2)から、xでの2階微分がわかれば良い訳ですが、さっき2階微分が連続になるように温度勾配を与えたので、(3)をξで（xで）2階微分すれば温度速度分布が得られます。

　勾配速度分布については、(2)をxで微分したものから得られるので、xでの3階微分です。これも(3)をξで3階微分すれば得られますが、(3)は3次関数なので一般には不連続です。でも(2)には、xでの3階微分は直接的には出てこないので、こんなところでどうだろう？と始めました。

　図-1の横軸は時間，縦軸は解析領域中央（x＝2）における温度変化です。d-ｵｲﾗｰ1次のラインは、通常の差分を行い1次のオイラー法で時間積分した結果で、ここでは厳密解とみなします。

　初期条件をC²級近似しシンプレクティック法に持ち込んだ結果を、s-離散-Non.Corで表します。御覧のようにいともあっさり発散します。t＝10では2.86×10²⁰⁴に達しました(^^;)。離散はLamped Massの意味。Consistent Massでもやってみましたが結果は変わりません。

　ここで気づいた事があります。式(1)でλ＝0とすると、

　　　　　　　　　　　　

 

 

 

ですが、これは正のフィードバック系です。

　自分が良く扱う梁の曲げ振動は、

　　　　　　　　　　

 

 

 

という形をしてますが、これは負のフィードバックです。右辺にマイナス符号があるので、Tが大きくなり過ぎると逆向きの運動をさせようという加速度（復元力）が働らくので、安定した振動解になります。しかし(2)は正のフィードバック系なので、Tがちょっとでも大きくなり過ぎるとどこまでも加速して行きます（どんだけぇ～(^^)）。つまり数値誤差にも非常に弱い、という事になります。

 

　という訳で修正項（安定化機構）が必要でないか？と考えました（←ドツボへの入り口）。そこで重み付き残差法のもともとの形に戻ると、

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　

 

でした。(4)右辺の2，3項目は物理量の連続性より0としました。でも今は初期条件を除きC¹級近似なので、数値上は2，3項目も0ではありません。(3)などでは物理量の連続性から、これらを0とする強制条件を与えた方が上手く行くのですが、(2)が非常に数値儀差に対して敏感な系だとわかったので、これらの値も考慮すれば安定するのでは？と考えました。(4)右辺の2，3項目は、要素端でのいわば不釣り合い力を表すからです。それらを考慮するという事は、数学的には偏微分方程式をかなり人為的に操作する事にもなるのですが、系が超発散しやすいので、仕方ないだろうという訳です。

 

　最初に言いますがみなさん、こんな事は絶対に考えてはいけませんよ。数値解法が誤差に非常に弱いとわかった瞬間にあきらめるべきです。でもやってみたかったんだよなぁ～(^^;)。という訳で時間ステップの各段で(4)の境界項を考慮した結果が、図-1の赤ライン：s-離散-Cor.です。

　結果はある意味で安定しました。温度はt＝0から時間に比例して下がり続け、t＝10では－6.038。2.86×10²⁰⁴と比較すれば大進歩です。・・・やっぱり不釣り合い力が影響してるのだ。という事は不釣り合い力が出ないようにしてやれば・・・と、ますますドツボへ。

　初期条件では温度分布からC²級近似してるから、(4)右辺の2階微分項は0になります。だったら時間ステップの各段で得られた温度分布から、同じように温度勾配，温度速度，勾配速度をC²級近似になるように決め直したら？。結果は水色点線のs-離散-Cor.-C2です。オ～～、s-離散-Cor.と全く同じ結果だ。やはり不釣り合い力だ！。

　(4)右辺の3階微分項も0になるようにC³級近似しよう。それには最低4次近似必要だが、(3)の形状関数に4階微分値の定数項を未知数として含めるだけだから、なんとかなるだろう。温度分布は事前に与えられるとして、温度勾配は最初と同じように2階微分の両端での値と要素間での連続性から決め、3階微分も要素間での連続性から決める。・・・とやってみると、なんと条件が1個足りなくなりました。3階微分の値の1個が不定になるのです。普通ならこの辺りでやめるのですが・・・(^^;)。

 

　ところがドツボ状態だと、都合の良い方へ都合の良い方へと確証バイアスが働きます。3階微分の値の1個が不定にならならそれは自由に選べるから、最適近似を選べるという事じゃないか。これは良い事だ！。最初に従来の方法でC²級近似を行い、得られた温度勾配分布との差が最も小さくなるように、最小2乗法で3階微分の値を選んでやればどうだろう？。やってみると、なかなか綺麗なC³級近似が得られます。その計算結果が、s-離散-C3です。t＝10で5.23×10⁹⁶。最初に戻っちゃった・・・(^^)。

 

　C³級近似すると、(4)の境界項はいつも0です。予想ですが、少なくとも空間方向の近似が厳密解と同等程度でなければ、この方法は上手く行かないのでしょう。1次元ならば、それは空間方向の境界要素法を併用する事でほぼ可能ですが、仮に上手く行ったとしても、非常に不安定な系を解いてる事に変わりはありません。ちょっと激しい変動をする熱源密度なんかを付けただけで発散しない保証は0です。また2次元以上では境界要素法でも、空間近似の影響が入ります。ここで万策尽きました。

　これだけ手をかけてもまともな応えが得られないという事は、放物型へのシンプレクティック法の使用は、明らかに失敗です。理論と数値計算は違うのだ、と言い訳しておきます。

 

　撤収ゅ～～！(^^;)。

 

コメント 2

　ネコ先生やばいっすよ。ついに順位が9位っすよ(^^)。
　　　　　　　　　　　　　　　　　　↑
　　　　　　　　　　　　　　（ニコニコするな）
by ddtddtddt (2018-01-15 18:41) 

ご心配、ありがとうございます。
また記事をコツコツと書き、順位があがるように努力します。

それはそれとしまして、インフルエンザにかかり、現在、大変なことになっております。
　――悪寒、筋肉痛、咳、頭痛などなど――
体質的に熱は出ないので、この点は助かっておりますが、ただでもない集中力が持続しないのがツライですね。だから、いま、数学の記事を書くのはちょっと無理ですね。

by nemurineko (2018-01-15 21:08)