Solving a large number of coupled non-linear equations

Question

I would like to numerically have a solution of the Gross-Pitaevskii equation for an impurity coupled with a 1D weakly interacting bosonic bath, given by:

\begin{align} i\frac{\partial \phi_0(x,t)}{\partial t} = \left\{-\frac{1 + m}{2} \partial_{x}^2 + |\phi_0(x,t)|^2 + G\delta(x) + i V(t) \partial_{x} \right\} \phi_0(x,t) \end{align} Where $V(t) = V_0 + \frac{i m}{\sqrt{\gamma}} \int dx \phi^*_0(x,t) \partial_{x} \phi_0(x,t)$ is to be solved at every time point $t$. Here, $i$ denotes the imaginary number $\sqrt{-1}$. To be able to solve this equation numerically, I discretized it using Conservative Finite Difference Scheme (CFDS) (For more details, check Eqn. (45) of https://www.mdpi.com/518624), which discretizes the original PDE into a series of non-linear coupled equations, where $n$ and $j$ denotes the discretized time and space index respectively:

\begin{align} i \frac{\phi^{n+1}_j - \phi^n_j}{\Delta t} &= -\frac{1+m}{2} \left[\frac{\phi^{n+1}_{j+1} - 2\phi^{n+1}_{j} + \phi^{n+1}_{j-1} }{2\Delta x^2} + \frac{\phi^{n}_{j+1} - 2\phi^{n}_{j} + \phi^{n}_{j-1} }{2\Delta x^2}\right] + \frac{i}{2} V^n \left[\frac{\phi^{n+1}_{j+1} - \phi^{n+1}_{j}}{\Delta x} + \frac{\phi^{n}_{j+1} - \phi^{n}_{j}}{\Delta x}\right]\nonumber\\ &+ \left[\frac{G\delta_{j}}{\Delta x} + \frac{1}{2}\left(|\phi^{n+1}_{j}|^2 + |\phi^{n}_{j}|^2 \right) \right] \frac{\phi^{n+1}_{j} + \phi^{n}_{j}}{2} \\ V^n &= V_0 + \frac{i m}{\sqrt{\gamma}} \sum\limits_{j=-L/2}^{L/2-1} \phi^{n*}_j \left(\phi_{j+1}^n - \phi_j^n \right) \end{align} In these equations, $\Delta x$ and $\Delta t$ are time steps that are to be determined by me and $V_0, m, G$, and $\gamma$ are positive parameters that I would like to provide. The initial conditions are given by $\phi(t=0,x) = \phi^0_j = 1 \ \forall j \in [-L/2, L/2)$ where $L$ is the length of the system. I would also like to implement periodic boundary condition by putting $\phi(t, -L/2) = \phi(t,L/2)$, or equivalently $\phi^n_{-L/2} = \phi^n_{L/2} \forall n \in [0, t_{\text{fin}}]$. I have two queries:

1. Can one solve the original partial differential equation numerically on Mathematica to obtain $\phi_0(x,t)$ with the given initial and boundary conditions?

2. Can the discretized non-linear coupled equations be solved in Mathematica numerically to obtain $\phi_j^n$ from the original initial and boundary conditions? If yes, what's the procedure to do so?

P.S : I don't require an analytical form of $\phi_0(x,t)$, only numerical solutions.

Answer 1

For this equation it could be better to use implicit difference scheme in combination with FEM as follows

sol[m_, gamma_, G_, V0_] := 
 Module[{L = 10, dt = 1/100, tmax = 100, eps = .05}, Psi[0][x_] := 1; 
  diracdelta[x_, eps] := Exp[-x^2/(2 eps)]/Sqrt[2  Pi  eps]; 
  V[0] = V0; 
  Do[Psi[i] = 
    Quiet@NDSolveValue[{I  (psi[x] - Psi[i - 1][x])/
          dt == -.5 (1 + m) Laplacian[
           psi[x], {x}] + (Abs[Psi[i - 1][x]]^2 + 
            G  diracdelta[x, eps]) psi[x] + V[i - 1] psi'[x], 
       PeriodicBoundaryCondition[psi[x], x == L/2, 
        Function[x, x - L]]}, psi, {x, -L/2, L/2}, 
      Method -> {"FiniteElement", InterpolationOrder -> {psi -> 2}, 
        "MeshOptions" -> {"MaxCellMeasure" -> eps}}]; 
   V[i] = V0 + 
     I  m/Sqrt[gamma] NIntegrate[
       Psi[i][x] Conjugate[Psi[i]'[x]], {x, -L/2, L/2}, 
       Method -> "LocalAdaptive", AccuracyGoal -> 5, 
       PrecisionGoal -> 4];, {i, 1, tmax}];]

Example of usage

sol[1, 1, 1, 1] // AbsoluteTiming

It takes about 5 sec on my Silver Pentium. Visualization

Plot[Evaluate[Table[Abs[Psi[i][x]], {i, 0, 100, 5}]], 
 Element[{x}, mesh], 
 PlotLegends -> Table[Row[{"t = ", i  .01}], {i, 0, 100, 5}], 
 PlotRange -> All, FrameLabel -> {"x", "|\[Psi]|"}, Frame -> True]