Python で1次元の熱方程式を2通りの方法で解いてみた

両端が断熱的な導線の熱伝導を記述する, 次の斉次 Neumann 境界付きの熱方程式の初期値境界値問題

\begin{equation}
\left\{
\begin{split}
&\partial_t u = \partial_x^2 u,\quad &t>0 , x\in (0,\pi),\\
&\partial_x  u(t,0) =\partial_x u(t,\pi) = 0 ,\quad& t>0,\\
&u(0,x) = f(x),\quad & x\in (0,\pi)
\end{split}
\right.
\end{equation}
を考えます*1.
ここで,  \partial_t  =\frac{\partial}{ \partial t},  \partial_x  =\frac{\partial}{ \partial x} を表します.
前進差分法とスペクトル法を用いる方法の2通りでの数値計算例のデモを作成しました. *2
コードがメインで, 数値解析の手法については多くを説明しません.
意味が通らない箇所があれば気軽にコメントをよろしくお願いします.
いずれの場合も初期値は

f(x)=
\begin{cases}
1, \quad \textrm{if $0 < x \le \frac{\pi}{2}$},\\
\cos{4x}, \quad \textrm{if $\frac{\pi}{2} < x < \pi$}\\
\end{cases}
により与え, 時刻  t  0 < t < 1 の範囲で解きます.

import matplotlib.pyplot as plt
#line domain
xmin = 0.0
xmax = np.pi
nx = 64
dx = (xmax-xmin)/nx

#time interval
tmin = 0.0
tmax = 1.0
nt = 1024
dt = (tmax-tmin)/nt

c = 1.0 #spreading coefficient

M = 51 #maximam fourier number

#initialize
F = np.zeros(nx)

T = np.linspace(tmin, tmax, nt, endpoint = False)
X = np.linspace(0, np.pi, nx, endpoint = False)+dx/2
F = np.cos(4*X)
for i in range(nx//2):
    F[i] = 1

1. 前進差分法で解こう

ここでは空間メッシュ数を  n_x =64 刻み幅を  h =\pi/n_x とし,
時間メッシュ数を  n_t = 1024 , 時間刻み幅を  \tau =\frac{1}{n_t} とします.
熱方程式は解  u がある程度滑らかだという仮定のもと, 時間微分を前進差分, 空間微分を中心差分に置き換えることで,

\frac{u(t+\tau,x) - u(t,x)}{\tau} = \frac{u(t,x+h)-2u(t,x)+u(t,x-h)}{h^2}
と離散化されます.

rsq = c*dt/(dx**2) #constant

#solution
v = np.zeros((nt,nx))

#initial data
v[0] = F
print(len(F))

#simulation
for t in range(nt-1):
    for x in range(1,nx-1):
        v[t+1,x] = (1-2*rsq)*v[t,x]+rsq*(v[t,x-1]+v[t,x+1])
    #Neumann condition
    v[t+1,0] = v[t+1,1]
    v[t+1,nx-1] = v[t+1,nx-2]

これを可視化すると次のような動画が得られます.

fig1 = plt.figure()
fig1.set_dpi(100)
ax = fig1.add_subplot(1,1,1)
def animate(t):
    ax.clear()
    plt.ylim([-1.2,1.5])
    
    p, = plt.plot(X,v[t])
    time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
    time_text.set_text('t = %.2f/%.2f,' 
                       % ((t+1)/nt,tmax))
    
anim = animation.FuncAnimation(fig1,animate,frames=nt,interval=60,repeat=False)
plt.show()

f:id:wakabame:20180604015116g:plain

2.スペクトル法で解こう

記事
trtmfile.com
で導入したフーリエ変換を用いて数値解を求めます. まずは初期値をフーリエ変換しましょう.

F_M = [np.zeros(nx) for m in range(M)]

#Fourier coefficients
A = [0]*M

for m in range(M):
    A[m] = (np.dot(np.cos(m*X),F))/nx

F_M[0] = np.cos(0*X)*A[0]
for m in range(1,M):
    F_M[m]=2*np.cos(m*X)*A[m]
fig1 = plt.figure()
fig1.set_dpi(100)
ax = fig1.add_subplot(1,1,1)
def animate(t):
    ax.clear()
    plt.ylim([-1.2,1.5])
    
    p, = plt.plot(X,F_M[t])
    time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
    time_text.set_text('t = %.2d/%.2d,' 
                       % (t,M-1))
    
anim = animation.FuncAnimation(fig1,animate,frames=M,interval=600,repeat=False)
plt.show()

f:id:wakabame:20180604015618g:plain
端点を含めて滑らかな関数なので, ギブス現象も起こらずによく近似されています.
スペクトル法を用いて解を求めましょう.

u = [[[0 for x in range(nx)] for m in range(M)] for t in range(nt)]  #u[t][m][x]
u[0]= F_M
for m in range(M):
    for t in range(nt):
        u[t][m] = np.dot(np.exp(-t*dt*m*m),u[0][m])
for m in range(1,M):
    for t in range(nt):
        u[t][m] += u[t][m-1]

これを可視化すると次のような動画が得られます.

fig1 = plt.figure()
fig1.set_dpi(100)
ax = fig1.add_subplot(1,1,1)
def animate(t):
    ax.clear()
    plt.ylim([-1.2,1.5])
    
    p, = plt.plot(X,u[10*t][M-1])
    time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
    time_text.set_text('t = %.2f/%.2f,' 
                       % ((t+1)/nt,tmax))
    
anim = animation.FuncAnimation(fig1,animate,frames=nt,interval=60,repeat=False)
plt.show()

f:id:wakabame:20180604033126g:plain

3. それぞれの解を比較しよう.

fig1 = plt.figure()
fig1.set_dpi(100)
ax = fig1.add_subplot(1,1,1)
def animate(t):
    ax.clear()
    plt.ylim([-1.2,1.5])
    
    p, = plt.plot(X,v[10*t])
    p, = plt.plot(X,u[10*t][M-1])
    time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
    time_text.set_text('t = %.2f/%.2f,' 
                       % ((t+1)/nt,tmax))
    
anim = animation.FuncAnimation(fig1,animate,frames=nt,interval=60,repeat=False)
plt.show()

f:id:wakabame:20180604033200g:plain
2つの数値解はそこそこ近いものになっていますね!
教えてくれたけますさん, ありがとうございます.
コードはこちらにございます.
github.com

*1: u の第1引数を時間変数, 第2引数を空間変数としています.

*2:有限要素法でもやりたかった