pythonで、if文の中にあるwhile文が上手く動作していないので、どうすればよいか教えていただけますでしょうか。

Question

ルンゲクッタ法を用い、化学反応の数値計算を行うプログラムです。
while文の中のif文で、t==2のときにAに値を代入したいのですが、できません。回答よろしくお願いいたします。

#中間体の生成を無視した場合
from sys import getsizeof
import numpy as np
import matplotlib.pyplot as plt

k1 = 1
k2 = 1
k3 = 1
k4 = 1

def reaction_equations(concentrations):
# Define the system of ordinary differential equations
# Modify this function according to your specific reaction system
A, B, C, D, E, F = concentrations
dAdt = - k1 * A * k4 * F
dBdt = (k1 * A) * (k4 * F) - k2 * B * k2 * B * k2 * C
dCdt = - k2 * C * k2 * B * k2 * B
dDdt = k2 * B * k2 * B * k2 * C
dEdt = k2 * B * k2 * B * k2 * C - k3 * E
dFdt = k3 * E - k4 * F * k1 * A

return [dAdt, dBdt, dCdt, dDdt, dEdt, dFdt]

def runge_kutta(dt, t_max, initial_concentrations):
# Perform the simulation using the fourth-order Runge-Kutta method
t = 0.0
concentrations = np.array(initial_concentrations)
time_points = [t]
concentration_points = [concentrations]

while t < t_max:
    if t == 2:
        A += 1.0
    else:
    # Compute concentration changes using the fourth-order Runge-Kutta method
      k1 = reaction_equations(concentrations)
      k2 = reaction_equations(concentrations + 0.5 * dt * np.array(k1))
      k3 = reaction_equations(concentrations + 0.5 * dt * np.array(k2))
      k4 = reaction_equations(concentrations + 0.5 * dt * np.array(k3))

      concentrations += (dt / 6.0) * (2 * np.array(k1) + 2 * np.array(k2) + 2 * np.array(k3) + 2 * np.array(k4))
    t += dt
    
    time_points.append(t)
    concentration_points.append(concentrations.copy())


return time_points, np.array(concentration_points)

Example usage

dt = 1e-3 # Time step size
t_max = 5 # Maximum simulation time
initial_concentrations = [1.0, 0.0, 1.0, 0.0, 0.0, 1.0] # Initial concentrations of species A, B, C, D, E, F, G
time_points, concentration_points = runge_kutta(dt, t_max, initial_concentrations)
print(concentration_points[-1])

Plotting the concentration curves

plt.plot(time_points[1:], concentration_points[1:, 0], label='R1X', color=(200/255, 0/255, 100/255))
plt.plot(time_points[1:], concentration_points[1:, 1], label='R1PdX', color=(100/255, 0/255, 200/255))
plt.plot(time_points[1:], concentration_points[1:, 2], label='R2', color=(0/255, 100/255, 200/255))
plt.plot(time_points[1:], concentration_points[1:, 3], label='R1R2', color=(200/255, 100/255, 0/255))
plt.plot(time_points[1:], concentration_points[1:, 4], label='HPdX', color=(100/255, 200/255, 0/255))
plt.plot(time_points[1:], concentration_points[1:, 5], label='Pd(0)', color=(200/255, 200/255, 0/255))

plt.xlabel('Time')
plt.ylabel('Concentration')

plt.xscale('log')

plt.yscale('log')

plt.legend()
plt.show()

Comments

[melian]

dt = 1e-3 # Time step size

とされていますが、10進数の 1e-3(0.001) は2進数では循環小数になります。そのため、while loop で t の値が 2.0 になることはありません。

　　:
t = 1.9999999999998905
t = 2.0009999999998906
　　:

t の値ではなく、整数のループカウンタを使うほうがよいかと思います。

[winterboum]

pythonはインデント命ですから、<code> を使って 載せ直してください

[winterboum]

今回の修正を元に戻し、
「中間体の生成を無視した場合」 の行の上と 「Example usage」の行の上に ```(逆コーテーション３つ）だけの行を入れて。
んでもって、プレビュー確認してね

Answer 1

コメントにも書きましたが、dt(0.001) は2進数では循環小数になってしまいますので、numpy.arange で予め time_points を作成しておいて、time_points をループカウンタとして使う方がよいかと思います。
※ A はおそらく concentrations[0] ではないでしょうか

def runge_kutta(dt, t_max, initial_concentrations):
    # Perform the simulation using the fourth-order Runge-Kutta method
    t = 0.0
    concentrations = np.array(initial_concentrations)
    time_points = np.arange(t, t_max+dt, dt)
    concentration_points = [concentrations]

    for ti in time_points[:-1]:
        if ti == 2.0:
            concentrations[0] += 1.0
        else:
            # Compute concentration changes using the fourth-order Runge-Kutta method
            k1 = reaction_equations(concentrations)
            k2 = reaction_equations(concentrations + 0.5 * dt * np.array(k1))
            k3 = reaction_equations(concentrations + 0.5 * dt * np.array(k2))
            k4 = reaction_equations(concentrations + 0.5 * dt * np.array(k3))

            concentrations += (dt / 6.0) * (2 * np.array(k1) + 2 * np.array(k2) + 2 * np.array(k3) + 2 * np.array(k4))
    
        concentration_points.append(concentrations.copy())

    return time_points, np.array(concentration_points)

Comments

[DAYANVICIEDO]

解決いたしました。本当にありがとうございます🙇

Answer 2

いくつか気になったことがあります。

while t < t_max:
    if t == 2:
        A += 1.0
    else:
    # Compute concentration changes using the fourth-order Runge-Kutta method
      k1 = reaction_equations(concentrations)
      k2 = reaction_equations(concentrations + 0.5 * dt * np.array(k1))
      k3 = reaction_equations(concentrations + 0.5 * dt * np.array(k2))
      k4 = reaction_equations(concentrations + 0.5 * dt * np.array(k3))

      concentrations += (dt / 6.0) * (2 * np.array(k1) + 2 * np.array(k2) + 2 * np.array(k3) + 2 * np.array(k4))
    t += dt
    
    time_points.append(t)
    concentration_points.append(concentrations.copy())

・Aに値を入れていますが、どこでも使っていない。
・tが2のときに、k? やconcentrations の値を更新していない。
・そもそもAを初期化していないので、このコードはエラーになると思う。