ルンゲクッタ法を用いた回路の数値解析でforループ内が混同してしまう

Question

ルンゲクッタ法を用いて、回路に交流電流をバイアスした時の電流と電圧の特性を一々の値とともにファイルに保存できるコードを作成しています。ところが、ルンゲクッタ法によって増やしていく時間に対して、時間に依存する交流電流を０から３まで増やしていくというコードを作ろうとしたら、forループの中身があやふやになってしまいました。案の定、結果のグラフを見てみると、本来電圧が０の時に電流は０となるグラフのはずですが、なぜか電圧が０の時に電流が発生してしまうグラフになってしまいました。
（本来であれば、電流と電圧がともにステップを刻みながらお互いに原点を通って増加していくグラフになります）

どのような考え方をすれば、forループの中身が混在しないのでしょうか。
また、コードで間違っている点を、解説とともに改善策をお教えいただければ非常にありがたいです。

コード
```import numpy as np
import matplotlib.pyplot as plt

pi = np.pi

faia11 = float()
faia21 = float()
faia12 = float()
faia22 = float()
faib1 = float()
faib2 = float()
tau = float()
I = float()

A = 1
w = 0.01
i_dc = 0
I = i_dc + A * np.sin(w * tau)

def f1_1(tau,faia21):
    return faia11

def f2_1(tau,faia11,faia21):
    return (I/2 - j - np.sin(faia11) - faia21) / BC

def f1_2(tau,faia22):
    return faia21

def f2_2(tau,faia12,faia22):
    return (I/2 - j - np.sin(faia12) - faia22) / BC

def f1_3(tau,faib2):
    return faib2

def f2_3(tau,faib1,faib2):
    return (I/2 + j - np.sin(faib1)*0.5 - faib2) / BC

BC = 0.1
BL = 1

a = 0
b = 200
N = 20000
h = (b-a)/N

faia11 = float()
faia21 = float()
faia12 = float()
faia22 = float()
faib1 = float()
faib2 = float()
tau = float()
s = float()
v = float()
fai_a = 0

with open("(バイアス電流０から３まで）φa=0の交流バイアスステップ（βc=0.1,βL=1).txt", "w") as f:

    for I in np.arange(0, 3.001, 0.01):

        for k in range(0, 200000, 1):
            tau = k * h

            j = (faia11 + faia12 - faib1 - 2 * pi * fai_a)/(BL * pi)

            k1_1 = h * f1_1(tau,faia21)
            d1_1 = h * f2_1(tau,faia11,faia21)
            k1_2 = h * f1_2(tau,faia22)
            d1_2 = h * f2_2(tau,faia12,faia22)
            k1_3 = h * f1_3(tau,faib2)
            d1_3 = h * f2_3(tau,faib1,faib2)
            j = (faia11 + k1_1 + faia12 + k1_2 - (faib1 + k1_3) - 2*pi*fai_a)/(BL*pi)

            k2_1 = h * f1_1(tau + d1_1/2,faia21 + d1_1/2)
            d2_1 = h * f2_1(tau + k1_1/2,faia11 + k1_1/2, faia21 + d1_1/2)
            k2_2 = h * f1_2(tau + d1_2/2,faia22 + d1_2/2)
            d2_2 = h * f2_2(tau + k1_2/2,faia12 + d1_1/2, faia22 + d1_1/2)
            k2_3 = h * f1_3(tau + d1_3/2,faib2 + d1_3/2)
            d2_3 = h * f2_3(tau + k1_3/2,faib1 + k1_3/2, faib2 + d1_3/2)
            j = (faia11 +  k2_1 + faia12 + k2_2 - (faib1 + k2_3) - 2*pi*fai_a)/(BL*pi)

            k3_1 = h * f1_1(tau + d2_1/2,faia21 + d2_1/2)
            d3_1 = h * f2_1(tau + k2_1/2,faia11 + k2_1/2, faia21 + d2_1/2)
            k3_2 = h * f1_2(tau + d2_2/2,faia22 + d2_2/2)
            d3_2 = h * f2_2(tau + k2_2/2,faia12 + d2_1/2, faia22 + d2_1/2)
            k3_3 = h * f1_3(tau + d2_3/2,faib2 + d2_3/2)
            d3_3 = h * f2_3(tau + k2_3/2,faib1 + k2_3/2, faib2 + d2_3/2)
            j = (faia11 +  k3_1 + faia12 + k3_2 - (faib1 + k3_3) - 2*pi*fai_a)/(BL*pi)

            k4_1 = h * f1_1(tau + d3_1/2,faia21 + d3_1/2)
            d4_1 = h * f2_1(tau + k3_1/2,faia11 + k3_1/2, faia21 + d3_1/2)
            k4_2 = h * f1_2(tau + d3_2/2,faia22 + d3_2/2)
            d4_2 = h * f2_2(tau + k3_2/2,faia12 + d3_1/2, faia22 + d3_1/2)
            k4_3 = h * f1_3(tau + d3_3/2,faib2 + d3_3/2)
            d4_3 = h * f2_3(tau + k3_3/2,faib1 + k3_3/2, faib2 + d3_3/2)
            j = (faia11 +  k4_1 + faia12 + k4_2 - (faib1 + k4_3) - 2*pi*fai_a)/(BL*pi)

            faia11 += 1/6 * (k1_1 + 2 * k2_1 + 2 * k3_1 + k4_1) 
            faia21 += 1/6 * (d1_1 + 2 * d2_1 + 2 * d3_1 + d4_1) 
            faia12 += 1/6 * (k1_2 + 2 * k2_2 + 2 * k3_2 + k4_2) 
            faia22 += 1/6 * (d1_2 + 2 * d2_2 + 2 * d3_2 + d4_2)
            faib1 += 1/6 * (k1_3 + 2 * k2_3 + 2 * k3_3 + k4_3) 
            faib2 += 1/6 * (d1_3 + 2 * d2_3 + 2 * d3_3 + d4_3)

            if k > 160000:
                s += (faia21 + faia22 + faib2)/2

        v = s / 40000

        print("{:.2f} {:.10f} {:.3f} {:.10f} {:.10f} {:.10f} {:.10f} {:.10f} {:.10f} {:.10f}".format(fai_a, v, I ,j, faia11, faia21, faia12, faia22, faib1, faib2), file=f)

        s = 0

Answer 1

ChatGPTの回答（内容、動作未確認）

ルンゲクッタ法を使用して回路の特性を計算するコードを作成する際に、forループ内で混乱が生じる可能性があります。このような混乱を避けるために、以下の点に注意してください。

1. 各ループの計算に使用する変数の初期化：forループの内側で使用する変数（例：faia11、faia21、faia12など）は、各反復の開始時に適切に初期化される必要があります。これにより、前の反復の結果が次の反復に影響を与えることがなくなります。

2. 反復ごとに変数を初期化する必要がある場合：ループごとに特定の変数を初期化する必要がある場合は、ループの最初で変数を初期化してください。この場合、変数をループの前に移動させるか、ループの前に明示的な初期化ステップを追加します。

3. 変数の再代入：ルンゲクッタ法では、一時変数を使用して各ステップの計算結果を保持する必要があります。各ステップの計算結果を変数に再代入する前に、一時変数に結果を格納してください。

以下は、修正されたコードの例です。修正点はコメントで示しています。

# 他のコード

with open("(バイアス電流０から３まで）φa=0の交流バイアスステップ（βc=0.1,βL=1).txt", "w") as f:

    for I in np.arange(0, 3.001, 0.01):
        # 反復ごとに変数を初期化
        faia11 = 0.0
        faia21 = 0.0
        faia12 = 0.0
        faia22 = 0.0
        faib1 = 0.0
        faib2 = 0.0
        tau = 0.0
        s = 0.0
        v = 0.0
        fai_a = 0.0

        for k in range(0, 200000, 1):
            tau = k * h

            # 他の計算

            # k1_1などの一時変数を追加し、計算結果を格納
            k1_1 = h * f1_1(tau,faia21)
            d1_1 = h * f2_1(tau,faia11,faia21)
            k1_2 = h * f1_2(tau,faia22)
            d1_2 = h * f2_2(tau,faia12,

Comments

[Ackngawe.-]

回答ありがとうございます。
コードを回してみましたが、結果は変わらなかったです。
修正前と修正後の結果を見てみると、電圧の結果がおかしくなっているので、Iのforループは関係ないのかなと思います。
直流をバイアスしたコードの結果は正しく出ているのですが、交流をバイアスするとおかしくなってしまうみたいです、、