ソリトン分裂に関して(6)

格子間隔5m、最大水深10m、最低水深1m、α=0.5、β=0.9で概ね安定して計算が出来た。つまり、この場合は時間ステップを半分にしてやれば修正係数は不要になる。が、格子間隔2m、最大水深90m、最低水深1mではやはり数値発散してしまった。
 
安定するまでα、βの値を調整するとα=0.005、β=0.01となってしまった。つまり、時間ステップをさらに0.005倍して計算すればいい筈だがこれは途轍もない計算時間を要する。時間ステップをそれ程小さくしたくない。(小さくしずぎると別の問題、数値分散誤差が無視できなくなる)すると、これは分散項を無視しているのと変わりない。
 
もう一度、原因について振り返ってみると主な原因は時間微分項の精度不良だと考えていたが良く考えてみるとある程度近似されていると考えられるからそれ程悪くは無いと思う。そう考えると問題の式は
イメージ 1
なのだからh/Δx のような因子がやはり怪しい。これはこれで何となく分かるがこれは計算間違いのしようが無い因子だからF.G,Hに主な原因を求めていた。
ここで思考停止。もう分けが分からない。
 
これでは先に進めないので問題は保留にしよう。仕方が無いのでh/Δxに問題があると仮定して考えると、つまり分散項の評価はh/Δxが十分小さな水域でしか適用出来ないという制約があるという事になる。
では、h/Δxがどれ位だと計算が安定するのだろうか?色々と条件を変えてみると水深と格子間隔に関して分散項の計算安定性は極めて大きな関係がある事が分かった。
ん、、、これはやはり計算の仕方が悪いのか、、、。
色んな文献を調べてみてもそんな関連は見当たらなかった。
そもそもそんな関係は変な気がするし、それでは格子間隔が小さい精度の高いデータで水深が制限される事になってしまう。
結論から言うと Δx  <  ζh 1.5 <ζ<2 という関係があった。
 
少しキーワードを変えて調べてみると次の論文が見つかった。

非線形分散波理論の津波数値解析への適用性と新しい数値モデルの提案」
鴫原良典、今村文彦
 イメージ 2
 
私の解き方はどれにも該当しないので単純な比較は出来ないが、計算の安定性に
Δx  <  ζh という関係が示されている。
こうなると非線形分散波方程式(修正ブシネスク方程式(Madsen and Sorensen1992))が深海から極浅海域までカバーするという意味が良く分からない。
 
確かに領域を分けて浅い海域では小さい格子、深海域では大きめの格子で対応できるし、実際の津波解析では普通に行われているので間違ってはいない。しかしこれは例えば分散項を含まなければそういうわけ方でも問題は無い。しかし、分散項を含んでいる状態で上記の制限があるならこのわけ方は基本的には出来ない話だ。
分散項を含まなければそういうわけ方が出来るというのは単純な話だ。
基本的には沿岸部に近いところを精密に評価したいわけだからそれで良い。そして沿岸部に近いところは基本的には、良いですか基本的には、水深は浅いし、外洋に行けば深くなるという理屈に矛盾しない。沿岸部に深い場所や外洋に極浅の箇所があっても問題ない。しかし、分散項を含んだ形ではその理屈は相応しくない。つまり、沿岸部に深い場所存在した時点で計算が出来なくなる。
分散項を含んだ形では格子間隔と水深に制約があって良いとするならそういう制約関数を導入して方程式を次のように出来る。
イメージ 3
r(Δx, h)が制約関数でΔx  <  ζh で1.0 Δx  ≧ζh なら水深に応じて小さくなるような関数。一応、修正係数α、βは残しておいた。さらに、移流項を数値分散を押させるために3次精度風上差分に変更した。
 
問題は制約関数をどの様に定義するかだ。試行錯誤した結果
イメージ 4
イメージ 5
ただし、γは1以下の値で0.90とした。また、以下の条件では1.0となる関数である。
関数としては以下のような形状だ。(4m格子の場合)
イメージ 6
つまり、制約内では1.0でありそこから外れるほど小さくなるような関数である。これが旨く機能するなら以前導入した補正係数α、βは不要になる。補正係数は計算安定圏内でも一律に値を押さえ込んでしまうため分散項を一律にカットしてしまうが上記の制約関数を導入すれば安定圏内では人為的な補正は入らない事になるし安定圏外でも安定圏内に近いほど無駄にカットされないという利点がある。
 
テストしてみた。
殆どあらゆる海底地形で安定した計算が出来た。しかし、問題は非線形長波理論と分散項を含んだ非線形分散波理論でまずは有意な差がでるのか?そして問題はソリトン分裂を再現出来るかである。