Pythonにて配列の計算をしたい

Question

・目的
Pythonを使って一次元非定常熱伝導方程式を差分法で表したFTCS式で計算し、0℃の試料（深さ方向の1次元のみ考える）に一定時間レーザー光を照射したときの温度分布をグラフに示したいと思っています。
・問題
コード中の# レーザー光照射による吸収熱の設定で得た計算Qの配列を＃FTCS式の配列に足したいのですが、複数のエラーが表示され計算できませんでした。
(エラー内容）
TypeError: only size-1 arrays can be converted to Python scalars
The above exception was the direct cause of the following exception:
Traceback (most recent call last):
File "Q.py", line 55, in <module>
temp_new[i] = Q+temp[i] + dt * alpha * (temp[i+1] - 2.0 * temp[i] + temp[i-1]) / (dx * dx)
ValueError: setting an array element with a sequence.

最後のValueErrorに関しては要素の種類が異なる場合に起こるとありましたが、具体的な修正箇所,方法までは理解できませんでした。
また何かおかしいところなどあれば教えていただきたいです。

・FTCS式について補足
FTCS式では右辺は全て一つ前の時間のデータで計算する必要があるため、左辺をtemp_new[i]と別の配列に格納しています。
FTCS法のforループが終わると、先程格納した temp_new を temp の配列にコピーして temp を更新します。＃updateでforループを回し temp_new の各要素を temp の同じ番号の要素に代入しているような状態です。

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

[norapomu]

出たエラーの内容も掲載しましょう
どういうエラーで困っているのかがわかりません

[ab-3]

ご指摘ありがとうございました。

Answer 1

• 初歩的な質問で申し訳ないのですが、今回のように別々の配列同士を足すことは可能なのでしょうか？

可能です。

• また何かおかしいところなどあれば教えていただきたいです。

配列の要素temp_new[i]はfloatを一個しか入れることができないのに、右辺の計算結果であるfloatが101個入っている配列を入れようとしているからです。

Comments

[ab-3]

回答ありがとうございます。右辺に関して、tempと同様にQ[i]とすればよいでしょうか？101個のfloatが入ったQの配列からtemp[i]に合わせて一つのfloatを抜き出すにはどのようにコードで指定するべきでしょうか。知識不足で申し訳ございません、回答いただけるとありがたいです。

[ppaul]

そのご質問の答えは、何をやりたいかによります。
tempと同様にQ[i]とすればよいでしょうか？に答えるためには、
ab-3さんがやりたいことを理解して、ab-3さんの書いたコードを理解しなければならず、それには情報がたりません。
やりたいことを、Pythonコードではなく、数式なり文章で表現できて、それを質問に載せれば、誰かが回答できるかもしれません。

[ab-3]

回答ありがとうございます。forループのiですが、試料を表面i＝0から裏面に向かって深さ方向(1次元)に101等分した要素番号としています。表面i＝0は境界条件を指定しているためforループはそのひとつ奥のi＝1から始まるようにしています。FTCS式では左辺(位置iにおける時刻n+1の温度)＝右辺(位置iおよびi−1.i+1の時刻nの温度)で計算されており、次の時刻n+1の位置iにおける温度は1つ前の時刻nの位置iとその両隣りi+1、i−1の温度によって決まるということです。時刻については#timeloopで時間の要素数nに対してforループを回しています。ここで自分が行いたい事は深さ方向iに対して指数関数的に減少するレーザー光の吸収熱温度Qを設定しFTCSの右辺に足したいということになります。左辺(次の時刻n+1の位置iの温度)を求めるためには1つ前の時刻nの位置iにおけるQ[i]を取り出して右辺に足すべきかと考えています。長文失礼致しました。何か1つでもご理解頂けるようでしたら回答頂けると助かります。

[ppaul]

左辺(次の時刻n+1の位置iの温度)を求めるためには1つ前の時刻nの位置iにおけるQ[i]を取り出して右辺に足すべきかと考えています。

1つ前の時刻nの位置iにおけるQ[i]を取り出して右辺に足す、と決めれば、「右辺に関して、tempと同様にQ[i]とすればよい」のでしょう。
「右辺に足すべきかと考えています。」を誰かに決めてください、というのであれば、プログラミングの問題ではなく熱伝導方程式自体を教えてほしいと言っているように見えます。その場合は、teratailではなく物理学の質問として別の場で質問される方がよろしいかと思います。

[ab-3]

回答ありがとうございます。すみません、こちらの表現が曖昧でした。意味としては最初ににおっしゃって頂いた通り、1つ前の時刻nの位置iにおけるレーザ光吸収熱による温度Q[i]を取り出して右辺に足すということでした。更新したコード通り、FTCS式の右辺に足すQをQ[i］として計算しましたが、IndexError: list index out of rangeが表示されました。リストのサイズを超えたインデックスの要素を参照しようとした時に発生するとのことですが、今回の場合どこが問題なのでしょうか？
度々すみません、回答いただけるとありがたいです。

[ppaul]

新しいエラーメッセージはどの行で出ていますか？

[ab-3]

＃FTCS内の temp_new[i] = Q[i]+temp[i] ・・・式で出ています。

[ab-3]

すみません、あれからいろいろ調べてみたのですがよく分かりませんでした…エラーの要因について何か分かりましたでしょうか？