pythonで初期値境界値問題を差分法で解きグラフを表示するシステム

Question

前提・実現したいこと

pythonで初期値境界値問題を差分法で解くシステムを作っています。

区間 (0,1) 上で、拡散方程式　∂u/∂t = D*(∂^2)*u/∂x^2 について,
境界条件はノイマン条件で初期値は u(x,t) = x 。
初期値境界値問題を D = 1.0 として差分法によって解、t = 0 や t = 0.1 など様々な時刻tにおける解 u(x,t) のグラフを描きたいです。

該当のソースコード

import matplotlib.pyplot as plt
import numpy as np

T = 1.0
M = 60000
L = 1.0
N = 100
D = 1.0

h = L/N
τ = T/M
α = D*τ/(h**2)

x = np.zeros(N+2)
u = np.zeros(N+2)
new_u = np.zeros(N+2)

for j in range(1,N,1):
    x[j] = (j - 0.5)*h
    u[j] = x[j]
   
u[0] = u[1]
u[N+1] = u[N]

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)
lines, = ax.plot(x[1:N+1], u[1:N+1])

for k in range(0,M-1, 1):
    for j in range(1,N,1):
        new_u[j] = α*u[j-1] + (1 - 2*α)*u[j] + α*u[j+1]
        
    for j in range(1,N,1):
        u[j] = new_u[j]

    #境界条件
    u[0] = u[1]
    u[N+1] = u[N]

    lines.set_data(x[1:N],u[1:N])
    plt.pause(τ)     
plt.show()

このように、アニメーションとして描くシステムは作ることができました。
この内容を少し変えたら t を指定してグラフが描けると思うのですが、その方法がわかりません。
何日いろいろ試しても全然分からなかったため質問させていただきました。よろしくお願いいたします。

補足情報（FW/ツールのバージョンなど）

anacondaのspyderで作業しています。

Answer 1

「tを指定して」ということなので

lines.set_data(t[0:M], u[1:N])

のように書きたいということだと思います(上の例は明らかに配列サイズが違うのでNGです)。行うべきことは3つです。

1. 変数tを用意する
2. 変数uを2次元化する
3. 計算と描画の処理を分離する

変数tを用意する

変数x同様、tを定義して下さい。

t = np.zeros(M)

tは時間方向のインデックスですので、計算ループの中で

t[k] = k * τ

とでもしておきます。

変数uを2次元化する

ソースコードを確認するとu(t, x)という関数は

u = np.zeros(N+2)

という一次元配列で表現されています。Nは空間方向の分割数と思われますので、時間の情報が保存されていません。数式ではtとxを持っているのに、プログラム上ではxしか持っていないのです。そこで、まず時間の情報を含められるように、二次元配列にする必要があります。

u = np.zeros(shape=(M, N+2))

uへの代入部分は初期化の部分はu[0, j] = x[j]のように、時間の部分に0を入れ、ループの中ではu[k, j]とループの時間に対するインデックスkを用います。またnew_uは次のステップのために用いる値でしたが、全ての計算結果を保存するようになりましたので不要です。new_u[j]の代わりにu[k+1, j]として計算結果を代入します。

計算と描画の処理を分離する

ソースコードでは「計算しつつ都度描画する」ということを行っていますが、時間発展の計算(ある時点での計算結果がないと次のステップの計算ができない)状態ですので、最後まで計算しないと横軸をtとするグラフを描けません。したがって

計算を行うループと描画を行うループ(x方向に描画を切り替えるループ)を別々にする

ことが必要になります。

Comments

[yuuui]

ご回答ありがとうございます。
私の説明が下手で申し訳なかったのですが、例えば「 t = 0.1 の時の u(x,0.1) のグラフを(縦軸は u 、横軸は x )」を描きたいです。
そうしたい場合はどのようにすればよいのでしょうか。

[yymmt]

t=0.1はインデックスkでは6000に相当しますので、for文のfor k in range(0, M-1, 1)の部分をfor k in range(0, 6000, 1)と変えてループを6000で止めて下さい。そうすると変数uにt=0.1に相当する値が残ります。その後lines.set_dataとplt.pauseの行は削除して、plt.show()の直前にインデントなしで
ax.plot(x[1 : N + 1], u[1 : N + 1])
を追記すれば良いです。

[yuuui]

ご返信ありがとうございます。
試してみたところ、それっぽいグラフは出たのですが、今までアニメーションで見えていたグラフに、x= 1 の端から直線が付け足されたようなグラフになりました。この直線が付いているのが正しいグラフなのでしょうか。
---
import matplotlib.pyplot as plt
import numpy as np

T = 1.0
M = 60000
L = 1.0
N = 100
D = 1.0

h = L/N
τ = T/M
α = D*τ/(h**2)

x = np.zeros(N+2)
u = np.zeros(N+2)
new_u = np.zeros(N+2)

for j in range(1,N,1):
x[j] = (j - 0.5)*h
u[j] = x[j]

u[0] = u[1]
u[N+1] = u[N]

fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

for k in range(0,600, 1):
for j in range(1,N,1):
new_u[j] = α*u[j-1] + (1 - 2*α)*u[j] + α*u[j+1]

for j in range(1,N,1):
u[j] = new_u[j]

#境界条件
u[0] = u[1]
u[N+1] = u[N]

ax.plot(x[1 : N + 1], u[1 : N + 1])

plt.show()
---
どこかにミスがありますでしょうか？何度も申し訳ございません。

[yymmt]

xの配列がN+2個要素があるにもかかわらず、for文でx[j] = (j-0.5) * hではN回しか処理していませんので、最後の2つがゼロのままです。変数 x は以後更新されませんので値が戻ってしまいます。対処方法としては最後の2点を描画しない(ax.plot(x[1:N], u[1:N]とする)か、あるいはx, uの初期化用のforループの後にx[N] = (N - 0.5) * h, X[N+1] = (N+1 - 0.5)*hとして最後まで値を含めることでしょうか。

[yuuui]

なるほど、確かにxの要素が最後0のままでした。
グラフもきちんと描けました。
ありがとうございます。