Raman model equations using RK4

Question

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 = dzfunPs(Ps, Pp)
    k2 = dzfunPs(Ps+k1/2, Pp+k1/2)
    k3 = dzfunPs(Ps+k2/2, Pp+k2/2)
    k4 = dzfunPs(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)(PsNs + (hwS/2)B0)*Pp
    return a + b
def funNs(Ps, Pp, Ns, Np):
    a = (-alphaRS - alphaS + (gR/Aeff)Pp)Ns
    b = (gR/Aeff)(np.sqrt(2PpNp))Ps
    c = (alphaRS + alphaS + (gR/Aeff)Pp)(hwS/2)B0
    return a + b + c
def funNp(Ps, Pp, Np):
    a = -(alphaRP + alphaP)Np
    b = -(wS/wP)(gR/Aeff)(Ps)np.sqrt(2PpNp)
    c = (alphaRP + alphaP)(hwS/2)B0
    d = (wS/wP)(gR/Aeff)(Ps)(2Ppnp.sqrt(hwS/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 = raman_freq.readlines()
 raman_real = raman_real.readlines()
 raman_imag = raman_imag.readlines()
 #d_lamda = d_lamda.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)))
 #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 = nzdz # 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 = 2math.piFrepump
 wS = 2math.piraman_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)

Answer 1

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.

Answer 2

While writing the method you should have gotten doubts:

• k1,k2,... updates Ps
• r1,r2,... updates Pp
• q1,q2,... updates Ns
• v1,v2,... updates Np

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.