Solving Nonlinear Schrödinger Equation in Python

L’équation de Schrödinger non linéaire( $\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
$$

montre comment l’enveloppe et la phase d’une impulsion optique changent lorsqu’elle se propage à travers une fibre optique, en tenant compte de l’atténuation $\alpha$, de la dispersion de la vitesse de groupe (GVD) $\beta_2$ et de la modulation de phase auto-induite (SPM) $\gamma$ causée par la non-linéarité.

Cet article démontrera l’utilisation de la méthode de Fourier à pas divisé (Split-Step Fourier method) pour résoudre l’équation de Schrödinger non linéaire en utilisant Python.

Cet article est également disponible en English, 简体中文.

L’idée fondamentale de la Split-Step Fourier method

Les impulsions optiques sont généralement représentées par l’enveloppe $A$ et l’onde porteuse $E$.

Parce que l’effet de la modulation de phase auto-induite $\gamma$ sur l’impulsion n’a qu’une solution analytique dans le domaine temporel.

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

L’effet de la dispersion de la vitesse de groupe $\beta_2$ sur l’impulsion n’a qu’une solution analytique dans le domaine fréquentiel:

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

Il est donc nécessaire de calculer d’abord l’effet de la modulation de phase auto-induite (SPM) dans le domaine temporel, puis d’effectuer une transformation de Fourier vers le domaine fréquentiel pour calculer l’effet de la dispersion de la vitesse de groupe (GVD) et de l’atténuation, et enfin d’utiliser la transformée de Fourier inverse pour revenir au domaine temporel.

L’idée fondamentale peut être exprimée comme suit :

$$
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]]
$$

Affichage des résultats (avertissement de plusieurs figures)

Le code source peut être trouvé dans mon dépôt GitHub. Cliquez sur l’icône de la lune dans le coin supérieur droit pour activer le mode sombre afin d’améliorer votre expérience.

Les effets de l’atténuation, de la dispersion et de la modulation de phase auto-induite sur l’impulsion sont affichés séparément, suivis des résultats de leur action combinée.

Effets de l’atténuation uniquement

L’atténuation typique des fibres optiques est de 0.2 dB/km@1550 nm. Les résultats suivants montrent la variation de puissance de l’impulsion de sortie dans le domaine temporel lorsque l’impulsion gaussienne d’entrée pénètre dans une fibre de 1 km en tenant compte uniquement des effets d’atténuation.

Les résultats suivants montrent la variation de puissance de l’impulsion de sortie dans le domaine fréquentiel lorsque l’impulsion gaussienne d’entrée pénètre dans une fibre de 1 km en tenant compte uniquement des effets d’atténuation.

À partir des changements observés dans les impulsions d’entrée et de sortie, on peut constater que seule une partie de l’énergie est atténuée dans la fibre de 1 km.

Effets de la dispersion uniquement

Le paramètre de dispersion de la vitesse de groupe $\beta_2$ est réglé à $+100 \times 10^3$ fs²/m. Les résultats suivants montrent la variation de puissance de l’impulsion de sortie dans le domaine temporel lorsque l’impulsion gaussienne d’entrée pénètre dans une fibre de 1 km en tenant compte uniquement des effets de dispersion.

Les résultats suivants montrent la variation de l’intensité de l’impulsion dans le domaine temporel à mesure que la distance de transmission augmente après l’entrée d’une impulsion gaussienne dans une fibre optique de 1 km en tenant compte uniquement de l’effet de dispersion chromatique.

Les résultats suivants montrent la variation de la chirp de l’impulsion de sortie dans le domaine temporel lorsque l’impulsion gaussienne d’entrée pénètre dans une fibre de 1 km en tenant compte uniquement des effets de dispersion.

La dispersion provoque un élargissement de l’impulsion dans le domaine temporel, avec la lumière de basse fréquence s’accumulant sur le front d’attaque de l’impulsion et la lumière de haute fréquence s’accumulant sur le bord de fuite (l’axe du temps négatif correspond au front d’attaque de l’impulsion et l’axe positif correspond au bord de fuite de l’impulsion). Cela correspond à un chirp négatif (rouge) ($\beta_2 > 0$) (dispersion normale).

Effets de la modulation de phase auto-induite uniquement

Le coefficient non linéaire $\gamma$ est réglé à $10 \times 10^{-3} , \text{rad/(W} \cdot \text{m)}$. Les résultats suivants montrent la variation de l’impulsion de sortie dans le domaine fréquentiel lorsque l’impulsion gaussienne d’entrée pénètre dans une fibre de 1 km en tenant compte uniquement des effets de modulation de phase auto-induite.

Les résultats suivants montrent la variation de fréquence de l’impulsion dans le domaine fréquentiel lorsque l’impulsion gaussienne d’entrée pénètre dans une fibre optique de 1 km en tenant compte uniquement des effets de modulation de phase auto-induite, avec l’augmentation de la distance de transmission.

Les résultats suivants montrent une vue 3D de la variation de la fréquence de l’impulsion dans le domaine fréquentiel avec l’augmentation de la distance de transmission, lorsque l’impulsion gaussienne d’entrée pénètre dans une fibre optique de 1 km en tenant compte uniquement des effets de modulation de phase auto-induite.

Les résultats ci-dessus montrent que la modulation de phase auto-induite modifie le spectre de l’impulsion. La puissance de la fréquence centrale diminue et se répartit de manière élargie de chaque côté, se séparant en deux pics symétriques. Avec l’augmentation de la distance de transmission, deux petites valeurs de crête symétriques apparaissent entre les deux pics principaux.

En prenant en compte simultanément l’atténuation, la dispersion et la modulation de phase auto-induite

Les paramètres de la fibre optique sont définis comme suit :

• La longueur de la fibre optique est de $1000m$
• L’atténuation est de $0.2 dB/km@1550nm$
• Le paramètre de dispersion de vitesse de groupe $\beta_2$ 为 $+100*10^3fs^2/m$
• Le coefficient de non linéarité $\gamma$ est réglé à $10*10^{-3}rad/(W\cdot m)$

Les paramètres de simulation sont configurés comme suit :

• Dans le domaine temporel.
  • La largeur de l’impulsion [-60, 60]ps
  • La plage de coupure d’intensité -30dB
• Dans le domaine fréquentiel
  • La plage spectrale ~[-120, 120]GHz
  • La plage de coupure d’intensité -30dB

Les résultats dans le domaine temporel

Les résultats dans le domaine fréquentiel

Annexe

Veuillez trouver ci-joint un extrait du code source générant les images pour prendre en compte simultanément l’atténuation, la dispersion et la modulation de phase auto-induite : lien pour référence.

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')