

連立微分方程式を下記のようにルンゲクッタ法で解きたいのですが、実行するとEに関する複素数の項が常に0となってしまいます。  
Eの配列のデータタイプを複素数に設定してみたりしたのですが、出力は変わりません。  
なにか改良案ございましたらご教授いただきたく思います。

[code]

    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 )
    
[/code]


CC BY-SA 4.0
2022-05-09T08:19:59.307
88727
2022-05-09T12:06:21.293
2022-05-09T12:06:21.293
3060
52538
question
python numpy 数値解析
ルンゲクッタ法による連立微分方程式の解法