# Solving a large number of coupled non-linear equations

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.

• Yes, it is possible to solve a system of nonlinear equations using NDSolve. Commented Jul 4 at 17:17
• See algorithm for this kind of PDEs on mathematica.stackexchange.com/questions/265919/… Commented Jul 4 at 17:36
• @AlexTrounev Thank you for your response! Could you please elaborate how one can approach solving the non-linear equations on Mathematica numerically? Commented Jul 8 at 8:39
• If you want to solve this equation numerically, then why not determine the numerical parameters $m,\gamma, G, V0, L$? Commented Jul 8 at 13:44

## 1 Answer

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]


• Thank you for your response! I have one question: Why do you multiply the time index in plotlabel with dt (in this case, 0.01)? The differential equation is simulated for t = 100 or t=100*dt = 1? Commented Jul 10 at 9:45
• Equation is simulated for $0\le t \le 1$, dt is the time step and tmax is the number of steps. Commented Jul 10 at 11:05