Python求解Nonlinear Schrödinger Equation

非线性薛定谔方程( $\alpha=0$ )：

$$
\frac{\partial A}{\partial z}=-\frac{\alpha}{2}A-i \frac{\beta_2}{2} \frac{\partial^2 A}{\partial T^2}+i\gamma|A|^2A
$$

显示了光脉冲在通过光纤传输时的包络和相位如何变化，考虑了衰减 $\alpha$ 、群速度色散 $\beta_2$ （GVD）和非线性导致的自相位调制 $\gamma$ （SPM）。

本文将展示在 Python 中使用分步傅里叶法求解非线性薛定谔方程。

查看多语言版本 English，Français。

分步傅里叶法的基本思想

光学脉冲一般用 $A$ 表示包络， $E$ 表示载波。

由于自相位调制 $\gamma$ 对脉冲的影响只在时域有解析解：

$$
A(t, z+\Delta z)=A(t, z)e^{i\gamma \Delta z|A(t, z)|^2}
$$

群速度色散 $\beta_2$ 的影响只频域有解析解：

$$
\tilde{A}(\omega, z+\Delta z)=\tilde{A}(\omega, z)e^{i\frac{\beta_2}{2}(\omega-\omega_0)^2\Delta z}
$$

因此需要先在时域计算 SPM 的影响，再通过傅里叶变换到频域，计算 GVD 和衰减的影响 ，最后通过傅里叶逆变换回时域。

基本思想可以表示为：

$$
A(t, z+\Delta z)= \underbrace{iFFT}_4 [\underbrace{e^{i\frac{\beta_2}{2}(\omega-\omega_0)^2\Delta z-\frac{\alpha}{2}\Delta z}}_3 \underbrace{FFT}_2[\underbrace{e^{i\gamma|A(T, z)|^2\Delta z} A(T, z)}_1]]
$$

结果展示（多图预警）

在我的 GitHub 仓库查看源代码。点击右上方的小月亮开启深色模式，体验更佳。

分别展示只考虑衰减、色散和自相位调制对脉冲的影响，最后展示共同作用的结果。

只考虑衰减的影响

通常光纤的衰减为 $0.2 dB/km@1550nm$。以下结果显示了输入高斯脉冲进入只有衰减效应的 1km 光纤后在时域上输出脉冲的功率变化。

以下结果显示了输入高斯脉冲进入只有衰减效应的 1km 光纤后在频域上输出脉冲的功率变化。

从输入脉冲和输出脉冲的变化可以看出，1km 的光纤只衰减了部分能量。

只考虑色散的影响

群速度色散参量 $\beta_2$ 设置为 $+100*10^3fs^2/m$。以下结果显示了输入高斯脉冲进入只有色散效应的 1km 光纤后在时域上输出脉冲的功率变化。

以下结果显示了输入高斯脉冲进入只有色散效应的1km光纤后随着传输距离的增加在时域上脉冲的强度变化。

以下结果显示了输入高斯脉冲进入只有色散效应的 1km 光纤后在时域上脉冲的啁啾变化。

色散导致脉冲在时域上展宽，脉冲前沿聚集低频光，脉冲后沿聚集高频光（时间负轴对应脉冲前沿，正轴对应脉冲后沿）。对应 $\beta_2>0$ 的负（红）啁啾（正常色散）。

只考虑自相位调制的影响

非线性系数 $\gamma$ 设置为 $10*10^{-3}rad/(W\cdot m)$ 。以下结果显示了输入高斯脉冲进入只有自相位调制效应的 1km 光纤后在频域上输出脉冲的变化。

以下结果显示了输入高斯脉冲进入只有自相位调制效应的 1km 光纤后随着传输距离的增加在频域上脉冲的频率变化。

以下结果显示了输入高斯脉冲进入只有自相位调制效应的1km光纤后随着传输距离的增加在频域上脉冲的频率变化的 3D 视图。

从以上结果可以看出，自相位调制改变了脉冲的频谱。中心频率功率减小并向两边进行展宽，分别分开成两个对称的峰，随着传输距离的增加，两个峰之间又出现两个对称的小峰值。

综合考虑衰减、色散和自相位调制

光纤的参数设置为：

• 光线长度为 $1000m$
• 衰减为 $0.2 dB/km@1550nm$
• 群速度色散参量 $\beta_2$ 为 $+100*10^3fs^2/m$
• 非线性系数 $\gamma$ 设置为 $10*10^{-3}rad/(W\cdot m)$

仿真参数设置为：

• 时域
  • 脉冲宽度 [-60, 60] ps
  • 强度截止范围 -30dB
• 频域
  • 频谱范围~[-120, 120]GHz
  • 强度截止范围 -30dB

时域的结果

频域的结果

附录

附上部分生成综合考虑衰减、色散和自相位调制图片的源代码作为参考。

plot.pySSFM.ipynb

import numpy as np
from scipy.fftpack import fft, ifft, fftshift, ifftshift, fftfreq
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import cm

global pi; pi=np.pi 

def plotFirstAndLastPulse(matrix,fiber:Fiber_config,sim:SIM_config, nrange:int,**kwargs):
  t=sim.t[int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]*1e12
  plt.figure()
  plt.title("Initial pulse and final pulse")
  plt.plot(t,getPower(matrix[0,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]),label="Initial Pulse")
  plt.plot(t,getPower(matrix[-1,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]),label="Final Pulse")
  #plt.axis([-duration*5*1e12,duration*5*1e12,0,1.05*amplitude**2])
  plt.xlabel("Time [ps]")
  plt.ylabel("Power [W]")
  plt.legend()
  saveplot('first_and_last_pulse',**kwargs)
  plt.show()  
def plotPulseMatrix2D(matrix,fiber:Fiber_config,sim:SIM_config, nrange:int, dB_cutoff,**kwargs):
  #Plot pulse evolution throughout fiber in normalized log scale
  fig, ax = plt.subplots()
  ax.set_title('Pulse Evolution (dB scale)')
  t = sim.t[int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]*1e12
  z = fiber.zlocs 
  T, Z = np.meshgrid(t, z)
  P=getPower(matrix[:,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]  )/np.max(getPower(matrix[:,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]))
  P[P<1e-100]=1e-100
  P = 10*np.log10(P)
  P[P<dB_cutoff]=dB_cutoff
  surf=ax.contourf(T, Z, P,levels=40)
  ax.set_xlabel('Time [ps]')
  ax.set_ylabel('Distance [m]')
  cbar=fig.colorbar(surf, ax=ax)
  saveplot('pulse_evo_2D',**kwargs) 
  plt.show()

def plotPulseMatrix3D(matrix,fiber:Fiber_config,sim:SIM_config, nrange:int, dB_cutoff,**kwargs):
  #Plot pulse evolution in 3D
  fig, ax = plt.subplots(1,1, figsize=(10,7),subplot_kw={"projection": "3d"})
  plt.title("Pulse Evolution (dB scale)")

  t = sim.t[int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]*1e12
  z = fiber.zlocs 
  T_surf, Z_surf = np.meshgrid(t, z)
  P_surf=getPower(matrix[:,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]  )/np.max(getPower(matrix[:,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]))
  P_surf[P_surf<1e-100]=1e-100
  P_surf = 10*np.log10(P_surf)
  P_surf[P_surf<dB_cutoff]=dB_cutoff
  # Plot the surface.
  surf = ax.plot_surface(T_surf, Z_surf, P_surf, cmap=cm.viridis,
                        linewidth=0, antialiased=False)
  ax.set_xlabel('Time [ps]')
  ax.set_ylabel('Distance [m]')
  # Add a color bar which maps values to colors.
  fig.colorbar(surf, shrink=0.5, aspect=5)
  saveplot('pulse_evo_3D',**kwargs)
  plt.show()


def plotPulseChirp2D(matrix,fiber:Fiber_config,sim:SIM_config, nrange:int,**kwargs):
  #Plot pulse evolution throughout fiber in normalized log scale
  fig, ax = plt.subplots()
  ax.set_title('Pulse Chirp Evolution')
  t = sim.t[int(sim_config.number_of_points/2-nrange):int(sim_config.number_of_points/2+nrange)]*1e12
  z = fiber.zlocs 
  T, Z = np.meshgrid(t, z)

  Cmatrix=np.ones( (len(z),len(t))  )*1.0

  for i in range(fiber.ntraces):
    Cmatrix[i,:]=getChirp(t/1e12,matrix[i,int(sim_config.number_of_points/2-nrange):int(sim_config.number_of_points/2+nrange)])/1e9


  for kw, value in kwargs.items():
    if kw.lower()=='chirpplotrange' and type(value)==tuple:
      Cmatrix[Cmatrix<value[0]]=value[0]
      Cmatrix[Cmatrix>value[1]]=value[1]


  surf=ax.contourf(T, Z, Cmatrix,levels=40,cmap='RdBu')

  ax.set_xlabel('Time [ps]')
  ax.set_ylabel('Distance [m]')
  cbar=fig.colorbar(surf, ax=ax)
  cbar.set_label('Chirp [GHz]')
  saveplot('chirp_evo_2D',**kwargs) 
  plt.show()


def plotEverythingAboutPulses(pulseMatrix,fiber:Fiber_config,
                              sim:SIM_config, 
                              nrange:int, 
                              dB_cutoff,
                              **kwargs):

  print('  ')
  plotFirstAndLastPulse(pulseMatrix,fiber,sim, nrange,**kwargs)
  plotPulseMatrix2D(pulseMatrix,fiber,sim,nrange,dB_cutoff,**kwargs)
  plotPulseChirp2D(pulseMatrix,fiber,sim,nrange,**kwargs) 
  plotPulseMatrix3D(pulseMatrix,fiber,sim,nrange,dB_cutoff,**kwargs)

  print('  ')  


def plotFirstAndLastSpectrum(matrix,fiber:Fiber_config,sim:SIM_config, nrange:int,**kwargs):
  f=sim.f[int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]/1e9
  plt.figure()
  plt.title("Initial spectrum and final spectrum")
  plt.plot(f,getPower(matrix[0,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)])*1e9,label="Initial Spectrum")
  plt.plot(f,getPower(matrix[-1,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)])*1e9,label="Final Spectrum")
  plt.xlabel("Freq. [GHz]")
  plt.ylabel("PSD [W/GHz]")
  plt.legend()
  saveplot('first_and_last_spectrum',**kwargs)
  plt.show()

def plotSpectrumMatrix2D(matrix,fiber:Fiber_config,sim:SIM_config, nrange:int, dB_cutoff,**kwargs):
  #Plot pulse evolution throughout fiber in normalized log scale
  fig, ax = plt.subplots()
  ax.set_title('Spectrum Evolution (dB scale)')
  f = sim.f[int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]/1e9 
  z = fiber.zlocs 
  F, Z = np.meshgrid(f, z)
  Pf=getPower(matrix[:,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]  )/np.max(getPower(matrix[:,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]))
  Pf[Pf<1e-100]=1e-100
  Pf = 10*np.log10(Pf)
  Pf[Pf<dB_cutoff]=dB_cutoff
  surf=ax.contourf(F, Z, Pf,levels=40)
  ax.set_xlabel('Freq. [GHz]')
  ax.set_ylabel('Distance [m]')
  cbar=fig.colorbar(surf, ax=ax) 
  saveplot('spectrum_evo_2D',**kwargs) 
  plt.show()

def plotSpectrumMatrix3D(matrix,fiber:Fiber_config,sim:SIM_config, nrange:int, dB_cutoff,**kwargs):
  #Plot pulse evolution in 3D
  fig, ax = plt.subplots(1,1, figsize=(10,7),subplot_kw={"projection": "3d"})
  plt.title("Spectrum Evolution (dB scale)")

  f = sim.f[int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]/1e9 
  z = fiber.zlocs 
  F_surf, Z_surf = np.meshgrid(f, z)
  P_surf=getPower(matrix[:,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]  )/np.max(getPower(matrix[:,int(sim.number_of_points/2-nrange):int(sim.number_of_points/2+nrange)]))
  P_surf[P_surf<1e-100]=1e-100
  P_surf = 10*np.log10(P_surf)
  P_surf[P_surf<dB_cutoff]=dB_cutoff
  # Plot the surface.
  surf = ax.plot_surface(F_surf, Z_surf, P_surf, cmap=cm.viridis,
                        linewidth=0, antialiased=False)
  ax.set_xlabel('Freq. [GHz]')
  ax.set_ylabel('Distance [m]')
  # Add a color bar which maps values to colors.
  fig.colorbar(surf, shrink=0.5, aspect=5)
  saveplot('spectrum_evo_3D',**kwargs) 
  plt.show()


def plotEverythingAboutSpectra(spectrumMatrix, # 光谱矩阵
                               fiber:Fiber_config, # 光纤属性
                               sim:SIM_config,  # 仿真参数
                               nrange:int, # 画图范围
                               dB_cutoff, # 画图范围
                               **kwargs):

  print('  ')  
  plotFirstAndLastSpectrum(spectrumMatrix,fiber,sim, nrange,**kwargs)
  plotSpectrumMatrix2D(spectrumMatrix,fiber,sim,nrange,dB_cutoff,**kwargs)
  plotSpectrumMatrix3D(spectrumMatrix,fiber,sim,nrange,dB_cutoff,**kwargs)

  print('  ')  
#Plot pulses
nrange=600 # 时域的范围 -60，60ps
cutoff=-30 # 时域强度截止值dB

plotEverythingAboutPulses(pulseMatrix,fiber,sim_config,nrange,cutoff,savename='gaussian',chirpPlotRange=(-40,40))

#Plot spectra
nrange=400 # 频域的范围
cutoff=-30 # 频域强度截止值dB
plotEverythingAboutSpectra(spectrumMatrix,fiber,sim_config,nrange,cutoff,savename='gaussian')