5.1 陽解法

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

heat1d-v0.py

#!/opt/local/bin/python2
# -*- coding: utf-8 -*-
# heat1d-v0.py

import numpy as np
import matplotlib.pyplot as plot
import matplotlib.animation as animation

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)

plot.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:
        plot.plot(x,newu)
    u=newu.copy()

plot.show()

少し改善したバージョン。初期値 $f$ の計算に numpy の vectorize() を使うとか，newu[] を省略するとか， 計算しながら表示するとか，range() の代わりに xrange() を使うとか。

heat1d-e.py

#!/opt/local/bin/python2
# -*- coding: utf-8 -*-
# heat1d-e.py
# http://d.hatena.ne.jp/Megumi221/20080306/1204770689 を参考にした

import numpy as np
import matplotlib.pyplot as plot

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)

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

for n in xrange(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)
    plot.draw()

(細かい話ですが，某マシンで range() では84.822秒， xrange() ではユーザー時間$80.218$ 秒, 経過時間$2$ 分$38.58$ 秒。)