対流項を持つ偏微分方程式の数値解法２ 風上差分の高精度化をはかる

対流項を持つ偏微分方程式の数値解法２　風上差分の高精度化をはかる

 

fは何度でも微分できる関数とすると、

（１）、（２）式からf''(x)を消去し、f'(x)について解くと、

となる。

ここで、

とおくと、

となり、２次精度の差分による近似式を得ることができる。

これは、の３点を用いて、曲線y=f(x)を２次関数で近似し、この近似式から微分係数を求めていることと同じこと。

だから、一般的に、の２点を用いて、

と後退差分を用いて求めた微分係数の近似値よりは精度がよい。

実際、f(x)=x²とすると、

となり、f(x)が１次関数、２次関数のとき、（３）式は正確に微分係数を求めることを確かめることができる。

 

前回、

の空間微分を

と風上差分（上流差分）で近似し、

という漸化式を得た。

ここで、

で、Cはクーラン数である。

 

（８）式は１次精度しか持っていないので、新たに得た２次の精度をもつ（左側）差分を用いて、

高精度化をはかることにする。

すると、（７）式は

したがって、

という漸化式を得ることができ、これを用いて逐次的に偏微分方程式（７）の数値解を求めることができる。

 

この方針にしたがって、c=0.5、Δt=Δx=1、つまり、クーラン数

について解いたものが次のグラフ。

比較参照のために、厳密解と１次精度の上流差分を用いて計算した結果は次の通り。

厳密解である垂直衝撃の前方では２次精度の風上差分の計算結果の方が厳密解に近いのだけれど、下流のx=4以降で負の値を持ち、その後、数値解が振動するなど、２次精度の上流差分は物理的にありえない数値解を求めてしまう。

 

また、c=0.5、Δx=Δt=1とし、クーラン数C=0.5とすると、時間の経過と共に、垂直衝撃波の下流での振動が激しくなる。１次精度の風上差分は

だから、風上差分の高精度にすると、安定に計算できるクーラン数の範囲が狭くなるようだ。

 

コメント 3

　あんまり関係ないかも知れないですけど、自分の業界では弾性体の動的挙動（波動伝播など）を解析するために、運動方程式に対してニューマークのβ法ってのを良く使います・・・っていうか、この方法が偏愛されまくってます（シンプレクテック法やルンゲクッタも使えよ！(^^;)）。

　ニューマークのβ法は2次精度の陰解法で、パラメータβの設定によって加速度の扱いが異なります。

　β＝1/6は線形加速度法と同じで、前進差分にほぼ対応します。つまりテーラー展開そのままです。
　β＝1/4は平均加速度法と同じで、中央差分にほぼ対応し、このとき無条件安定となってるのですが、クーラン条件の制約からはのがれられません。

　無減衰の場合は、β＝1/6でcΔt/Δx≦1/10程度にしないと経験的に危険と言われてます（たぶん理論もある）。
　β＝1/4でもcΔt/Δx≦1/6です。ちなみにシンプレクテック法，ルンゲクッタでも1/6以下にしないと、けっこう発散します。

　減衰付きの場合は1/3程度でけっこう行けます。現実の構造物に減衰がない訳ないので、うちの業界はかろうじて市販PCでやってける事になります(^^)。

　という訳で、クーラン数をもうちょっと小さくしたらどうかな？、って思いました。はずれの可能性大ですが・・・(^^;)。
by ddtddtddt (2017-11-30 12:47) 

こんばんは。そして、コメントありがとうございます。

空間微分に２次以上の差分を用いると、１階の波動方程式の数値解は不連続面（衝撃波）の前後で数値解が振動するらしいですね。
明日、数学の記事として上げる予定の２次精度のLax-Wendorffスキームも不連続面の上流側で振動します。
この振動を抑える方法は幾つは存在するようですが、これは空間微分に高次の差分をもちいる方法の宿命的なもので避けられないらしいです。

今回用いたシンプルな上流差分による計算は散々な結果でしたが、おそらく、もっと安定した２次の上流差分があると思います。

今回用いた上流差分で離散化した漸化式の係数を見ると、
1−3C/2
というものがあり、この係数が負になるとマズいので――このとき、数値的に安定でも解が負になる場合がある。優対角条件もどき(^^ゞ――
1−3C/2>0
ノイマンの安定性条件とは別に、数値解が意味を持つために
0<C<2/3
という制約があらたにつくのでしょうね。
数値的な安定性とは別に、対流（移流）項をもつ空間微分には条件がつくようなので、大変なようですね。

２階の波動方程式
∂²u/∂t²=c²∂²u/∂x²
は、制限は厳しいながら、安定的に精度よく解く陽解法が存在するのに、２階の波動方程式の一つの解にすぎない
∂u/∂t+c∂u/∂x=0
は、数値的にうまく解く方法がなかなかないのだから面白いですね。
もっともクーラン数が１のとき、この問題は、１次の上流差分で誤差なく正確に解くことができますが(^^)

それにしても、放物型、楕円型の偏微分方程式の数値解法とは異なり、波動方程式などの双曲型の偏微分方程式の数値計算は厄介ですね。

by nemurineko (2017-11-30 18:08) 

先のコメントの付け足しです。
今日取り上げたLax法は、クーラン数Cが0<C≦1のとき（数値開放的には）安定ですが、uの空間分布を見ると、ジグザクしていて、振動の様相を呈しています。
Lax法は、空間微分に２次の中心差分を用いていますので、このために数値解が振動の様相を示しているのだと思います。

by nemurineko (2017-11-30 18:28)