actually i'm realizing a project in python and i'm relatively new in this language. I want to solve the brillouin scattering's equations with Runge-Kutta 4 method and make a simulation of the behavior (i attach the document where are the equations which have complex values and are accoupled).
https://www.usna.edu/Users/physics/mungan/_files/documents/Publications/JQE.pdf
The libraries that i use are the following:
import math
import cmath
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
First i'm going to explain the physic system: i have a fiber of lenght L and at one of its ends there are a electromagnetic field provided by a laser that oscillates in time. I need to solve the equations in space and time for the electromagnetic field (El), the scattering field (Es) produced by the contact of the laser with the fiber and the pressure wave (rho) which produces the said scattering.
My attempt to solve the equations is the following: i suppose a initial distribution for the electric field, the scattering field and the pressure wave along the fiber and then i calculate the space derivative for the three equations with the method of finite differences; my initial conditions are: a gaussian for the electric field, a gaussian multiply by 10^-10 (to simulate a zero) for the scattering field and a constant for the pressure wave. The boundary conditions are: value of zero outside the fiber for the pressure wave and the scattering field, and for the electromagnetic field: the value of zero at the end of the fiber and the value of the electric field at the beggining. I made and array for each initial condition and for each space derivative.
This is the method of finite differences:
def finitas(Elp,Esp):
for o in range(0,Nz-1,1):
dzEl[o]=(-3*Elp[o]+4*Elp[o+1]-Elp[o+2])/(2*dz)
dzEs[o] =(-3*Esp[o]+4*Esp[o+1]-Esp[o+2])/(2*dz)
dzEl[Nz-1]=(-3*Elp[Nz-1]+4*Elp[Nz])/(2*dz)
dzEl[Nz]=-3*Elp[Nz]/(2*dz)
dzEs[Nz - 1] = (-3 * Esp[Nz - 1] + 4 * Esp[Nz]) / (2 * dz)
dzEs[Nz] = -3 * Esp[Nz] / (2 * dz)
return 0
And these are my initial conditions and my initialized arrays:
#CONDICIONES INICIALES (Funciones en t=0)
#Campo electrico del láser
El=np.array(np.ones(Nz+1),dtype=np.complex128)
for i in range(0,Nz+1,1):
El[i] = math.exp(-(i*dz - Nz*dz/ 2) ** 2)*El0 #Gaussiana
#El[i] = El0*math.exp(-(i * dz)) *math.sin(i*dz) #Exponencial decreciente
#Campo electrico de Brillouin
Es=np.array(np.ones(Nz+1),dtype=np.complex128)
for i in range(0,Nz+1,1):
Es[i] = (math.exp(-(i * dz - Nz*dz / 2) ** 2))*Es0 #Gaussiana
#Es[i] = Es0*math.exp(-(i * dz)) *math.sin(i*dz) #Exponencial dececiente
#GRAFICAR LA PARTE REAL E IMAGINARIA POR SEPARADO
#Densidad
rho=np.array(np.ones(Nz+1),dtype=np.complex128)
for i in range(0,Nz+1,1):
rho[i] = rho0
#Diferenciales
dzEl=np.array(np.zeros(Nz+1),dtype=np.complex128)
dzEs=np.array(np.zeros(Nz+1),dtype=np.complex128)
#Nuevos
Elnuevo=np.array(np.zeros(Nz+1),dtype=np.complex128)
Esnuevo=np.array(np.zeros(Nz+1),dtype=np.complex128)
rhonuevo=Esnuevo=np.array(np.zeros(Nz+1),dtype=np.complex128)
All of these is to solve for the time derivative and solve this with the runge-kutta method. The code that i made it's the following:
def rkEs(ecEl,ecEs,ecrho):
K1 = np.array(np.zeros(3), dtype=np.complex128)
K2 = np.array(np.zeros(3), dtype=np.complex128)
K3 = np.array(np.zeros(3), dtype=np.complex128)
K4 = np.array(np.zeros(3), dtype=np.complex128)
for k in range(0,Nz,1):
x = np.array([El[k], Es[k], rho[k]], dtype=np.complex128)
Respuesta = np.array(np.zeros(3), dtype=np.complex128)
Runge=x
K1[0]=(ecEl(Runge[0],Runge[1],Runge[2],dzEl[k]))
K1[1]=(ecEs(Runge[0],Runge[1],Runge[2],dzEs[k]))
K1[2]=(ecRho(Runge[0],Runge[1],Runge[2]))
Runge[0]=x[0]+K1[0]*dt/2 #Primera evolucion de El
Runge[1]=x[1]+K1[1]*dt/2 #Primera evolucion de Es
Runge[2] = x[2] + K1[2] * dt / 2 #Primera evolucion de rho
K2[0]=(ecEl(Runge[0],Runge[1],Runge[2],dzEl[k]))
K2[1] = (ecEs(Runge[0],Runge[1],Runge[2], dzEs[k]))
K2[2] = (ecRho(Runge[0],Runge[1],Runge[2]))
Runge[0] = x[0] + K2[0]*dt/2 #Segunda evolucion de El
Runge[1] = x[1] + K2[1] * dt / 2 #Segunda evolucion de Es
Runge[2] = x[2] + K2[2] * dt / 2 # Segunda evolucion de rho
K3[0]=(ecEl(Runge[0],Runge[1],Runge[2],dzEl[k]))
K3[1] = (ecEl(Runge[0],Runge[1],Runge[2], dzEs[k]))
K3[2] = (ecRho(Runge[0],Runge[1],Runge[2]))
Runge[0] = x[0] + K3[0]*dt #Tercera evolucion de El
Runge[1] = x[1] + K3[1] * dt #Tercera evolucion de El
Runge[2] = x[2] + K3[2] * dt # Segunda evolucion de rho
K4[0] = (ecEl(Runge[0],Runge[1],Runge[2],dzEl[k]))
K4[1] = (ecEl(Runge[0],Runge[1],Runge[2], dzEs[k]))
K4[2] = (ecRho(Runge[0],Runge[1],Runge[2]))
Respuesta[0]=x[0]+(dt/6)* (K1[0]+2*K2[0]+2*K3[0]+K4[0])
Respuesta[1] = x[1] + (dt / 6) * (K1[1] + 2 * K2[1] + 2 * K3[1] + K4[1])
Respuesta[2] = x[2] + (dt / 6) * (K1[2] + 2 * K2[2] + 2 * K3[2] + K4[2])
Elnuevo[k]=Respuesta[0]
Esnuevo[k]=Respuesta[1]
rhonuevo[k]=Respuesta[2]
return 0
And the defined equations are the following:
def ecEl(Elp,Esp,rhop,dzElp):
dtEl=P3*Elp/P2 + (P/P2)*Esp*rhop*1j - P1*dzElp/P2
return dtEl
def ecEs(Elp,Esp,rhop,dzEsp):
dtEs=P3*Esp/P2 + P*Elp*rhop.conjugate()*1j/P2 + P1*dzEsp/P2
return dtEs
def ecRho(Elp,Esp,rhop):
dtRho= 1j*LAMBD*Elp*Esp.conjugate()/P1 - P4*rhop/P1
return dtRho
There is a lot of parameters, i took the values of the document atached earlier but i think that they can be equal by one for simplicity and i made the "f" in the equations of the document above equal by zero, if i am wrong please correct me, also i don't know how to treat this parameter. The used parameters are the following:
#DATOS PROBLEMA:
#Indice de refraccion
n=1.4447
#Ancho de banda(MHz)
b=20
#Ganancia del laser
z=0
#coeficiente de perdida de la fibra(db/m)
a=0.2e-3
#velocidad de la luz(m/micro s)
c=2.9970e2
#c=1
#gamma-coeficiente de electrostriccion
g=0.9002
#Densidad promedio(kg/m^3)
rho0=2210
rho0=complex(rho0)
#Longitud de onda del láser(m)
laser=1.55e-6
#Polarizacion
M=1.5
#Parametro de acomplamiento optico
P= math.pi * g / (2 * n * rho0 * laser * M)
#Permitividad electrica del vacio(A^2*micro s^4/kg*m^3)
#e0=1
e0=8.854e12
#Velocidad del sonido(m/micro s)
v=5960*1e-6
#Acoplamiento acustico
LAMBD= math.pi * n * e0 * g / (laser * v)
#Radio de la fibra(m)
r=13.75e-6
#Area modal de la fibra
A=math.pi*r**2
#Potencia inicial del laser(W-vatios)
P0=1e-3
#Longitud de la fibra(m)
L=20
#fast chirp rate(1/micro s)
beta=2e10
#CONDICIONES
t=0
tf=tf=n*L/c
#PARÁMETROS ECUACIONES
P1=1
P2=n/c
P3=(z-a)/2
P4=math.pi*b
#DIFERENCIALES Y DIMENSION DE ARREGLOS
dz=0.001
dt=n*dz/c
Nz=L/dz
Nz=int(Nz)
Nt=tf/dt
Nt=int(Nt)
#VALORES INICIALES
Es0=1*10e-10+0j
El0=math.sqrt(2*P0/(n*c*e0*A))+0j
And for running my code i use the next, where i make all this process in a loop
for i in range(0,Nt,1):
# Condicion de frontera
El[0] = El0 * (math.cos(math.pi * beta * (i * dt + n * L / c) ** 2) + 1j * math.sin(
math.pi * beta * (i * dt + n * L / c) ** 2))
Es[Nz] = 0
finitas(El, Es)
rkEs(ecEl, ecEs, ecrho)
print('El')
El=Elnuevo
Es=Esnuevo
"finitas" its the function to make the finite differences method and "rkEs" the function defined above.
My problem is* when i execute the code i obtain results that grow up in orders of 20 each iteration, and grow super super fast making appear errors of this type:
RuntimeWarning: overflow encountered in cdouble_scalars
and i obtain results with nan if i print the solutions.
According to my logic this should work and i have tried for several weeks find how can i make this work, but i don't know what is exactly happen.
I expect that the electric field decrease while the scattering field increase and the pressure wave oscilate (correct me if i'm wrong please)
Please, if someone know how can i solve my problem and answer questions like the error is due to the code? or the initial condition? the step of my runge-kutta? or if i need to change the order of my parameters? or it's because the complex numbers? i will be very very gratefully
Thanks for the attention.
