Fortran

【1次元熱伝導方程式】Fortranからgnuplotでグラフの1次元アニメーション作成

こんにちは(@t_kun_kamakiri)(‘◇’)ゞ

Fortranを学生時代に使っていたのですが、最近になって使い始めています。
完全に趣味です(^^)/

それはそうと、Fortranは数値計算する上では直感的でとてもわかりやすい言語だと個人的には思っているのですが、めんどくさい点があるのです。

それは、Fortranってグラフ作成などの可視化ツールがない・・・・(_;

他のプログラミング言語に詳しくないから何とも言えないですが、Pythonとか対話型のプログラミング言語なら可視化のライブラリが用意されているのに・・・

Fortranでよく使うグラフ作成ツールと言えば、gnuplot(ニュープロット)かなと思っています。

カマキリ

gnuplotはタダだし・・・

そこで次の疑問が・・・
疑問

1次元グラフをどうやってgnuplotでアニメーションするのか?

正直、↑これは学生の時からやり方わからなくて放置していたのですが、社会人になって改めて考えるとあっさりできたので記事にしておこうと思います。
※Fortranのコードとスクリプトは本記事にアップされています。
※FortranのコードはFortran90/95の書き方にしているつもりですが正確に書けているかは不明です。
※FortranはWSLを使ってインストールしています。
※gnuplotはこのサイトを参考にすれば良いかと思います。
僕の使用環境はWindows10でWindows Subsystem for Linux(WSL)を使って、UbuntuからFortranやgnuplotをインストールして使用しています。
カマキリ

説明がくどいな~って思った方は「まとめ」へGo!!

スポンサーリンク

1次元熱伝導方程式をFortranで解く

まずは解くべき方程式を明らかにしておきましょう。

1次元熱伝導方程式
\begin{align*}\frac{\partial T}{\partial t}=\kappa\frac{\partial^2T}{\partial^2 x}\tag{1}\end{align*}
熱伝導方程式は、初期状態を与えれば(1)式に従って時間発展をする方程式です。
今回は初期状態を以下のようにしておきます。
中央だけ中央だけ温度が高く、両端は温度0(K)の条件となっています。
ですので、この初期状態から時間発展させると高い温度は低い温度へ熱移動が生じるので、下記のような変化を示します。
これをgnuplotを使ってアニメーションを作ろうということになります。

全体の流れ

アニメーション作成までの全体の流れを押さえておくと少し気持ちのハードルは下がるかもしれないので、先に全体の流れを示しておきます。

  1. 1次元熱伝達方程式は偏微分方程式であるため、それを数値的に解くためには離散化してあげる必要があります。
    今回は、「時間に対しては陽解法」「空間に対しては2次精度」で解くことを考えています。
  2. 次にForranを使って離散化した方程式を数値的に解くことを考えます。
    (※コードは後ほど)
  3. Fortranで出力したデータをgnuplotで順番に読み込んでいくためのスクリプトをFortranを使って作成します。

これらのことができれば、あとはgnuplotで「3.」で作成されたスクリプトを実行するだけアニメーションを作成することができます。

※1次元熱伝導方程式を題材にしたのは、簡単な方程式の方がアニメーション作成の方法を理解しやすいと思ったからです。
ゆえに、離散化も特にこだわってはいません。

1次元熱伝導方程式の離散化

1次元熱伝導方程式
\begin{align*}\frac{\partial T}{\partial t}=\kappa\frac{\partial^2T}{\partial^2 x}\tag{1}\end{align*}
この偏微分方程式を離散化するにあたって、
  • 時間微分の項(左辺)
  • 空間微分の項(右辺)

をどのように扱うのかが問題です。

いっちばん手っ取り早い離散化で進めていきましょう。

ある地点(x_{i})での時間微分の項

\begin{align*}\frac{\partial T}{\partial t}\simeq \frac{T^{n+1}(x_{i})-T^{n}(x_{i})}{\Delta t}\tag{2}\end{align*}
微分を差分に置き換えただけです。
数値計算では時間をステップ数で数えるため、時間変数の代わりにステップ\(n\)として文字の頭に添えています。
ある時間\(t\)での空間微分の項
\begin{align*}\kappa\frac{\partial^2T}{\partial^2 x}\simeq \kappa\frac{T^{n}(x_{i}+\Delta x)-2T^{n}(x_{i})-T^{n}(x_{i}-\Delta x)}{\Delta x^2}\tag{3}\end{align*}
これは空間に対して2次の精度となります。
2次の精度とは関数をテーラー展開してときに何字の項まで残したかというものです。

【補足説明です】
テーラー展開については下記の記事をご参考ください。
前回の記事はこちら
\begin{align*}T(x_{i}+\Delta x)=T(x_{i})+\frac{\partial T}{\partial x}|_{(x_{i})}\Delta x+\frac{1}{2}\frac{\partial^2T}{\partial x^2}|_{(x_{i})}\Delta x^2+o(\Delta ^3)\end{align*}
\begin{align*}T(x_{i}-\Delta x)=T(x_{i})-\frac{\partial T}{\partial x}|_{(x_{i})}\Delta x+\frac{1}{2}\frac{\partial^2 T}{\partial x^2}|_{(x_{i})}\Delta x^2+o(\Delta ^3)\end{align*}
この両方の差から、\(\frac{\partial^2 T}{\partial x^2}\)の2次精度(3)式を得ることができます。
補足説明終了です。

(2)(3)式より、数値計算で解くべき離散化の式がわかりました(↓以下です。)
\begin{align*}T^{n+1}(x_{i})=T^{n}(x_{i})+\kappa dt\frac{T^{n}(x_{i}+\Delta x)-2T^{n}(x_{i})-T^{n}(x_{i}-\Delta x)}{\Delta x^2}\tag{4}\end{align*}
では、これをFortranで解きましょう(^^)/

Fortranで解く

Fortranで解くと言ってもプログラムを書かないことには何も始まらないので、プログラムを書きます。

以下がプログラムです。

プログラムだけ載せられてもよくわからんという人もいるかもしれないので、説明を加えておきます(‘ω’)

↑変数の定義をモジュールとしてまとめておきます。
↑ここが今回の最も重要な部分で、アニメーションを作るためには各ステップ(出力したいステップ)のデータファイルを順に出力する必要があります。
そのためのサブルーチンを書いています。
↑ここからがメインのプログラム部分です。
ここではメインプログラムで使うための変数を指定しています。
↑ここで初期状態を作成しています。
↑ここで(4)式の離散化した式を解いています。
100ステップごとに温度の空間分布を出力するようにif文で制御しています。
↑ちなみに「data」という名前のディレクトリを自分で作っておく必要があります。
以上がプログラムの中身です。
ここまでできたら、以下のコマンドで計算を実行させます。

実行ができれば「data」という名前のディレクトリの中に、以下のように大量のデータが保存されていることが確認できます。

これらのデータが指定したタイムステップで出力された温度の空間分布のデータです。

あとは、これらのデータを順番にgnuplotで読み込めばアニメーションが完成ということになります。

そのためには、gnuplotのスクリプトが必要になります。

スクリプトは書くのが少々面倒なので、スクリプト自体をFortranで書くことにします。

アニメーション用のgnuplotのスクリプトをFortranで作成

gnuplotのスクリプトを作成するためのFortranのプログラムがこちらです。

これはあまり説明を必要としないのですが、簡単に説明すると「do文」のところで先ほどの「data」という名前のディレクトリ内の大量のデータを順番に読み込む作業を行っています。

では、以下のコマンドを実行してみます。

「./a.out >movie」とするとmovieというファイルが出来上がり、その中にFortranの「print文」の内容を全て書き出してくれます。

この「movie」というファイルがgnuplotのスクリプトとなります。

以上のことが終えれば、あとはgnuplotで「movie」という名前のスクリプトを走らせるだけです。

アニメーションを作成してみよう(^^)/

「movie」というスクリプトは2通りの方法で実行できます。

方法1

↑このようにコマンドを打っても良いですし、

方法2

↑このように「gnuplot」と打って、「load “movie”」とスクリプトファイルをloadしても良いです。

アニメーションはこんな感じ

カマキリ

やったー

まとめ

はじめての人にもできるように、わかりやすくまとめてみました(‘◇’)ゞ

が、説明がくどく長い記事になってしまったことをお許しください(‘ω’)

まだまだFortranもgnuplotも勉強不足ですが、何とかできたかなっていう感じです。

ごちゃごちゃ説明してしまいましたが、下記のzipファイルを解凍して下記のコマンドを打ってお試しください。

20190428_thermal1D

COMMENT

メールアドレスが公開されることはありません。 が付いている欄は必須項目です