1

实际上我正在用python实现一个项目，而且我在这种语言方面相对较新。我想用 Runge-Kutta 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

有很多参数，我取了前面附加文档的值，但我认为为简单起见，它们可以等于 1，如果我错了，我将上述文档等式中的“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获得结果。

根据我的逻辑，这应该可行，我已经尝试了几个星期，找到了如何使这项工作，但我不知道到底发生了什么。

我希望电场减小而散射场增加并且压力波振荡（如果我错了请纠正我）

请，如果有人知道我该如何解决我的问题并回答诸如错误是由于代码引起的问题？还是初始条件？我的runge-kutta的步骤？或者如果我需要更改参数的顺序？还是因为复数？我将非常非常感激

感谢您的关注。

4

0 回答 0