7.1 陽解法

最初の素朴なバージョン。全部計算して，描画していって，最後に一気に表示する。

heat1d-v0.py

#!/usr/bin/env python3
# heat1d-v0.py

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    y=x.copy()
    for i in range(0,len(y)):
        if y[i] > 0.5:
            y[i] = 1-y[i]
    return y

N=50

x=np.linspace(0.0, 1.0, N+1)
u=f(x)
newu=np.zeros(N+1)

h=1.0/N
lam=0.5  # lambda は予約語で使えない？
tau=lam*h*h
dt=0.01
skip=int(dt/tau)

Tmax=1.0
nmax=int(Tmax/tau)

plt.plot(x,u)

for n in range(1,nmax):
    for i in range(1,N):
        newu[i]=(1-2*lam)*u[i]+lam*(u[i-1]+u[i+1])
    if n%skip == 0:
        plt.plot(x,newu)
    u=newu.copy()

plt.show()

少し改善したバージョン。 初期値 $f$ の計算に numpy の vectorize() を使うとか， 内側の for を取って newu[] を省略するとか， 計算しながら表示するとか。

heat1d-e.py

#!/usr/bin/env python3
# heat1d-e.py --- 1次元熱方程式
# http://d.hatena.ne.jp/Megumi221/20080306/1204770689 を参考にした

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return min(x,1-x)

N=50
x=np.linspace(0.0, 1.0, N+1)
# vf=np.vectorize(f)
# u=vf(x)
u=np.vectorize(f)(x)

h=1.0/N
lam=0.5  # lambda は予約語で使えない？
tau=lam*h*h
Tmax=1.0
nmax=int(Tmax/tau)
dt=0.001
skip=int(dt/tau)

plt.ion()
line,=plt.plot(x,u)

for n in range(1,nmax):
    u[1:N]=(1-2*lam)*u[1:N]+lam*(u[0:N-1]+u[2:N+1])
    if n%skip == 0:
        line.set_ydata(u)
        plt.pause(0.00001) # 除くと動かない。必要な理由を理解できてない