Nonlinear Schrodinger equation

[GenB31415]

I used the py-pde package (author David Zwicker) to numerically investigate the dynamics of solutions (solitonic and decaying) to the nonlinear 1D Schrödinger equation (NSE). As the base schema I started from the pde_schroedinger.py that is extended by including the non-linear part, see en.wikipedia.org/wiki/Nonlinear_Schrodinger_equation. The attached files correctly describe the expected dynamic.

[GenB31415]

In this post I show the numerical solution (with py-pde) of 1D nonlinear Schrödinger equation with Double-well potential (see //en.wikipedia.org/wiki/Double-well_potential).



. I will show the *.py files soon.

[Elliotepsteino]

Can you please post the code used to solve the nonlinear schrödinger equation?

[GenB31415]

The required code is attached. It corresponds to the soliton dynamics for NSchEq in a double-well potential V(x).
In line 32 parameter g = 5 can be variated: g=5>0 is a soliton state, g=-5<0 is a decaying state. Other parameters of the potential can be variated too.

Nonlinear Schrödinger's Equation, 1D case

This example implements a complex PDE using the :class:~pde.pdes.pde.PDE. We here
chose the Schrödinger equation <https://en.wikipedia.org/wiki/Schrödinger_equation>_
and nonlinear Schredinger equation https://en.wikipedia.org/wiki/Nonlinear_Schr%C3%B6dinger_equation],
with a spatial potential V_0.V(x) in non-dimensional form:

.. math::
i \partial_t \psi = -\nabla^2 \psi - g |\psi|^2 \psi - V_0 V(x) \psi

Note that the example imposes Neumann conditions at the wall, so the wave packet is
expected to reflect off the wall.
Cases
(i) the soliton solution(g>0, waveguide propagation), and
(ii) decaying solution (g<0) for Nonlinear Schrodinger equation.
"""

# slightly extended variant of linear (and potential free) case (g=0), see pde_schroedinger.py, by David Zwicker

from math import sqrt, cosh
import numpy as np
from pde import PDE, CartesianGrid, MemoryStorage, ScalarField, plot_kymograph
import matplotlib.pyplot as plt

grid = CartesianGrid([[0, 2*20]], 128, periodic=False)  # generate grid

# create a (normalized) wave packet with a certain form as an initial condition
center = grid.axes_bounds[0][1]/2
iVel=0.9                 # can be changed
g = 5                    # can be variated: g=5>0 is a soliton state, g=-5<0 is a decaying state
ampPot_V0 = 0.15          # V_0 is amplitude of potential V_0*V(x)

sInitState = "exp(I * "+str(iVel)+" * x) * exp(-(x - " + str(center) + ")**2)"
#sInitState = "exp(I * "+str(iVel)+" * x) / cosh(x - " + str(center)+ ")" # variant with sch
print("\n\nInit state: "+ sInitState)
initial_state = ScalarField.from_expression(grid, sInitState)
initial_state /= sqrt(initial_state.to_scalar("norm_squared").integral.real)

#V_x = f""+str(ampPot_V0)+"*((x-"+str(center)+")**2)" # V_0 * V(x) = V0*(x-c)^2
V_x = f""+str(ampPot_V0)+"*(-10*(x-"+str(center)+")**2 + 0.03*(x-"+str(center)+")**4)" # V_0 * V(x) = V0*(x-c)^2, OK!
#V_x = f"t*"+str(ampPot_V0)+"*(-10*(x-"+str(center)+")**2 + 0.03*(x-"+str(center)+")**4)" # V_0 * V(x) = V0*(x-c)^2
#----------plot potential----
def plot2D(*args):
    plt.grid(1)
    xAxis = np.arange(args[1],args[2],args[3])
    def xfunction(x,input):
        return eval(input)    
    plt.title(args[4])
    plt.plot(xAxis,xfunction(xAxis,args[0]))
    plt.xlabel("x"); plt.ylabel("V(x,t)")
    plt.show(block=False)

print("Plot potential: V(x)=" + V_x)  # V_0 * V(x)
plot2D(V_x, grid.axes_bounds[0][0], grid.axes_bounds[0][1], 0.1, "V(x)=" + V_x)
#--------------
sFun0 = f"I*(laplace(ψ)+("+str(g)+")*conjugate(ψ)*ψ**2)" # define the pde
sFun = f"{sFun0}+I*({V_x})*ψ"   # add potential V(x)
eq = PDE({"ψ": sFun})           # define the pde

# solve the pde and store intermediate data
storage = MemoryStorage()
eq.solve(initial_state, t_range=4*2, dt=1e-5, tracker=[storage.tracker(0.02)])

# visualize the results as a space-time plot
sFun = sFun.replace("conjugate(ψ)*ψ**2","ψ|ψ|^2") # simplify the math in title
sFun = sFun.replace("laplace","Δ")
sFun = sFun.replace("**","^")
sTitle = "dψ/dt=" + sFun + ", iV=%g" % (iVel)+ ", cnt=%g" % (center)
plot_kymograph(storage, scalar="norm_squared", title=sTitle)
print("Done.\n\n")

[GenB31415]

Hi Elliotepsteino The required code is attached. It corresponds to the soliton dynamics for NSchEq in a double-well potential V(x).   In line 32 parameter g = 5  can be variated: g=5>0 is a soliton state, g=-5<0 is a decaying state.Other parameters of the potential can be variated as well. Regards On Sunday, April 3, 2022, 11:18:50 PM CDT, Elliotepsteino ***@***.***> wrote: Can you please post the code used to solve the nonlinear schrödinger equation? — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: ***@***.***>

[GenB31415]

The required code is attached. It corresponds to the soliton dynamics for NSchEq in a double-well potential V(x).
In line 32 parameter g = 5 can be variated: g=5>0 is a soliton state, g=-5<0 is a decaying state. Other parameters of the potential can be variated too.

Nonlinear Schrödinger's Equation, 1D case 
======================

This example implements a complex PDE using the :class:`~pde.pdes.pde.PDE`. We here
chose the `Schrödinger equation <https://en.wikipedia.org/wiki/Schrödinger_equation>`_ 
and nonlinear Schrödinger equation <https://en.wikipedia.org/wiki/Nonlinear_Schr%C3%B6dinger_equation]>,
with a spatial potential V_0.V(x) in non-dimensional form:

.. math::
    i \partial_t \psi = -\nabla^2 \psi - g |\psi|^2 \psi - V_0 V(x) \psi
    
Note that the example imposes Neumann conditions at the wall, so the wave packet is
expected to reflect off the wall.
Cases 
(i) the soliton  solution(g>0, waveguide propagation), and
(ii) decaying solution (g<0) for Nonlinear Schrodinger equation.
"""

# slightly extended variant of linear (and potential free) case (g=0), see pde_schroedinger.py, by David Zwicker

from math import sqrt, cosh
import numpy as np
from pde import PDE, CartesianGrid, MemoryStorage, ScalarField, plot_kymograph
import matplotlib.pyplot as plt

grid = CartesianGrid([[0, 2*20]], 128, periodic=False)  # generate grid

# create a (normalized) wave packet with a certain form as an initial condition
center = grid.axes_bounds[0][1]/2
iVel=0.9                 # can be changed
g = 5                    # can be variated: g=5>0 is a soliton state, g=-5<0 is a decaying state
ampPot_V0 = 0.15          # V_0 is amplitude of potential V_0*V(x)

sInitState = "exp(I * "+str(iVel)+" * x) * exp(-(x - " + str(center) + ")**2)"
#sInitState = "exp(I * "+str(iVel)+" * x) / cosh(x - " + str(center)+ ")" # variant with sch
print("\n\nInit state: "+ sInitState)
initial_state = ScalarField.from_expression(grid, sInitState)
initial_state /= sqrt(initial_state.to_scalar("norm_squared").integral.real)

#V_x = f""+str(ampPot_V0)+"*((x-"+str(center)+")**2)" # V_0 * V(x) = V0*(x-c)^2
V_x = f""+str(ampPot_V0)+"*(-10*(x-"+str(center)+")**2 + 0.03*(x-"+str(center)+")**4)" # V_0 * V(x) = V0*(x-c)^2, OK!
#V_x = f"t*"+str(ampPot_V0)+"*(-10*(x-"+str(center)+")**2 + 0.03*(x-"+str(center)+")**4)" # V_0 * V(x) = V0*(x-c)^2
#----------plot potential----
def plot2D(*args):
    plt.grid(1)
    xAxis = np.arange(args[1],args[2],args[3])
    def xfunction(x,input):
        return eval(input)    
    plt.title(args[4])
    plt.plot(xAxis,xfunction(xAxis,args[0]))
    plt.xlabel("x"); plt.ylabel("V(x,t)")
    plt.show(block=False)

print("Plot potential: V(x)=" + V_x)  # V_0 * V(x)
plot2D(V_x, grid.axes_bounds[0][0], grid.axes_bounds[0][1], 0.1, "V(x)=" + V_x)
#--------------
sFun0 = f"I*(laplace(ψ)+("+str(g)+")*conjugate(ψ)*ψ**2)" # define the pde
sFun = f"{sFun0}+I*({V_x})*ψ"   # add potential V(x)
eq = PDE({"ψ": sFun})           # define the pde

# solve the pde and store intermediate data
storage = MemoryStorage()
eq.solve(initial_state, t_range=4*2, dt=1e-5, tracker=[storage.tracker(0.02)])

# visualize the results as a space-time plot
sFun = sFun.replace("conjugate(ψ)*ψ**2","ψ|ψ|^2") # simplify the math in title
sFun = sFun.replace("laplace","Δ")
sFun = sFun.replace("**","^")
sTitle = "dψ/dt=" + sFun + ", iV=%g" % (iVel)+ ", cnt=%g" % (center)
plot_kymograph(storage, scalar="norm_squared", title=sTitle)
print("Done.\n\n")