python ルンゲクッタ法

Question

import numpy as np
import matplotlib.pyplot as plt
from numpy.lib.type_check import isrealobj

#関数の定義
def f1(fai2):
    return fai2

def f2(fai1,fai2):
    return (i-np.sin(fai1)-fai2)/bc

#初期値
a = 0.0
b = 2.0
N = 100
h = (b-a)/N
j = 0

bc = 2.0

t = np.arange(a,b,h)
fai1_0 = 0.0
fai2_0 = 0.0

#配列
fai1 = np.empty(N)
fai2 = np.empty(N)
fai1[0] = fai1_0
fai2[0] = fai2_0

#iを増加させたとき
for i in np.arange(0,2.0,0.01):
    
    #ルンゲクッタ法
    for j in range (N-2):

        k1 = h * f1(fai2[j])
        d1 = h * f2(fai1[j] , fai2[j])

        k2 = h * f1(fai2[j] + d1/2)
        d2 = h * f2(fai1[j] + k1/2,  fai2[j] + d1/2)

        k3 = h * f1(fai2[j] + d2/2)
        d3 = h * f2(fai1[j] + k2/2 , fai2[j] +d2/2)

        k4 = h * f1(fai2[j] + d3)
        d4 = h * f2(fai1[j] + k3 , fai2[j] + d3)     

        fai1[j+1] = fai1[j] + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) 
        fai2[j+1] = fai2[j] + 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) 

        if j == N-1:
          break

        j += 1
    
        #足した値を平均化する
        v = np.average(fai2[6000:])

    plt.scatter (i, v, marker='.')

#iを減少させたとき
for i in np.arange(2.0,0,-0.01):

    #ルンゲクッタ法
    for j in range (N-2):

        k1 = h * f1(fai2[j])
        d1 = h * f2(fai1[j] , fai2[j])

        k2 = h * f1(fai2[j] + d1/2)
        d2 = h * f2(fai1[j] + k1/2,  fai2[j] + d1/2)

        k3 = h * f1(fai2[j] + d2/2)
        d3 = h * f2(fai1[j] + k2/2 , fai2[j] +d2/2)

        k4 = h * f1(fai2[j] + d3)
        d4 = h * f2(fai1[j] + k3 , fai2[j] + d3)     

        fai1[j+1] = fai1[j] + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4) 
        fai2[j+1] = fai2[j] + 1/6 * (d1 + 2 * d2 + 2 * d3 + d4) 

        if j == N-1:
          break

        j += 1
    
        #足した値を平均化する
        v = np.average(fai2[6000:])

    plt.scatter (i, v, marker='.')

plt.xlabel("v")
plt.ylabel("i")
plt.show()
```ルンゲクッタ法を用いたコードなのですが、横軸をv、縦軸をiとして,ルンゲクッタ法を用いて計算した6000以降の値を平均した値を、iを０から2.0までの間で繰り返し、それぞれ得た値をプロットしていくプログラムなのですが、実行してもがvが発散してnanとなり、グラフにプロットされません。
iをループせずに固定して、その値ひとつをプロットしたときはできたのですが、i をループさせた途端にできなくなってしまいました。
また、iの範囲は０から2.0までと指定しているのに、結果として出てくるグラフは―0.04から0.04になってしまいます。どうすればよいでしょうか。![イメージ説明](3b0313782abbebc7a36cff90e99e922b.png)

Comments

[jbpb0]

plt.scatter (i, v, marker='.')

のすぐ上に

print(i)
print(v)

を(インデントを「plt.scatter」に合わせて)追加して実行してみて、値が正しく計算されてるのかを確認してみたらいかがでしょうか

[Ackngawe.-]

確認してみたところ、iは正しくループできていましたが、vが発散してnanとなってしまっていました。
iをループせずに固定して、vの点をひとつプロットしたときはできたのですが、なぜiをループさせたらできないのでしょうか、、

[jbpb0]

> vが発散してnanとなってしまっていました。

それでしたら、

> なぜか点がグラフにプロットされません。

という質問ではなく、「v」がうまく計算できません、という質問ですね
上記分かったことを、質問を編集して追記することをお勧めします

[Ackngawe.-]

わかりました。
ありがとうございます！

[melian]

N = 100 としていますが、これですと np.average(fai2[6000:]) で範囲外アクセスになってしまいます。

[jbpb0]

> np.average(fai2[6000:]) で範囲外アクセス

は、

> iをループせずに固定して、その値ひとつをプロットしたときはできた

の場合もダメなはずですよね

[Ackngawe.-]

N ＝10000にしたらできました！
ありがとうございます！