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(i,fai1,fai2,bc):
    return (i-np.sin(fai1)-fai2)/bc

a = 0.0
b = 100
N = 10000
h = (b-a)/N
j = 0
bc = 4.0
i = 1.5

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


while j < N-2:

    print(f'{j=}\t{fai2[j]=:.6f}\t{t[j]=:.3f}')

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

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

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

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

    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

plt.plot (t,fai2)
plt.xlabel("τ")
plt.ylabel("φ2")
plt.show()
```RCSJモデル
i = sinφ　+ dφ/dτ + βc*d^2φ/dτ^2
を、連立1階常微分方程式に変換し、ルンゲクッタ法を用いてグラフ化するコードです。本来であればサインカーブのようなグラフが出てくるはずなのですが、実際にこのコードをデバックすると、下のような全然違ったグラフが出てきてしまい、完全に行き詰まってしまいました。
どこがおかしいのか自分ではわからず、ご指導いただきたいです。どうかよろしくお願いします。
![イメージ説明](231cef9acbd37f163b3a5e5f9e26e3fb.png)

Answer 1

f2 関数の引数の順序を間違えている様です(それと、i と bc は定数なので引数に含める必要はないとは思いますが)。

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

以下が「正しい」結果かどうかは不明です。

Comments

[Ackngawe.-]

できました！
いつもありがとうございます！