2021年5月4日こんにちは(@t_kun_kamakiri)。今回は、Pythonを使って「1次元の移流方程式」のクーラン条件のお話をします。クーラン条件は数値的安定性とも関係があるので、必ず学んでおいた方が良いですね。本記事では、わざと数値的に不安定な状態での数値計算もやりますので、解の振る舞いがどのようになるのかを確認しながら見ていただければと思います。本件の基本的な内容はこちらのサイトにそってやっていきます。この記事ではこんな人を対象にしています。Colaboratoryでアニメーション作成。前回は「1次元の移流方程式」のアニメーション作成をしましたが、今回はそのアニメーションを使って数値的安定性について考えていきたいと思います。今回の内容はこちら移流方程式における安定性条件(クーラン条件)を考える。Colaboratoryを使ってコードを書きながら理解を深めていきたいと思います。1次元移流方程式を後退差分で解くプログラムをPythonで書きました。x}\big(u^{n}_{i}-u^{n}_{i-1}\big)\tag{1}\end{align*}です。数値的安定性(1)式を離散化して、時々刻々と変動する数値解を安定的求めるには、計算を進めるにあたって数値的な誤差が増大しないことが条件になります。移流方程式には厳密解が存在するので、厳密解をとして\(\tilde{u}^{n}_{i}\)、数値解を\({u}^{n}_{i}\)とします。すると、厳密解と数値解には誤差があるので誤差を\(\tilde{\epsilon}^{n}_{i}\)と置いておきしょう。「いかに数値誤差\(\tilde{\epsilon}^{n}_{i}を小さくすることができるか\)」です。(2)式を(1)式に代入します。x}\big(\tilde{u}^{n}_{i}-\tilde{u}^{n}_{i-1}\big)\end{align*}を満たしています。移流方程式が線形な方程式であるため、線形解である\(u^{n}_{i}=\tilde{u}^{n}_{i}+\tilde{\epsilon}^{n}_{i}\)のうちの右辺のどちらも移流方程式を満たすという事です。今、考えるべきは以下の式です。x}\big(\tilde{\epsilon}^{n}_{i}-\tilde{\epsilon}^{n}_{i-1}\big)\tag{4}\end{align*}(4)式は厳密解というよりは、そもそも誤差なので\(\tilde{\epsilon}^{n}_{i}=0\)の予定ですが、(4)式がある以上はこの式に従ってもらいながら発散されるととても困ることになります・・・・この誤差\(\tilde{\epsilon}^{n}_{i}\)が時間とともに増大するかどうかを考えて、増大するような数値解法だと解が不安定になるというわけです。カマキリどうやって増大するかどうかを調べるのか?例えば、以下のようにまとめると数値的安定性条件が出てきそうなのがわかります。0\end{align*}じゃないと、解が安定しないのではないかというのが予想できます。(4-1)式の左辺は「n+1ステップ」の「格子点i」での誤差に対して、右辺は「nステップ」の「格子点iとi-1」の誤差です。あるステップのときのある場所ではマイナスの値をとったりプラスの値を取ったりして解が不連続になりやすいと予想できます。すると、不連続点の位置での微分は値が発散するということになるので(今は離散化しているのでせいぜい大きな値をとるくらい)、解が不安定になりやすい・・・・とちょっぴり予想できます。今回使う方法は、「フォン・ノイマンのフーリエ解析法」と呼ばれるものです。名前だけで圧倒されそうな解法ですが、ポイントは以下となります。フーリエ級数展開)フーリエ変換は下記のようにある関数\(f(x)\)を、\(\sin\)波や\(\cos\)波や正弦波\(e^{ikx}\)などの様々な異なる波数の波に分解して波の空間(波数空間)に変換する方法のことです。どんな関数も異なる波数の\(\sin\)波の和で書いてやるのです。※以下の場合は周期関数\(f(x)\)に関して、例を挙げています。x\)(周期1/10)の波が10個か・・・」と、波のくせに数えれるという恩恵を受けることができます。フーリエ級数展開)した後の関数(波数空間での関数)のこと「1個、10個、3個」は各波数の波の振幅と考えることができます。結局複雑な波形でも、ある波数を持つ波の重み付けによる重ね合わせで表現できるのです。数値的に不安定な解とはこういうことです。「1個、10個、3個」などと書いた各波数の波の振幅が時間が経つにつれて、どんどん大きくなってしまう場合が不安定な数値解である。誤差\(\tilde{\epsilon}^{n}_{i}\)をフーリエ級数展開して、ある波数の波が増大する可能性があるのかどうかを調べることで、解の安定性条件が導くことができます。フーリエ変換について学びたい方は☟こちらからどうぞ!「フーリエ変換」の記事はこちら大学数学では、本題に行きます(^^)/時間が経つと波の振幅が大きくなる(大きくなると不安定)誤差\(\tilde{\epsilon}^{n}_{i}\)、時刻\(n\)ステップと格子点\(i\)で以下のようにフーリエ級数展開を行います。x)}\)が\(A^{n}_{k}\)で重み付けされた形になっています。※虚数は、格子点の\(i\)と紛らわしいので「虚数\(I\)」と表記しました。これを(4)式に代入します。(4)式は線形な方程式なので(5)式の和のほとんどは書かなくてもよくなります。x}\)と置いています。これを今、以下のようにまとめます。x+yI\end{align*}※実部と虚部に分けました。\begin{align*}|g|=\bigg(4\sin^2\frac{\phi}{2}\big(\alpha-\frac{1}{2}\big)-\sin^2\frac{\phi}{2}+1\bigg)\tag{7}\end{align*}となります。わかりにくいので絵でイメージすると以下のようになります。\frac{A^{n+1}_{k}}{A^{n}_{k}}\)が1より大きいという事は、ある波数\(k\)の重みがどんどん増大してやがて数値的には無視できないくらいの誤差になるということを意味しています。1\,\,\,\,(不安定)\end{matrix}\right.\end{align*}という判別ができます。1\)となる条件を導けばよいという事になります。1\tag{8}\end{align*}となります。1\)を満たしますね。先に予想していたような結果を得ることができました。x}\)をクーラン数といいます。クーラン数の値を変えてシミュレーションをしてみましょう。以下の3パターンで解析をします。?カマキリさて、どうなるでしょう?Pythonコードを書きます。html='jshtml')animここから\(\alpha\)の値を変えて、アニメーションを見てみます。=0.5\)この条件は、クーラン条件からすると安定的に解ける条件のはずですがどうなるでしょうか。☟動画を見てください。https://takun-physics.net/wp-content/uploads/2020/05/20200506_232031.mp4解が安定はしていますが、「+x方向に進みながら波形の形が変わっている」という結果でした。=1.0\)続いて・・・この条件は、クーラン条件からすると安定的に解ける条件のはずですがどうなるでしょうか。☟動画を見てください。=1.5\)これはどうでしょうか?この条件は、クーラン条件からすると不安定な解の振る舞いをする条件のはずですがどうなるでしょうか。☟動画を見てください。1\,\,\,\,(不安定)\end{matrix}\right.\end{align*}の条件と一致している数値計算の結果ですね。1\tag{8}\end{align*}となることを示しました。(8)式をクーラン条件といいます。ここで注意しておきたいのは、本記事の内容は「\(\alpha>0\)かつ後退差分をした際」の安定条件についてですが、\(\alpha<0\)である場合に関しては、後退差分ではいかなる条件も不安定な解の振る舞いをします。\(\alpha<0\)の条件で数値計算を行うのであれば前進差分で解く方が安定性が良いです。それは以下の記事でもまとめていますので、良ければ参考にしてください。コード自体はFortranを使用していますが、コードは無視して動画を見るだけでも勉強になると思います。上記の記事で以下の結果を得ています。次回の記事はこちら【第7回Python流体の数値計算】1次元拡散方程式を差分法で実装する。Pythonの完全初心者は書籍で学ぶとよいPython自体が不安だという方は、こちらの書籍から勉強しても良いかも知れません。☟こちらの本がPython初心者が挫折することなく勉強できる本です。04:45時点)Amazon楽天市場YahooAmazonの情報を掲載しています僕自身もこの書籍ではじめは勉強しました。小難しい余計なことが書いていないし、年末の一週間くらいで読むことができたので、初心者の方にはとてもお勧めです。流体の数値計算以前に、流体力学を勉強しないと!!って思っている方は以下の参考書からはじめてみるのが良いと思います。09:20時点)Amazon楽天市場YahooAmazonの情報を掲載しています工学よりに書かれていて実用的な内容も触れながら、流体の理論的な内容も易しく書かれていますので大学1年生で偏微分を学んだ方にとってはちょうどよい参考書です。数値計算は微分を離散化して代数的に解くため、どうしても誤差が発生します。それをいかに誤差が少なく、安定的に解くかがカギになっており、数値的解法は多く提案されています。以下の参考書のような、数値計算に特化した参考書を一冊手元に置きながら勉強を進めることで「流体の数値計算」の理解が深まると思います。kamakiri宇宙に入ったカマキリPython【第2回Python流体の数値計算】Numpyの使い方を確認する。し」と「値渡し」の違いを理解する。コメントをキャンセルメールアドレスが公開されることはありません。Colab...【第7回Python流体の数値計算】1次元拡散方程式を差分法で実装する。製造業コミュニティ製造業コミュニティは製造業で働く方や製造業へ就職を考えている学生の方が匿名で気軽に情報交換を出来るオンラインコミュニティです。本が聴けるアプリAudible(オーディブル)☟30日間の無料体験(最初の1冊は無料)プロフィールカマキリ(^^)大学の専攻は物性理論で、Fortranを使って数値計算をしていました。CAEを用いた流体解析は興味がありOpenFOAMを使って勉強しています。CAE解析(流体解析)notenoteで内容は主に「プログラミング言語」の勉強の進捗を書いています。「現在勉強中の内容」「日々思ったこと」も日記代わりに書き記しています。1DCAEの操作方法などを動画にしています。(音声はありません_(._.)_)製造業コミュニティをやっています。Facebookを使った実名登録。製造業ブロガーコミニティ案内オンライン家庭教師のご要望あればご連絡ください。Amazonのアソシエイトとして、当メディアは適格販売により収入を得ています。・当メディアはGoogleによるユーザーに合わせた自動広告が表示されます。Twitterカマキリ☟こちら日々の「思ったこと」「勉強している内容」を発信しています。4Pythonを使ってTwitterAPIで自作botを作りました。製造業コミュニティの方々と運営しています。botを作るのかを知りたい方はこちら。【超丁寧に解説】Pythonを使ってTwitterAPIでGoogleスプレッドシートに書いた内容を画像付きでツイート投稿する。宇宙に入ったカマキリ。