please I need help here, I am trying to implement RK4 method to solve Raman model and I am getting "RuntimeWarning: invalid value encountered in sqrt" in my code. This is mainly the sqrt that contain Pp and Np. I am not why this is happening because I have a sqrt with Pp and there was no error in that. Please help. This is my code:
import numpy as np
import math
# Solving Raman model equations using fourth order Runge Kutta methods
def RK4(Ps, Pp, Ns, Np):
k1 = dz*funPs(Ps, Pp)
q1 = dz*funNs(Ps, Pp, Ns, Np)
r1 = dz*funPp(Ps, Pp, Ns)
v1 = dz*funNp(Ps, Pp, Np)
k2 = dz*funPs(Ps+k1/2, Pp+r1/2)
q2 = dz*funNs(Ps+k1/2, Pp+r1/2, Ns+q1/2, Np+v1/2)
r2 = dz*funPp(Ps+k1/2, Pp+r1/2, Ns+q1/2)
v2 = dz*funNp(Ps+k1/2, Pp+r1/2, Np+v1/2)
k3 = dz*funPs(Ps+k2/2, Pp+r2/2)
q3 = dz*funNs(Ps+k2/2, Pp+r2/2, Ns+q2/2, Np+v2/2)
r3 = dz*funPp(Ps+k2/2, Pp+r2/2, Ns+q2/2)
v3 = dz*funNp(Ps+k2/2, Pp+r2/2, Np+v2/2)
k4 = dz*funPs(Ps+k3,Pp+r3)
q4 = dz*funNs(Ps+k3,Pp+r3, Ns+q3, Np+v3)
r4 = dz*funPp(Ps+k3,Pp+r3, Ns+q3)
v4 = dz*funNp(Ps+k3,Pp+r3, Np+v3)
Ps = Ps + (k1 + 2*k2 + 2*k3 +k4)/6
Ns = Ns + (q1 + 2*q2 + 2*q3 +q4)/6
Pp = Pp + (r1 + 2*r2 + 2*r3 +r4)/6
Np = Np + (v1 + 2*v2 + 2*v3 +v4)/6
return Ps, Ns, Pp, Np
def funPs(Ps, Pp):
a = -(alphaRS + alphaS)*Ps
b = (gR/Aeff)*Ps*Pp
return a + b
def funNs(Ps, Pp, Ns, Np):
a = (-alphaRS - alphaS + (gR/Aeff)*Pp)*Ns
b = (gR/Aeff)*(np.sqrt(2*Pp*Np))*Ps
c = (alphaRS + alphaS + (gR/Aeff)*Pp)*(h*wS/2)*B0
return a + b + c
def funPp(Ps, Pp, Ns):
a = -(alphaRP + alphaP)*Pp
b = -(wS/wP)*(gR/Aeff)*(Ps + Ns + (h*wS/2)*B0)*Pp
return a + b
def funNp(Ps, Pp, Np):
a = -(alphaRP + alphaP)*Np
b = -(wS/wP)*(gR/Aeff)*(Ps)*np.sqrt(2*Pp*Np)
c = (alphaRP + alphaP)*(h*wS/2)*B0
d = (wS/wP)*(gR/Aeff)*(Ps)*np.sqrt(2*Pp*(h*wS/2)*B0)
return a + b + c + d
# calculate Raman susceptibility
def chi_r_s(w):
m = (raman_real + 1j*raman_imag)
return m
def chi_r_i(w):
n = (raman_real + 1j*raman_imag)
return n[::-1]
pump = [1550]
for i in range(len(pump)):
raman_freq = open('Raman frequency.txt')
raman_real = open('Raman real part.txt')
raman_imag = open('Raman imag part.txt')
raman_freq = raman_freq.readlines()
raman_real = raman_real.readlines()
raman_imag = raman_imag.readlines()
raman_freq = np.array(list(map(float, raman_freq)))
raman_real = np.array(list(map(float, raman_real)))
raman_imag = np.array(list(map(float, raman_imag)))
raman_freq = np.array(raman_freq)/(10**12)
#--------------------------------------------------------------------------------
gR = 0.0096 #w-1/m-1 # nonlinear coefficient
dz = 0.01 # m # step size of fibre
c = 299792458 #speed of light(m/s)
nz = 100 # number of steps
fibre_length = nz*dz # 25 m
alphaRS = 4e-3
alphaS = 4.61e-3 # signal attenuation coefficient
alphaRP = 4e-3
alphaP = 6.9e-3 # pump attenuation coefficient
pumpzero = 1550 # nm pump wavelength
Frepump = 193.1e12 # Ghz
Freqsignal = 181.6e12
wP = 2*math.pi*Frepump
wS = 2*math.pi*raman_freq + wP #rad/ps
Aeff = 72.2e-6
h = 6.62e-34
B0 = 10e12
# Set initial values for pump, signal and noise power
Ps = 0.1 # in W
Pp = 0.3162 # in W
Ns = 1e-4 # in W
Np = 0.0032 # in W
z = 0 # propagated length
for j in range(int(nz)):
Ps, Pp, Ns, Np = RK4(Ps, Pp, Ns, Np)
z+=dz
lamda = c/(wS/2/math.pi*1e3 + c/pumpzero)
#plt.xlim(1551.4,1700)
#plt.ylim(0,11.5)
x = round(pump[i] - pumpzero,2)
x = str(x)
I tried decreasing dz and the initial values for Pp and Np.