拡散方程式の数値とグラフ

Question

前提

Pythonで１次元熱拡散方程式の数値計算を中心差分法で行い、グラフを書くプログラムを作っています。
エラーと正しい数値が出ない原因がわかりません。
初期条件はT(0,l)=0(0≦l<L)
1(l=L)
境界条件はT(t,0)=0, T(t,L)=1

実現したいこと

発生したエラーの解消とグラフを１本(t)ごとに色が変わるようにしたいです。
初期条件、境界条件が正しく出るようにしたいです。

発生している問題・エラーメッセージ

Traceback (most recent call last):
File "C:\Users\kucoo\Downloads\3班_配布物\1214.py", line 18, in <module>
T[i][N]=1.0
IndexError: index 11 is out of bounds for axis 0 with size 11

該当のソースコード

ソースコード
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
M=41
N=11
tend=20.0
L=1.0
k=0.01
dt=0.5
dl=0.1
T=np.zeros((M,N))
a=k*dt/(dl*dl)
for j in range(N-1):
  T[0][j]=0.0
T[0][N-1]=1.0
for i in range(M):
  T[i][0]=0.0
  T[i][N]=1.0
for i in range(M-1):
  for j in range(N-1):
     T[i+1][j]=T[i][j]+a*(T[i][j+1]-2*T[i][j]+T[i][j-1])
     print(i*dt,j*dl,T[i][j])
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.set_xlabel("x")
ax1.set_ylabel("T:Temperature ")
ax1.set_title(r'Temperature')
ax1.set_xlim(0, 1)
ax1.glid()
x = np.arange(0, L+dl, dl)
for t in range(tend):
  y = T[:, int(L/dl)]
  ax1.plot(x, y, label="t: {}".format(f"{x}"))
plt.legend()
plt.savefig("temperature.png")

### 試したこと

ここに問題に対して試したことを記載してください。

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

ここにより詳細な情報を記載してください。

Comments

[melian]

T[i][N-1]=1.0 でしょうね。

Answer 1

質問に掲載されているコードのうち，21行目のfor j in range(N-1):における範囲設定が致命的なエラーとなって正しく計算できていなかった上，31行目に用意した配列xの長さと33行目に用意した配列yの長さの不一致で描画不可です．NやMの意味の理解に苦しみましたが，熱拡散方程式の式4を参考に，上のコードを読み替えてコードを再構築し，グラフ描画しました．

import os
from glob import glob

import numpy as np
from PIL import Image
from matplotlib import pyplot as plt
from matplotlib import cm
from matplotlib.colors import Normalize

tend = 20.0  # time end [s]
length = 1.0 # length [m]
alpha = 0.01 # [m^2/s]
dt = 0.5     # [s]
dl = 0.1     # [m]
s = alpha * dt / (dl * dl)

t = np.arange(0, tend + dt, dt)
L = np.arange(0, length + dl, dl)

T = np.zeros((len(t), len(L))) # celsius degree
# initialization
T[0, :-2] = 0 # T(0, l) = 0 (0≦l<L)
T[0, -1] = 1  # T(0, L) = 1
T[:, 0] = 0   # T(t, 0) = 0
T[:, -1] = 1  # T(t, L) = 1

for n in range(len(t) - 1):
  for i in range(1, len(L) - 1):
     T[n + 1][i] = T[n][i] + s * (T[n][i + 1] - 2 * T[n][i] + T[n][i - 1])

norm = Normalize(0, tend)
plt.figure(figsize = (8, 6))
for i, _t in enumerate(t):
  plt.plot(L, T[i], color = cm.jet(norm(_t)))
plt.xlim(-0.01, 1.01)
plt.ylim(-0.01, 1.01)
plt.xlabel("x[m]")
plt.ylabel("T:Temperature")
plt.grid()
plt.colorbar(cm.ScalarMappable(norm = norm, cmap = cm.jet), label = "Time[s]")
plt.savefig(f"temperature.png", dpi = 80)

# 動画にする
os.makedirs("./fig", exist_ok = True)
for file in glob("./fig/*.png"):
  os.remove(file)

for i, _t in enumerate(t):
  plt.figure(figsize = (6, 6))
  plt.plot(L, T[i])
  plt.ylim(-0.01, 1.01)
  plt.xlim(-0.01, 1.01)
  plt.xlabel("x[m]")
  plt.ylabel("T:Temperature")
  plt.title(f"t: {_t} [s]")
  plt.grid()
  plt.savefig(f"./fig/{i:02}.png", dpi = 80)
  plt.close()

images = list(map(lambda file: Image.open(file), sorted(glob("./fig/*.png"))))
images[0].save("out.gif", save_all = True, append_images = images[1:], duration = int(dt * 1000), loop = 1)

maptlotlibは複数描画するときなるべく色が被らないようにしてくれますが，今回41本の線を描こうとすると色が被ったので無理やりcm.jetを利用して色を適用しています．カラーマップは他にもたくさんあるので，自分の好きなものを描画に利用されると良いでしょう．また，時間で推移するものを図示するときは，可能な限り動画で表現されることを推奨します．

Answer 2

エラーメッセージを読みましょう.
読みたくない気持ちもわかりますが, エラーメッセージは有益な情報を教えてくれます.

Traceback (most recent call last):
File "C:\Users\kucoo\Downloads\3班_配布物\1214.py", line 18, in <module>
T[i][N]=1.0
IndexError: index 11 is out of bounds for axis 0 with size 11

つまり, 18行目で指定しているインデックス11はサイズ11の0軸の外側ですという意味です.

もっとかみ砕いて言うと,

N = 11
T[i][N]

ではT[i]の12番目の要素を指定してしまっていることが原因です.

したがって,

T[i][N-1] = 1.0

に直しましょう.

Comments

[uqrt]

ご指摘ありがとうございます。そことその後の中心差分を直したら良いということですね。申し訳ありませんが、グラフに関しての質問にもお答えいただけると幸いです。