Pythonで1次元熱伝導方程式のグラフのアニメーションが簡単にできた

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

前回の記事で1次元熱伝導方程式をFortranで解いて、1次元のグラフのアニメーションをgnuplotで作成するというのをしました。

前回の記事はこちら

↑この記事でのアニメーション結果はこちら。

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

Fortranは可視化ツールがないので仕方なくgnuplotのスクリプトを書いてグラフのアニメーションを作りました。

カマキリ

めんどうだった・・・

本記事の内容

今回は、1次元熱伝導方程式をPythonで解いてグラフのアニメーションまで作成したいと思います。
Pythonの初学者のための記事です。

※Pythonの動作にはJupyter Notebookを使用しています。

Pythonの使用環境
1次元熱伝導方程式を離散化して解く
Pythonで解いて1次元のアニメーション
Pythonのアニメーションの結果
参考書

Pythonの使用環境

Pythonの使用環境は色々あると思います。

以下に代表的なものを紹介しておきます。

①Pythonの本家をインストール

②AnacondaをインストールしてJupyter notebookを使用

③PycharmをインストールしてPythonを使用

④visual studioをインストールしてPythonを使用

カマキリ

僕はAnacondaをインストールしてJupyter notebookを使用したいと思います。

1次元熱伝導方程式を離散化して解く

解きたい方程式は、下記の記事でも詳しく解説しています。

前回の記事はこちら

Fortranからgnuplotを使って1次元グラフのアニメーション作成

Pythonで解きたい方程式は↓こちらです。

1次元熱伝導方程式

\begin{align*}\frac{\partial T}{\partial t}=\kappa\frac{\partial^2T}{\partial^2 x}\tag{1}\end{align*}

熱伝導方程式は、初期状態を与えれば(1)式に従って時間発展をする方程式です。

今回は初期状態を以下のようにしておきます。

(1)式を離散化しないと数値的に解けないので、離散化します。

ある地点\(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次の精度となります。

\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*}

では、これをPythonで解いて1次元のアニメーションを作成しましょう(^^)/

Pythonで解いて1次元のアニメーション

Pythonのコードを載せておきましょう。

%matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

imax=100
ntmax=2500
dt=0.004
kp=1
xmax=10
dx=xmax/imax

def graph(a,f,x_label,y_label,title):
    plt.plot(a,f)
    plt.grid()
    plt.title(title)
    plt.xlabel(x_label)
    plt.ylabel(y_label)
    plt.show()

    
X = np.linspace(0, xmax, imax)

#initial state
temp=[]
for i in range(imax):
    x=dx*i
    if (i<imax/2):
        temp.append(x)
    else:
        temp.append(xmax-x)

graph(X,temp,"x","temperature(K)","Initial state")

fig = plt.figure()
ims=[]
for t in range(0,ntmax):
    for i in range(1,imax-1):
        temp[i]=temp[i]+kp*dt/(dx*dx)*(temp[i+1]-2*temp[i]+temp[i-1])
    if (t%100==0):
        im = plt.plot(X,temp, "r")
        ims.append(im)

plt.grid()
plt.xlabel("x")
plt.ylabel("temperature(K)")   
ims = animation.ArtistAnimation(fig, ims)
#plt.show()