別配列同士の足し算を同じforループ内で行えるのか

Question

・目的
Pythonを使って一次元非定常熱伝導方程式を差分法で表したFTCS式で計算し、0℃の試料（深さ方向iの1次元のみ考える）に一定時間レーザー光を照射したときの温度分布をグラフに示したいと思っています。
・問題
コード中の（#レーザー吸収熱温度 Qの設定）で得た計算Qの配列を＃FTCS式における別のtempという配列に足したいのですが、以下のエラーが表示され計算できませんでした。
(エラー内容）
＃FTCS内のemp_new[i] = Q[i]+temp[i] +・・・式にてIndexError:list index out of rangeのエラー。Qの配列とtempの配列を同じ要素iでforループを回して#FTCS内の式で計算するにはどのように設定するべきでしょうか？

・FTCS式について
forループのiですが、試料を表面i＝0から裏面に向かって深さ方向（一次元）に101等分した要素番号としています。表面i＝0と裏面i=nxは境界条件を指定しているためforループに入れていません。FTCS式では、左辺（位置iにおける時間n+1の温度）＝右辺（位置iおよびi-1、i+1の時間nの温度）で計算されています。時間については＃timeループで時間の要素nについて回しています。
FTCS式では右辺は全て一つ前の時間のデータで計算する必要があるため、左辺をtemp_new[i]と別の配列に格納しています。FTCS法のforループが終わると、先程格納した temp_new を temp の配列にコピーして temp を更新します。＃updateでforループを回し temp_new の各要素を temp の同じ番号の要素に代入しているような状態です。
今回はFTCS式の右辺にtempとは別の配列ですが、＃レーザー光照射による吸収熱の設定、で与えた配列の中で位置iにおけるQ[i]を足したいと思っています。

import matplotlib.animation as animation
import numpy as np

# input parameter(材料:Si、厚さ1㎜、微小間隔 dx：10㎛、計算期間tend:1μs、時間間隔dt:1ns、初期温度temp_init:0℃、空間分割数nx、時間分割数nt）
den = 2330.0
cp = 678.0
cond = 83.7
temp_init = 0.0
lx = 0.001
nx = 101
tend = 0.000001
dt =  0.000000001
tout =0.000001
a=0.064
alpha = cond / (den * cp)
dx = lx / (nx - 1)
nt = int(tend / dt)
nout = int(tout / dt)


# レーザー光照射による吸収熱温度の設定
Q=[2000000000000*dt*pow(2.718281828,-a*dx*i) for i in range(1, nx-1)]

#initial condition Temp
temp = np.full(nx, temp_init) 
time = 0.0

temp_new = np.zeros(nx)

# Boundary condition
temp[0] = temp[1] 
temp[nx-1] = temp[nx-2] # Neumann @ x=Lx

# graph data array
ims = []
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

gx = np.zeros(nx)
for i in range(nx):
    gx[i] = i * dx

# レーザー光照射時間の指定
t_intv= 0.00000005

# time loop
for n in range(1, nt+1):

    # FTCS(右辺第一項に吸収熱温度を設定/0～50nsの間だけ照射、それ以降はQ＝0)
    for i in range(1, nx-1):
        if time < t_intv :
            temp_new[i] = Q[i]+temp[i] + dt * alpha * (temp[i+1] - 2.0 * temp[i] + temp[i-1]) / (dx * dx)
        else :
            temp_new[i] = temp[i] + dt * alpha * (temp[i+1] - 2.0 * temp[i] + temp[i-1]) / (dx * dx)
             
    # update
    for i in range(1, nx-1):
        temp[i] = temp_new[i]

   
    # Boundary condition
    temp[0] = temp[1] 
    temp[nx-1] = temp[nx-2] # Neumann @ x=Lx

    time += dt

    if n % nout == 0:   
        print('n: {0:7d}, time: {1:8.1f}, temp: {2:10.6f}'.format(n, time, temp[nx-1]))
        im_line = ax.plot(gx, temp, 'b')
        im_time = ax.text(0, 110, 'Time = {0:8.1f} [s]'.format(time))
        ims.append(im_line + [im_time])  

# graph plot
ax.set_xlabel('x [m]')
ax.set_ylabel('Temperature [C]')
ax.grid()
anm = animation.ArtistAnimation(fig, ims, interval=50)
anm.save('animation.gif', writer='pillow')
plt.show()

Comments

[meg_]

> 今回の場合、おそらくQ[i]がエラーの要因かと思われます。Qとtempは別の配列ですが、どのようにすればエラーを解決してtempにQ[i]を足した計算ができますでしょうか？
Q[i]でIndexError:list index out of rangeが発生したのであれば(エラー発生時の)Q[i]は存在しないので「tempにQ[i]を足した計算」も不可能でしょう。

[ab-3]

ご指摘ありがとうございます。修正致しました。また何かおかしい所などあればお願い致します。

Answer 1

そもそもQのインデックスとQ[i]のiの範囲が別の変数ですから、範囲が合っているかどうかも不明ですね。

nx = 101
Q=[2000000000000*dt*pow(2.718281828,-a*dx*i) for i in range(1, nx-1)]

tend = 0.000001
dt =  0.000000001
nt = int(tend / dt)

for n in range(1, nt+1):
            temp_new[i] = Q[i]+temp[i] + dt * alpha * (temp[i+1] - 

どうでしょう？

Comments

[ab-3]

回答ありがとうございます。この場合、temp中のiにもQと同様にforループがかかっているのでしょうか？