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 = 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 = 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)
 lamda = c/(wS/2/math.pi*1e3 + c/pumpzero)
 x = round(pump[i] - pumpzero,2)
 x = str(x)
  • 1
    $\begingroup$ 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. $\endgroup$
    – whpowell96
    Commented Mar 25 at 17:57

2 Answers 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.

  • $\begingroup$ Hi Lutz, thank you for this response, does this mean I will update only k2, r2, q2 and v2 or all of the parameters. $\endgroup$
    – Nura Adamu
    Commented Mar 26 at 3:49
  • $\begingroup$ 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) $\endgroup$
    – Nura Adamu
    Commented Mar 26 at 4:24
  • $\begingroup$ 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|})$ $\endgroup$ Commented Mar 26 at 5:21
  • $\begingroup$ 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) $\endgroup$
    – Nura Adamu
    Commented Mar 27 at 14:14
  • $\begingroup$ 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. $\endgroup$
    – Nura Adamu
    Commented 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.