【プログラミング未経験者】Pythonにて時間ループの中で境界条件をif文を使って時間指定で切り替えたい

import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np

# input parameter
den = 8880.0
cp = 386.0
cond = 398.0
temp_bc = 100.0
temp_init = 0.0
lx = 1.0
nx = 101
tend = 20000.0
dt = 0.1
tout = 100.0

alpha = cond / (den * cp)
dx = lx / (nx - 1)
nt = int(tend / dt)
nout = int(tout / dt)

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

temp_new = np.zeros(nx)

# Boundary condition
temp[0] = temp_bc # Dirichlet @ x=0
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

# time loop
for n in range(1, nt+1):
    # FTCS
    for i in range(1, nx-1):
        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_bc # Dirichlet @ x=0
    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()
``````ここに言語を入力
コード
```### 前提・実現したいこと
プログラミング未経験の者です。以下リンク、科学技術計算講座3-熱伝導方程式シミュレーション-FTCS法のプログラミングのなかで
[リンク内容](http://cattech-lab.com/science-tools/simulation-lecture-3-4/)
指定したある時間幅だけに100℃の熱を加えるモデルを考えたとき、Pythonのプログラムはどのように設定するべきか教えていただきたいです。
具体的に、以下に示したソースコードでは、初期条件で表面温度（ここでは一次元的に考えてx=0の温度）temp_bc = 100.0で固定、計算期間が20000秒で0秒から0.1秒刻みで計算されています。計算期間をもっと短くしてある指定した時間にだけ表面に100℃の熱を与えたときの温度分布を出したいです。

サイトに問い合わせたところ特別に概要のみ以下のように教えて頂けました。

>　プログラムの時間ループ（# time loop）のなかで、
>   # Boundary condition
>   temp[0] = temp_bc # Dirichlet @ x=0  の所の境界条件をif文を使い時間指定で切り替えればよいと思います。


簡単にできるとのことでしたが、プログラミング経験が全くないためif文について調べてもいまいち時間指定での切り替えまでは理解できませんでした。


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


### 該当のソースコード



### 試したこと

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

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

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

y_waiwai

2021/07/13 05:35

このままではコードが読めないので、質問を編集し、<code>ボタンを押し、出てくる’’’の枠の中にコードを貼り付けてください

K_3578

2021/07/13 07:27

とりあえず既に指摘されているソースコードをMarkDownのコードブロックに入れる事ぐらいはしてください。わからないなら調べましょう。

ab-3

2021/07/13 07:46

ご指摘ありがとうございます。修正を行いましたが、また何か不備等あればご指定いただけると幸いです。

行動規範の内容に同意します

回答2件

ベストアンサー

やってみましたが、これで正しいかわかりません。
「計算期間をもっと短くしてある指定した時間にだけ表面に100℃の熱を与えたときの温度分布」という文章が具体的ではないからです。
（計算期間とはシミュレーションの開始から終了までの時間なのか、シミュレーションの制度を上げたいのか。指定した時間とは、開始時点から熱を加える時間なのか、途中からある期間熱を加えるのか。また、100度の熱を加えていない間はどのような熱を加えているのか。など）
プログラムは「うまいこと」はやってくれないため、やりたいことを具体的にするのはどうあがいても人間のしごとであり、プログラミングする人だけではなく、プログラムを人にお願いする人にも必須のスキルです。

python
1import matplotlib.pyplot as plt
2import matplotlib.animation as animation
3import numpy as np
4
5# input parameter
6den = 8880.0
7cp = 386.0
8cond = 398.0
9temp_bc = 100.0
10temp_init = 0.0
11lx = 1.0
12nx = 101
13tend = 20000.0
14
15# ↓ 変更箇所1：計算時間を短くする
16
17# dt = 0.1
18dt = 0.05
19
20# ↑ 変更箇所1：計算時間を短くする
21
22# ↓ 変更箇所２：「ある指定した時間にだけ」というご要求
23
24t_intv = 10.0 #10.0秒
25
26# ↑ 変更箇所２：「ある指定した時間にだけ」というご要求
27 
28tout = 100.0
29
30alpha = cond / (den * cp)
31dx = lx / (nx - 1)
32nt = int(tend / dt)
33nout = int(tout / dt)
34
35#initial condition
36temp = np.full(nx, temp_init) 
37time = 0.0
38
39temp_new = np.zeros(nx)
40
41# Boundary condition
42temp[0] = temp_bc # Dirichlet @ x=0
43temp[nx-1] = temp[nx-2] # Neumann @ x=Lx
44
45# graph data array
46ims = []
47fig = plt.figure()
48ax = fig.add_subplot(1, 1, 1)
49
50gx = np.zeros(nx)
51for i in range(nx):
52    gx[i] = i * dx
53
54# time loop
55for n in range(1, nt+1):
56    # FTCS
57    for i in range(1, nx-1):
58        temp_new[i] = temp[i] + dt * alpha * (temp[i+1] - 2.0 * temp[i] + temp[i-1]) / (dx * dx)
59
60    # update
61    for i in range(1, nx-1):
62        temp[i] = temp_new[i]
63
64    # Boundary condition
65    # ↓ 変更箇所3：「ある指定した時間にだけ」というご要求
66    if time < t_intv :
67        temp[0] = temp_bc # Dirichlet @ x=0
68    else :
69        temp[0] = 0.0
70    # ↓ 変更箇所3：「ある指定した時間にだけ」というご要求
71
72    temp[nx-1] = temp[nx-2] # Neumann @ x=Lx
73
74    time += dt
75
76    if n % nout == 0:   
77        print('n: {0:7d}, time: {1:8.1f}, temp: {2:10.6f}'.format(n, time, temp[nx-1]))
78        im_line = ax.plot(gx, temp, 'b')
79        im_time = ax.text(0, 110, 'Time = {0:8.1f} [s]'.format(time))
80        ims.append(im_line + [im_time])
81
82# graph plot
83ax.set_xlabel('x [m]')
84ax.set_ylabel('Temperature [C]')
85ax.grid()
86anm = animation.ArtistAnimation(fig, ims, interval=50)
87anm.save('animation.gif', writer='pillow')
88plt.show()
89
90

投稿2021/07/13 06:32

tain

総合スコア241

ab-3

2021/07/13 07:24

ご丁寧にご回答ありがとうございます。質問内容が不明確で申し訳ございません。ここでの計算期間ですが、シミュレーション開始から終了までの時間tend = 20000.0を示しています。また、指定した時間とは途中からある期間だけ表面に100℃の熱を加えることを示しています。100度の熱を加えていない間は初期パラメータの中でtemp_init = 0.0に示される通り0℃と考えております。例として、計算期間を20000秒から30秒に短縮しその中の10~12秒の間だけ、もともと全体が0℃の試料に対して表面に100℃の熱を加えたときの温度分布を求めたいということになります。修正していただいたコードをAnaconda Promptにて計算しグラフをみたところ、表面位置x=0が100℃にならず、0℃固定のままで0.2mのところに瞬間的に2.5℃程度のピークが出力されました。

行動規範の内容に同意します