0

連立微分方程式を下記のようにルンゲクッタ法で解きたいのですが、実行すると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 )
2
  • この記事 ルンゲ=クッタ法:常微分方程式をPythonで解く原理を解説 との比較だと「ルンゲ=クッタ法」というよりは「ロトカ・ヴォルテラ方程式」に見えます。そうだとするとforループの最後の行の先頭はE[i+1] = E[i]ではなくN[i+1] = N[i]なのでは? それでも未だ「解きたい方程式」の最初の行でWarningが発生しているようですが。
    – kunif
    2022年5月9日 12:10
  • コメントありがとうございます。 仰る通り最後の行はご指摘の方で合っていました。 WarningはEの初期値が0なのにEが分母にある項(0で割っている箇所)が存在しているのが問題のようですが、これは現在対処法を調べています。 2022年5月10日 5:36

0

この質問に回答するには、ログインする必要があります。

のタグが付いた他の質問を参照する。