8.1 陽解法

まず全部計算してから一気にグラフを描く。

#!/opt/local/bin/python3.5
# -*- 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()

計算しながらグラフを描く等、少し工夫する。

#!/opt/local/bin/python3.5
# -*- 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 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)
        plot.pause(0.0001)