プログラミングの記事を書きます。-korokoroこんにちは(@t_kun_kamakiri)(‘◇’)ゞたまに数値計算をして・・・時間刻みや空間刻みをどうやって決めれば安定して計算を進めてくれるかなって考えたりします。(時間刻みはこれくらいだろ!ってやってしまうこともたまにあります・・・ダメなやつだ(*_*;)そこで、本記事は以下の内容となります。本記事の内容はこちら数値計算の数値的安定性を考えてみよう(^^)/※1次元の熱伝導方程式を対象にします。カマキリ真面目に考えると面白いよ↓こちらの記事を題材にしたいと思います。前回の記事はこちらFortranからgnuplotを使って1次元グラフのアニメーション作成アニメーションはこんな感じhttps://takun-physics.net/wp-content/uploads/2019/04/20190429_185046.mp4カマキリ前回の記事を読んで遊んでみてくださいFortranからgnuplotを使って1次元グラフのアニメーション作成では、何をやっていたかを簡単におさらいしておきます。x}\tag{1}\end{align*}を次のようにして離散化して解きました。t}\end{align*}微分を差分に置き換えただけです。数値計算では時間をステップ数で数えるため、時間変数の代わりにステップ\(n\)として文字の頭に添えています。x^2}\end{align*}これは空間に対して2次の精度となります。よって、(1)式を離散化してものは次のようになります。x^2}\tag{2}\end{align*}このような方程式に限らず、解析解が見つかっていないような偏微分方程式を数値的に解く場合があると思います。精度:いかに元の方程式を忠実に解いているか?安定性:いかに安定して方程式を解いているか?以上の2点はとても重要です。x^2}\tag{2}\end{align*}のように、右辺を\(n\)時刻の変数の形にして表す離散化解法をオイラー陽解法といいます。時間微分に対して、1次精度で単純な差分の形に近似しているため、(2)式では直接\(n+1\)時刻の温度場を求めることができるのです。シンプルな数値的解法ではありますが、(2)式を式変形すると数値的に不安定な面が見えてきます。1次元熱伝導方程式のオイラー陽解法の安定性について※以降、\(T^{n}(x_{i})=T^{n}_{i}\)と書くこととします。x^2}\tag{3}\end{align*}これを、式変形します。x^2}\)と置きました。(4)式を眺めると、\(1-2c<0\)だと温度分布によっては、ある時刻で温度が負の値になるようなこともあるかなーーーってことに気付きます。0\)という条件が必要なのではないかなという気がするのです。\frac{1}{2}\tag{5}\end{align*}こうかな?というところまでわかりました。不安定な状況とは、はじめはおとなしい波を立てていたのに、時間が経つと振幅が大きい荒々しい波になってしまうような状況のことです。関数を波に分解したときに波が時間が経つにつれて増大してしまう状況が、数値的に不安定な状況であると言えます。それが偏微分方程式を離散化したときに起こってしまっては困るのです(‘_’)逆に、波の振幅が大きくてなってしまうような要因を数値解法の中に入っていないことが、安定性のポイントになるのです。カマキリ・・・ということをつい最近学びました。数値的に不安定な状況というのが何となく理解したところで、以下に数値的な安定性を調べる方法のポイントを整理しておきます。フーリエ級数展開)時間が経つと波の振幅が大きくなる(大きくなると不安定)このようなポイントで数値解の安定性を解析する方法を、ノイマンの安定性解析といいます。フーリエ変換は下記のようにある関数\(f(x)\)を、\(\sin\)波や\(\cos\)波や正弦波\(e^{ikx}\)などの様々な異なる波数の波に分解して波の空間(波数空間)に変換する方法のことです。どんな関数も異なる波数の\(sin\)波の和で書いてやるのです。※以下の場合は周期関数\(f(x)\)に関して、例を挙げています。x\)(周期1/10)の波が10個か・・・」と、波のくせに数えれるという恩恵を受けることができます。フーリエ級数展開)した後の関数(波数空間での関数)のこと「1個、10個、3個」は各波数の波の振幅と考えることができます。カマキリ↑上の説明はフーリエ変換というよりは、フーリエ級数展開したときの解説でありますが、似たようなものです(離散的か連続、周期関数か非周期関数かの違い)数値的に不安定な解とはこういうことです。「1個、10個、3個」などと書いた各波数の波の振幅が時間が経つにつれて、どんどん大きくなってしまう場合が不安定な数値解であるということです。A(k,t)e^{Ikx}\tag{5}\end{align*}これを、時刻\(n\)ステップと格子点\(i\)で書き換えます。x)}\tag{6}\end{align*}※虚数は\(I\)、格子点は\(i\)と書いています。紛らわしいですが・・・ちなみに\(\sum\)は波数\(k\)についての和です。x\)は離散化した際の格子幅です。数値的に安定とは、ある時刻\(n\)でのある格子点\(i\)の波数を考えた時に\(n+1\)時刻で増幅してしまう様な答えを導かないということです。(6)式を(4)式に代入して振幅\(A^{n}_{k}\)が大きくならないような条件を導こう(6)式を(4)式に代入するわけですが、元々の1次元の熱伝導方程式は線形の方程式であるので、各波数ごとに方程式が成り立ります。ある波数\(k\)についてのみ考えれば十分であるので、\(\sum\)はとっぱらって考えることにします。(6)式を(4)式に代入してみます。x)-1\big)\tag{8}\end{align*}となります。1\tag{9}\end{align*}が安定性の条件です。1\)という範囲に制限されているので、\(c\)にも条件が出てきそうであることがわかります。安定な数値解を得るためにはどのような波数\(k\)であっても、不等式(9)を満たしておく必要があるということです。x)=1\)の時です。x)=1\)のときにも不等式(9)に入っているための\(c\)である必要がありますが、この場合には\(c\)が消えてしまいます。ゆえに、どんな場合でも(9)式の上限側は成り立つことがわかります。x)=-1\)の時です。x)=-1\)のときでも、不等式(9)に入っているための\(c\)である必要があります。この不等式の「下限」の条件から、\(c\)の条件が決まります。\frac{1}{2}\tag{10}\end{align*}これが各波数の振幅が時間とともに増幅していってはいけないという条件です。不等式(10)が数値的に安定な条件ということになります。カマキリ(5)式と同じ結果が得られたまとめ今回は、下記の記事で数値計算した1次元の熱伝導方程式を題材にして数値的に安定な条件というのを導きました。空間微分に対しては2次精度で離散化すると、\begin{align*}T^{n+1}_{i}=cT^{n}_{i+1}+(1-2c)T^{n}_{i}+cT^{n}_{i-1}\tag{4}\end{align*}となります。x^2}\)ということになります。試しに↓下記の記事で遊んでみてください。離散化の際に「陰解法」など別の離散化手法にしたときに、数値的に安定な条件がどのように変わるかというのは、数値計算をする上では勉強しておいた方が良さそうですね。以下に参考文献を挙げておきます。【Fortranコードあり】Pythonで2次元の熱伝導方程式をアニメーションが簡単にできた【1階微分!中心差分の4次精度】テーラー展開から真面目に導出してみた。-【CAE解析の目的とは?】CAEは実機の代替えなのか?精度を無視してパラメータスタディすることは意味がないのか?【Udemyがお勧め理由5つ】プログラミングを勉強するなら動画学習おすすめ【まとめ】比熱の温度依存性について。古典論と量子論から理解する。【理想気体の比熱の温度依存性(7)】シュレディンガー方程式の極座標表記(2原子分子)から動径方向だけを抽出しようコメントを残すコメントをキャンセルプロフィールカマキリ(^^)大学以上の物理を勉強してブログに書いていきます。大学の専攻は物性理論で、Fortranを使って数値計算をしていました。CAE解析(流体解析)・毎日note更新しています。内容は、「現在勉強中の内容」「日々思ったこと」を書いています。1DCAEの操作方法などを動画にしています。youtube・Qiitaではプログラミング言語の基本的な内容をまとめています。(10)ブログをメールで購読メールアドレスを記入して購読すれば、更新をメールで受信できます。