プログラミングの記事を書きます。-korokoroこんにちは(@t_kun_kamakiri)(‘◇’)ゞFortranを学生時代に使っていたのですが、最近になって使い始めています。完全に趣味です(^^)/それはそうと、Fortranは数値計算する上では直感的でとてもわかりやすい言語だと個人的には思っているのですが、めんどくさい点があるのです。それは、Fortranってグラフ作成などの可視化ツールがない・・・・(*_*;他のプログラミング言語に詳しくないから何とも言えないですが、Pythonとか対話型のプログラミング言語なら可視化のライブラリが用意されているのに・・・Fortranでよく使うグラフ作成ツールと言えば、gnuplot(ニュープロット)かなと思っています。カマキリgnuplotはタダだし・・・そこで次の疑問が・・・1次元グラフをどうやってgnuplotでアニメーションするのか?正直、↑これは学生の時からやり方わからなくて放置していたのですが、社会人になって改めて考えるとあっさりできたので記事にしておこうと思います。※Fortranのコードとスクリプトは本記事にアップされています。※FortranのコードはFortran90/95の書き方にしているつもりですが正確に書けているかは不明です。※FortranはWSLを使ってインストールしています。※gnuplotはこのサイトを参考にすれば良いかと思います。Linux(WSL)を使って、UbuntuからFortranやgnuplotをインストールして使用しています。こちらの記事もどうぞ1次元熱伝導方程式をFortranで解くまずは解くべき方程式を明らかにしておきましょう。x}\tag{1}\end{align*}熱伝導方程式は、初期状態を与えれば(1)式に従って時間発展をする方程式です。今回は初期状態を以下のようにしておきます。中央だけ中央だけ温度が高く、両端は温度0(K)の条件となっています。この初期状態から時間発展させると高い温度は低い温度へ熱移動が生じるので、下記のような変化を示します。これをgnuplotを使ってアニメーションを作ろうということになります。全体の流れアニメーション作成までの全体の流れを押さえておくと少し気持ちのハードルは下がるかもしれないので、先に全体の流れを示しておきます。1次元熱伝達方程式は偏微分方程式であるため、それを数値的に解くためには離散化してあげる必要があります。今回は、「時間に対しては陽解法」「空間に対しては2次精度」で解くことを考えています。次にForranを使って離散化した方程式を数値的に解くことを考えます。(※コードは後ほど)Fortranで出力したデータをgnuplotで順番に読み込んでいくためのスクリプトをFortranを使って作成します。これらのことができれば、あとはgnuplotで「3.」で作成されたスクリプトを実行するだけアニメーションを作成することができます。※1次元熱伝導方程式を題材にしたのは、簡単な方程式の方がアニメーション作成の方法を理解しやすいと思ったからです。ゆえに、離散化も特にこだわってはいません。時間微分の項(左辺)空間微分の項(右辺)をどのように扱うのかが問題です。いっちばん手っ取り早い離散化で進めていきましょう。t}\tag{2}\end{align*}微分を差分に置き換えただけです。数値計算では時間をステップ数で数えるため、時間変数の代わりにステップ\(n\)として文字の頭に添えています。x^2}\tag{3}\end{align*}これは空間に対して2次の精度となります。2次の精度とは関数をテーラー展開してときに何字の項まで残したかというものです。【補足説明です】テーラー展開については下記の記事をご参考ください。x^2}\)の2次精度(3)式を得ることができます。(2)(3)式より、数値計算で解くべき離散化の式がわかりました(↓以下です。x^2}\tag{4}\end{align*}では、これをFortranで解きましょう(^^)/Fortranで解くFortranで解くと言ってもプログラムを書かないことには何も始まらないので、プログラムを書きます。以下がプログラムです。thermal_explicitプログラムだけ載せられてもよくわからんという人もいるかもしれないので、説明を加えておきます(‘ω’)↑変数の定義をモジュールとしてまとめておきます。↑ここが今回の最も重要な部分で、アニメーションを作るためには各ステップ(出力したいステップ)のデータファイルを順に出力する必要があります。そのためのサブルーチンを書いています。↑ここからがメインのプログラム部分です。ここではメインプログラムで使うための変数を指定しています。↑ここで初期状態を作成しています。↑ここで(4)式の離散化した式を解いています。100ステップごとに温度の空間分布を出力するようにif文で制御しています。↑ちなみに「data」という名前のディレクトリを自分で作っておく必要があります。以上がプログラムの中身です。ここまでできたら、以下のコマンドで計算を実行させます。1D_Thermal_Eq.f90./a.out実行ができれば「data」という名前のディレクトリの中に、以下のように大量のデータが保存されていることが確認できます。これらのデータが指定したタイムステップで出力された温度の空間分布のデータです。あとは、これらのデータを順番にgnuplotで読み込めばアニメーションが完成ということになります。そのためには、gnuplotのスクリプトが必要になります。スクリプトは書くのが少々面倒なので、スクリプト自体をFortranで書くことにします。アニメーション用のgnuplotのスクリプトをFortranで作成gnuplotのスクリプトを作成するためのFortranのプログラムがこちらです。animation_scriptこれはあまり説明を必要としないのですが、簡単に説明すると「do文」のところで先ほどの「data」という名前のディレクトリ内の大量のデータを順番に読み込む作業を行っています。以下のコマンドを実行してみます。>movie」とするとmovieというファイルが出来上がり、その中にFortranの「print文」の内容を全て書き出してくれます。この「movie」というファイルがgnuplotのスクリプトとなります。以上のことが終えれば、あとはgnuplotで「movie」という名前のスクリプトを走らせるだけです。アニメーションを作成してみよう(^^)/「movie」というスクリプトは2通りの方法で実行できます。“movie”」とスクリプトファイルをloadしても良いです。アニメーションはこんな感じhttps://takun-physics.net/wp-content/uploads/2019/04/20190429_185046.mp4カマキリやったーまとめはじめての人にもできるように、わかりやすくまとめてみました(‘◇’)ゞが、説明がくどく長い記事になってしまったことをお許しください(‘ω’)まだまだFortranもgnuplotも勉強不足ですが、何とかできたかなっていう感じです。説明してしまいましたが、下記のzipファイルを解凍して下記のコマンドを打ってお試しください。-【1階微分!中心差分の4次精度】テーラー展開から真面目に導出してみた。-【CAE解析の目的とは?】CAEは実機の代替えなのか?精度を無視してパラメータスタディすることは意味がないのか?【Udemyがお勧め理由5つ】プログラミングを勉強するなら動画学習おすすめ【まとめ】比熱の温度依存性について。古典論と量子論から理解する。【理想気体の比熱の温度依存性(7)】シュレディンガー方程式の極座標表記(2原子分子)から動径方向だけを抽出しようコメントを残すコメントをキャンセルプロフィールカマキリ(^^)大学以上の物理を勉強してブログに書いていきます。大学の専攻は物性理論で、Fortranを使って数値計算をしていました。CAE解析(流体解析)・毎日note更新しています。内容は、「現在勉強中の内容」「日々思ったこと」を書いています。1DCAEの操作方法などを動画にしています。youtube・Qiitaではプログラミング言語の基本的な内容をまとめています。(10)ブログをメールで購読メールアドレスを記入して購読すれば、更新をメールで受信できます。