# Raman model equations using RK4

I am trying to solve below ODE equations for Raman model but I am having errors, mostly overflow in multiply and add. Please I need your help. Below is the code I have written so far. I am new to Python thought. Thanks

import numpy as np
import math
import matplotlib.pyplot as plt

# Solving Raman model equations using fourth order Runge Kutta methods
def RK4(Ps, Pp, Ns, Np):
k1 = dz*funPs(Ps, Pp)
k2 = dz*funPs(Ps+k1/2, Pp+k1/2)
k3 = dz*funPs(Ps+k2/2, Pp+k2/2)
k4 = dz*funPs(Ps+k3,Pp+k3)

r1 = dz*funPp(Ps, Pp, Ns)
r2 = dz*funPp(Ps+r1/2, Pp+r1/2, Ns+r1/2)
r3 = dz*funPp(Ps+r2/2, Pp+r2/2, Ns+r2/2)
r4 = dz*funPp(Ps+r3,Pp+r3, Ns+r3)

q1 = dz*funNs(Ps, Pp, Ns, Np)
q2 = dz*funNs(Ps+q1/2, Pp+q1/2, Ns+q1/2, Np+q1/2)
q3 = dz*funNs(Ps+q2/2, Pp+q2/2, Ns+q2/2, Np+q2/2)
q4 = dz*funNs(Ps+q3,Pp+q3, Ns+q3, Np+q3)

v1 = dz*funNp(Ps, Pp, Np)
v2 = dz*funNp(Ps+v1/2, Pp+v1/2, Np+v1/2)
v3 = dz*funNp(Ps+v2/2, Pp+v2/2, Np+v2/2)
v4 = dz*funNp(Ps+v3,Pp+v3, Np+v3)

Ps = Ps + (k1 + 2*k2 + 2*k3 +k4)/6
Pp = Pp + (r1 + 2*r2 + 2*r3 +r4)/6
Ns = Ns + (q1 + 2*q2 + 2*q3 +q4)/6
Np = Np + (v1 + 2*v2 + 2*v3 +v4)/6
return Ps, Pp, Ns, Np

def funPs(Ps, Pp):
a = -(alphaRS + alphaS)*Ps
b = (gR/Aeff)*Ps*Pp
return a + b

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 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 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)*(2*Pp*np.sqrt(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')
#d_lamda = open('lamda.txt')
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)))
#d_lamda = np.array(list(map(float, d_lamda)))
raman_freq = np.array(raman_freq)/(10**12)
#b_x = np.zeros(5000)
#b_y = np.zeros(5000)
d_lamda = np.zeros(5000)
#--------------------------------------------------------------------------------
gR = 0.0096 #w-1/m-1 # nonlinear coefficient
dz = 1  # m  # step size of fibre
c = 299792458  #speed of light(m/s)
nz = 5000 # 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 inial 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)

• dz=1 seems suspect. The stability of the method depends on having a small enough time step. Unless you really need a custom time-stepping routine, just use an existing library. Mar 25 at 17:57

While writing the method you should have gotten doubts:

This is realized in your code in the naming of the derivative functions and the final assemblage of the component updates.

Thus a line like

    r2 = dz*funPp(Ps+r1/2, Pp+r1/2, Ns+r1/2)


is semantically wrong, apples and oranges. It should be

    r2 = dz*funPp(Ps+k1/2, Pp+r1/2, Ns+q1/2)


Now the logistical problem is that q1 is not computed at this stage. This can be easily solve by reordering the lines of the computation, so that k1,r1,q1,v1 are computed first, then k2,r2,q2,v2 etc.

An implementation using vectors and vector-valued functions is explained in Runge-Kutta 4th order with 4 ODEs or in Energy conservation in RK4 integration scheme in C++ for a mechanical second order system. The advantage is the use of less code lines for the method and more flexibility if the state vector dimension changes.

• Hi Lutz, thank you for this response, does this mean I will update only k2, r2, q2 and v2 or all of the parameters. Mar 26 at 3:49
• Hi Lutz, I am really thankful for your help. Most of the overflow errors have been resolve based on your suggestion. But I am still getting errors from the square root that contain Pp and Np. I updated k2, r2, q2 and v2 as shown below after defining k1, r1, v1 and q1 first: k2 = dzfunPs(Ps+k1/2, Pp+r1/2) r2 = dzfunPp(Ps+r1/2, Pp+r1/2, Ns+q1/2) q2 = dzfunNs(Ps+q1/2, Pp+q1/2, Ns+q1/2, Np+v1/2) v2 = dzfunNp(Ps+v1/2, Pp+v1/2, Np+v1/2) Mar 26 at 4:24
• The associations of variable and update have to be corrected everywhere, not just in the example I gave. The square root, if it is still a problem,needs to be extended and mollified, for instance using $x/\max(\epsilon,\sqrt{|x|})$ Mar 26 at 5:21
• Hi Lutz, I did updated all the parameters as shown below:k2 = dzfunPs(Ps+k1/2, Pp+r1/2) q2 = dzfunNs(Ps+k1/2, Pp+r1/2, Ns+q1/2, Np+v1/2) r2 = dzfunPp(Ps+k1/2, Pp+r1/2, Ns+q1/2) v2 = dzfunNp(Ps+k1/2, Pp+r1/2, Np+v1/2) k3 = dzfunPs(Ps+k2/2, Pp+r2/2) q3 = dzfunNs(Ps+k2/2, Pp+r2/2, Ns+q2/2, Np+v2/2) r3 = dzfunPp(Ps+k2/2, Pp+r2/2, Ns+q2/2) v3 = dzfunNp(Ps+k2/2, Pp+r2/2, Np+v2/2) k4 = dzfunPs(Ps+k3,Pp+r3) q4 = dzfunNs(Ps+k3,Pp+r3, Ns+q3, Np+v3) r4 = dzfunPp(Ps+k3,Pp+r3, Ns+q3) v4 = dzfunNp(Ps+k3,Pp+r3, Np+v3) Mar 27 at 14:14
• I am still getting this RuntimeWarning: invalid value encountered in sqrt b = (gR/Aeff)*(np.sqrt(2*PpNp))*Ps and b = -(wS/wP)*(gR/Aeff)*(Ps)*np.sqrt(2*PpNp), which are the square root that contained Pp and Np. But I have square root that contained Pp only but there was no error. I think the problem might be from multiplying those two parameters. Please help. Mar 27 at 14:18

Don't write your own ODE solver. There are perfectly good libraries (e.g https://pypi.org/project/diffeqpy/) that use better algorithms (e.g. Tsit5), with adaptive time-stepping and error control, that will also prevent you from making typos in your ODE solver implimentation.