Как я могу решить уравнения рассеяния Бриллюэна с помощью Rge Kutta 4 с python? - PullRequest
Новая
       92

Как я могу решить уравнения рассеяния Бриллюэна с помощью Rge Kutta 4 с python?

1 голос
/ 19 марта 2020

На самом деле я реализую проект в 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

Сначала я собираюсь объяснить физиологию * Система 1064 *: У меня есть волокно длины L, и на одном из его концов находится электромагнитное поле c, создаваемое лазером, который колеблется во времени. Мне необходимо решить уравнения в пространстве и времени для электромагнитного поля c (El), поля рассеяния (Es), возникающего при контакте лазера с волокном, и волны давления (rho), которая вызывает указанное рассеяние.

Моя попытка решить уравнения заключается в следующем: Я предполагаю начальное распределение для электрического c поля, поля рассеяния и волны давления вдоль волокно, а затем я вычисляю пространственную производную для трех уравнений методом конечных разностей; мои начальные условия: гауссиан для поля electri c, гауссиан, умноженный на 10 ^ -10 (для имитации нуля) для поля рассеяния и константа для волны давления. Граничными условиями являются: значение нуля вне волокна для волны давления и поля рассеяния, а также для поля электромагнита c: значение нуля на конце волокна и значение поля c при начало. Я сделал и массив для каждого начального условия и для каждой производной пространства.

Это метод конечных разностей: 

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

И для запуска моего кода я использую следующий, где я делаю весь этот процесс в al oop

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 , если я печатаю решения.

Согласно моим логам c это должно сработать, и я пытался в течение нескольких недель найти, как я могу сделать эту работу, но я не знаю, что именно происходит.

Я ожидаю, что поле электричества c уменьшается, пока поле рассеяния увеличивается, и волна давления колеблется (поправьте меня, если я ошибаюсь, пожалуйста)

Пожалуйста, , если кто-то знает, как я могу решить мою проблему и ответить на вопросы, такие как ошибка из-за код? или начальное состояние? шаг моей ступеньки-кутты? или если мне нужно изменить порядок моих параметров? или это потому что комплексные числа? я буду очень очень благодарен

Спасибо за внимание.

...