I am facing stack overflow error when I am trying to solve numerically the fractional nonlinear Schrödinger equation using finite difference method in python.

This nonlinear partial differential equation (PDE) models how the envelope and phase of light pulse changes when propagating through a single mode optical fiber, when taking power attenuation ($alpha$), group velocity dispersion ($beta_2$)and waveguide nonlinearity ($gamma$) causing self-phase modulation (SPM) into account.

The complete error code can be found here:

c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:82: RuntimeWarning: overflow encountered in multiply
  return - 1j * fiber.gamma * getPower(pulseVector) * pulseVector
c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:108: RuntimeWarning: invalid value encountered in multiply
  K2 = (fiber.dz)**alpha * RightHandSide(fiber,sim, pulseVector + a11 * K1)
c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:73: RuntimeWarning: overflow encountered in divide
  return ddy / (sim.time_step ** 2)
c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:79: RuntimeWarning: invalid value encountered in multiply
  return 1j * fiber.beta2 / 2 * SecondTimeDerivative(sim,pulseVector)
c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:24: RuntimeWarning: overflow encountered in square
  return np.abs(amplitude) ** 2
c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:82: RuntimeWarning: invalid value encountered in multiply
  return - 1j * fiber.gamma * getPower(pulseVector) * pulseVector
c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:109: RuntimeWarning: invalid value encountered in multiply
  K3 = (fiber.dz)**alpha * RightHandSide(fiber,sim, pulseVector + a22 * K2 + a21 * K1)
c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:107: RuntimeWarning: invalid value encountered in multiply
  K1 = (fiber.dz)**alpha * RightHandSide(fiber,sim, pulseVector)
c:UsersBalazsDesktopkutatásfractional nonlinear Schrödinger equationsimulation_v1.7functions.py:111: RuntimeWarning: invalid value encountered in multiply
  return pulseVector + w1 * K1 + w2 * K2 + w3 * K3

My full code can be found here.

The another problem is that, I cannot set the plotting functions, since I am unable getting the solution.

Before using the finite difference method, I have successfully solved numerically this PDE using Split-Step Fourier method. I used SI units. The time step was 0.01 femtosecond, the spatial step was around 50 micron. Unfortunately I can not set the fractional integrator (or any other kind of integrator), when I am using this approach. The next step was to find another solution for my problem. Namely the finite difference method.

My naive (and first) approach was to solve numerically in the easiest way. I faced stack overflow error for the first time. I used the Euler and the fourth order Runge-Kutta integrators. I used SI units. The time step was 0.01 femtosecond, the spatial step was around 1 nm. I did not use adaptive time step. I saved the solution matrix to csv to inspect the raw data. I noticed that the numbers were growing monotonically very fast till I got inf an NAN values. This means that the solution was exploded in finite time. This happened probably because the time and/or the spatial step was too small.

My second approach was to scale the quantities into atomic units using finite difference method. I could change the units between SI and atomic. I did not use adaptive time step. Well, I got the solution in some cases, but generally it did not solve the situation.

My third approach was to solve the PDE using pseudo spectral Fourier method. I could change the units between SI and atomic. I did not use adaptive time step. Well, I got the solution in some cases, but generally it did not solve the situation.

My fourth approach was to return using the finite difference method. I realized that I can reshape the whole PDE into dimensionless form. I could change the units between SI and dimensionless. I did not use adaptive time step. I also implemented the 3-stage fractional order Runge-Kutta integrator. With the other two integrators, now the PDE can be solved without stack overflow error. However the fractional-order ingegrator I still get stack overflow error no matter how I refine the time and/or the space grid. And also the results looks strange to me, compared to the expected results.

The expected results can be found here.

I’m completely out of ideas right now. Any new approaches are welcomed.