連立微分方程式を下記のようにルンゲクッタ法で解きたいのですが、実行するとEに関する複素数の項が常に0となってしまいます。
Eの配列のデータタイプを複素数に設定してみたりしたのですが、出力は変わりません。
なにか改良案ございましたらご教授いただきたく思います。
import numpy as np
import matplotlib.pyplot as plt
# 時間軸の設定
T = 1*10**-9
n = 100
h = T / n
t = np.arange(0,T,h)
#定数
eta = 3.6
har = 6.626*10**-34
f0 = 375*10**12
S = (eta*eta)/(2*har*f0)
# 解きたい方程式
f = lambda E,N,t=0 : E*cp1*(((G0*N-G0*N0)/(1+ep*abs(E)*abs(E)*S))-1/tau+((C2*beta*N*N)/(S*abs(E)*abs(E))))/2
g = lambda E,N,t=0 : (1+np.sin(2*np.pi*fm*t))/e - N/(C1+C2*N+C3*N*N) - (G0*(N-N0)*(S*abs(E)*abs(E)))/(1+ep*abs(E)*abs(E)*S)
# 初期値
E_0 = complex('0+0j')
N_0 = 0
# 定数
cp1 = complex('1+5j')
G0 = 9.7*10**-13
N0 = 1.1*10**24
ep = 0.05*10**-23
tau = 2*10**-12
C1 = 2*10**8
C2 = 2*10**-16
C3 = 0
beta = 1*10**-5
fm = 1*10**9
e = 1.6*10**-19
# 結果を返す配列
E = np.empty(n,dtype="complex_")
N = np.empty(n)
E[0] = E_0
N[0] = N_0
# 反復計算
for i in range(n-1):
k_1 = h * f(E[i],N[i],t[i])
k_2 = h * f(E[i] + k_1 /2 , N[i],t[i] + h/2 )
k_3 = h * f(E[i] + k_2 /2 , N[i],t[i] + h/2 )
k_4 = h * f(E[i] + k_3 , N[i], t[i] + h )
E[i+1] = E[i] + 1/6 * (k_1 + 2*k_2 + 2*k_3 + k_4 )
j_1 = h * g(E[i],N[i],t[i])
j_2 = h * g(E[i], N[i]+ j_1 /2 ,t[i] + h/2 )
j_3 = h * g(E[i], N[i]+ j_2 /2 ,t[i] + h/2 )
j_4 = h * g(E[i], N[i]+ j_3 , t[i] + h )
E[i+1] = E[i] + 1/6 * (j_1 + 2*j_2 + 2*j_3 + j_4 )
for
ループの最後の行の先頭はE[i+1] = E[i]
ではなくN[i+1] = N[i]
なのでは? それでも未だ「解きたい方程式」の最初の行でWarningが発生しているようですが。