Как я могу решить уравнения Бриллюэна рассеяния с помощью runge kutta 4 с помощью python?
На самом деле я реализую проект на Python, и я относительно новичок в этом языке. Я хочу решить уравнения рассеяния Бриллюэна с помощью метода Рунге-Кутта 4 и провести симуляцию поведения (я прилагаю документ, в котором указаны уравнения, которые имеют комплексные значения и связаны).
https://www.usna.edu/Users/physics/mungan/_files/documents/Publications/JQE.pdf
Я использую следующие библиотеки:
import math
import cmath
import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
Сначала я собираюсь объяснить физическую систему: у меня есть волокно длиной L, и на одном из его концов находится электромагнитное поле, создаваемое лазером, который колеблется во времени. Мне нужно решить уравнения в пространстве и времени для электромагнитного поля (El), поля рассеяния (Es), создаваемого контактом лазера с волокном, и волны давления (rho), которая вызывает указанное рассеяние.
Моя попытка решить уравнения заключается в следующем: я предполагаю начальное распределение электрического поля, поля рассеяния и волны давления вдоль волокна, а затем я вычисляю пространственную производную для трех уравнений с помощью метода конечных разностей; мои начальные условия: гауссианский для электрического поля, гауссовский умноженный на 10^-10 (для имитации нуля) для поля рассеяния и константа для волны давления. Граничными условиями являются: значение нуля вне волокна для волны давления и поля рассеяния, а для электромагнитного поля: значение нуля на конце волокна и значение электрического поля в начале. Я сделал массив для каждого начального условия и для каждой производной по пространству.
Это метод конечных разностей:
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
И это мои начальные условия и мои инициализированные массивы:
#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)
Все они предназначены для определения производной по времени и решения этого метода с помощью метода Рунге-Кутты. Код, который я сделал, следующий:
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
И определенные уравнения следующие:
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
Здесь много параметров, я взял значения из документа, прикрепленного ранее, но я думаю, что они могут быть равны единице для простоты, и я сделал "f" в уравнениях документа выше равным нулю, если я ошибаюсь поправьте меня пожалуйста, я тоже не знаю, как лечить этот параметр. Используются следующие параметры:
#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
И для запуска своего кода я использую следующий, где я делаю весь этот процесс в цикле
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" - это функция, которая делает метод конечных разностей, а "rkEs" - функцию, определенную выше.
Моя проблема в том, что * когда я выполняю код, я получаю результаты, которые растут на 20 порядков на каждой итерации и растут очень быстро, создавая ошибки этого типа:
RuntimeWarning: overflow encountered in cdouble_scalars
и я получаю результаты с nan, если распечатаю решения.
По моей логике это должно сработать, и в течение нескольких недель я пытался найти, как мне это сделать, но я не знаю, что именно происходит.
Я ожидаю, что электрическое поле уменьшится, а поле рассеяния увеличится и волна давления будет колебаться (поправьте меня, если я ошибаюсь, пожалуйста)
Пожалуйста, если кто-то знает, как я могу решить свою проблему и ответить на такие вопросы, как будто ошибка связана с кодом? или начальное состояние? шаг моей рунге-кутты? или если мне нужно изменить порядок моих параметров? или это из-за комплексных чисел? я буду очень очень признателен
Спасибо за внимание.