Pythonによる拡散方程式の数値計算

下記の拡散方程式をPythonで数値計算します．

$$
\cfrac{\partial f}{\partial t} = D \cfrac{\partial^2 f}{\partial x^2}
\tag{1}
$$

ここで, $f$は拡散する量，$x$は流下方向，$t$は時間, $D$は拡散係数です．
$x$軸方向に$\Delta x$で区切って上流から下流に向かって$1,2,3,....,n_x$と番号をつけ，ある時刻$t$における$i$番目の$f$の値を$f(i,j)$とします．
ただし，$j$は時間方向の番号で, 時刻が$t$のときは$j$, $t+\Delta t$の時は$j+1$となります．

拡散方程式直観的に中央差分が最適で，差分式は次式で表わされます.
$$
f(i,j+1)=f(i,j)+\cfrac{D \Delta t}{{\Delta x}^2}
\left[ f(i+1,j)-2f(i,j)+f(i-1,j) \right] \tag{2}
$$

ただ，(2)式のような差分式にすると$(i,j)$の2次元配列になり，特に$j$方向の配列が$\Delta t$毎に必要になるので，使用するメモリーがどんどん大きくなってしまいます．
そこで，メモリー節約のために，時間方向$j$方向の古い値は直前のもの以外は必要ないのでどんどん捨てて行く方式を採用し，下記の差分式とします．

$$
f_{new}(i+1)=f_{old}(i)+\cfrac{D \Delta t}{{\Delta x}^2}
\left[ f_{old}(i+1)-2f_{old}(i)+f_{old}(i-1) \right] \tag{3}
$$

この(3)式を$i=1 \sim n_x$で計算し，その後に

$$
f_{old}(i)=f_{new}(i) \tag{4}
$$

を$i=1 \sim n_x $で計算し，(3)式と(4)式を計算後に時刻を$\Delta t$だけ進めていくと，1次元のメモリーで時間発展の計算をどんどん進めることが可能となります．

ではこれをPythonで計算してみましょう．
初期条件は以前にやった移流方程式の時と同じく三角形分布としますが，拡散方程式の場合は影響が上流・下流の両方に及ぶので，初期三角形分布の頂点は対象領域の中央に置きます．また，拡散係数は任意に与えることが可能ですが，この例では$D=0.01$とします．

# coding: shift_jis

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from matplotlib.animation import PillowWriter
from matplotlib._version import get_versions as mplv
from scipy.stats import gaussian_kde

fig=plt.figure()

xl=100; nx=201; d=0.01
xb=5;xp=50;fp=0.5
etime=1000;dt=5

x=np.linspace(0,xl,nx)
dx=x[2]-x[1]
f=[0]*nx; fn=[0]*nx
ims=[]

for i in range(0,nx):
    if(x[i]>=(xp-xb) and x[i]<=xp):
        f[i]=fp/xb*(x[i]-xp+xb)
    elif(x[i]>xp and x[i]<(xp+xb)):
        f[i]=fp/xb*(xp+xb-x[i])
    else:
        f[i]=0

im=plt.title("Diffusion")
im=plt.xlabel("x")
im=plt.ylabel("f")

time=0
while(time<=etime):
    for i in range(1,nx-1):
        fn[i]=f[i]+d*(f[i+1]-2*f[i]+f[i-1])*dt/dx**2

    f=fn
    time=time+dt
    im=plt.plot(x,f,'r');
    ims.append(im)

ani = animation.ArtistAnimation(fig, ims, interval=10)
#plt.show()
ani.save('diffusion.gif',writer='imagemagick')
#ani.save('diffusion.mp4',writer='ffmpeg')

このコードを実行すると出来るgifアニメはこちらです．
https://i-ric.org/yasu/qiita/diffusion.gif"