\(\require{physics}\)

拡散方程式の計算

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from matplotlib import animation
plt.rcParams['font.family'] = 'Times New Roman'

# setup
L  = 1.0         # 領域の長さ
n  = 22          # CV(セル)の数,BC2個分を含む
ts = 200         # 全ステップ数(ts*dtは全時間)
dx = L/(n-2)     # セルの幅
dt = 1000        # 時間間隔(秒)
 temp = [20.0]*n  # 温度用配列の初期化 
K = 0.14*1.0E-06 # 熱拡散率

x = [(i-0.5)*dx for i in range(n)] # セル中心の登録
for i in [9,10,11,12]: # 初期条件
  temp[i]=50.0
temp_new = temp[:] # 次ステップの温度用配列
fig = plt.figure()
img0, = plt.plot(x,temp,"k:") # 初期値のプロット
imgs=[[img0,img0]]
# 時間進行
for m in range(ts): # 時間のループ
  for i in range(1, n-1):   # 空間のスイープ
    temp_new[i]=temp[i]+K*dt/dx**2 \
    *(temp[i+1]-2*temp[i]+temp[i-1])
  # 境界条件の設定(ノイマン条件)
  temp_new[0]   = temp_new[1]
  temp_new[n-1] = temp_new[n-2]
  temp=temp_new # 未来の値を現在の値に
  if m%10==0:img, = plt.plot(x,temp,"k", lw=2) #  10ステップ毎に表示
  plt.gca().set_xlim(0,1.0)
  plt.gca().set_ylim(10,60)
  plt.xlabel("$x$ (m)")
  plt.ylabel("Temperature (deg)")
  imgs.append([img0,img])

ani = animation.ArtistAnimation(fig, imgs, interval=20, blit=True,repeat_delay=1000)
ani.save('diffusion.gif', writer="imagemagick")